A correlation propagation model for nonlinear fourier transform of second order solitons

Inverse scattering transform or nonlinear Fourier transform (NFT) has been proposed for optic communication to increase channel capacity beyond the well known Shannon limit. Within NFT, solitons, as discrete outputs of the transform, can be a type of resource to carry information. Second-order solitons as the most basic higher order solitons show correlations among their parameters in the nonlinear Fourier domain as they propagate along a fibre. In this work, we report, for the first time, a correlation propagation model for second-order soliton pulses in the nonlinear Fourier domain. The model can predict covariance matrices of soliton pulses at any propagation distance using only the covariance matrices calculated at the input of the fibre with different phases in the nonlinear Fourier domain without the need of propagating the pulses.

www.nature.com/scientificreports/ a constant signal power throughout the length of the fibre. The study of correlation properties when there is significant power loss in the signal is not in the scope of this work. The main body of this report is structured as following: we first introduce a few fundamental equations of NFT that are relevant to our discussion, then we present our correlation propagation model and discuss the parameter space of second-order soliton pluses. Finally we conclude our work.
Nonlinear Fourier transform. The evolution of a noisy signal inside a lossless optical fibre in the nonlinear regime can often be expressed using the following stochastic nonlinear Schrödinger equation where κ is the noise strength, and N(t, z) is a Gaussian white noise term that follows the rule When noise is zero ( κ = 0 ), Eq. (1) can be solved analytically using the inverse scattering method 19 , also called Nonlinear Fourier Transform (NFT), where a corresponding scattering problem associated with Eq. (1) can be found as with the boundary condition Two scattering coefficients can be defined as The NFT is defined as a transformation from q(t, z) → Q c ( , z) and Q d ( k , z) , where the continuous part of the spectrum Q c ( , z) = b( ) a( ) for real and the discrete part of the spectrum Q d ( k , z) = b( ) a ′ ( ) for discrete points of k in the upper complex plane that is corresponding to the roots of a( ) . In the nonlinear Fourier domain, the propagation of the nonlinear spectra of a pulse can be expressed as: where |Q c ( , 0)| and |Q d ( k , 0)| are the initial NFT spectral magnitude and ∠Q c ( , 0) and ∠Q d ( k , 0) are the initial phase at the beginning of the fibre. In this work, we focus on the discrete part of the spectrum and use subscripts |Q d k | and ∠Q d k to represent the discrete spectral magnitude and phase associated with the discrete eigenvalue k . Using Eq. (7b), one can find the NFT spectral magnitude and phase at position z as, Considering a pulse with two eigenvalues 1 and 2 propagating along a fibre, due to the differences in 1 and 2 , ∠Q d 1 and ∠Q d 2 increase at a different rate as the propagation distance increase. We define the difference between the two phases as the "NFT phase difference" , where Previously, we discovered that the correlation between eigenvalues in multi-eigenvalue systems is influenced by the NFT phase difference due to propagation, i.e., 4zRe( 2 1 − 2 2 ) 16 . In general, Eqs. (8) indicate that in the NFT domain, propagating a pulse noiselessly and losslessly is equivalent to scale its spectral magnitude and change its initial phase.
Correlation propagation model. Optical solitons are pulses with only discrete eigenvalues. In this paper, we investigate the "correlation" between eigenvalues and the spectral magnitudes or phases associated with the eigenvalues for two-eigenvalue solitons (also as known as second order solitons).
Previous work 16 shows that the correlation of the imaginary part of the eigenvalues of a collection of 2sech(t) soliton pulses change periodically as the pulse propagates along a fibre. In general, both eigenvalues and , i is the variance of i and σ i is the standard deviation of i. It is worth mentioning that the covariance matrix is defined with respect to a given q(t) with the addition of a white Gaussian noise.
To quantify the correlation, we consider the correlation coefficient and to simplify our notation, a two-eigenvalue pulse will be denoted in time and NFT domain, respectively, by Propagation and noise addition. Previous work 16 suggests that there are different types of perturbations from different sources, .i.e, transmitter, amplifier or receiver. The perturbation model for eigenvalues, which was developed in the previous work 16 , led to two main conclusions: (1) the overall perturbation in eigenvalues can be approximated as the linear addition of all the perturbation from different sources and (2) the distributed noise along a fibre can be approximated by a set of point noise. Here, we review the main points of the model and expand it to include spectral magnitudes and phases. A noisy fibre of length L is divided into M equally separated segments as shown in Fig. 1. Within each segment, the noise distributed throughout the segment is approximated by a single point noise n m (t) added at the end of the segment, while the segment itself becomes noiseless. A soliton pulse q m (t) propagates through the segment m of the fibre noiselessly. The perturbation in eigenvalues and Q d k within the segment is then evaluated using the pulse q ′ m (t) = q m (t) + n m (t). The influence of noise of the previous segment has negligible effects on the current segment 16 . The perturbations between segments are considered independently. Hence, the total perturbation in and Q d k at the end of the fibre can be approximated as the sum of the perturbation in each segment. Any other point noise sources, i.e. transmitter noise and receiver noise, can be added to the corresponding segments where the noise sources sit. As a consequence, the correlation matrix of the 8 degrees of freedom of a two-soliton pulse at the end of the propagation can be written as the summation of the correlation matrices that only have noise in one of the segments along the fibre. In other words, where the subscript m denotes the segment in which the noise is added.
Since the eigenvalues remain unchanged during noiseless and lossless propagation, km (L) can be replaced by km (z m ) where z m is the distance from the input of the fibre to the end of segment m. For spectral amplitude Q d k , its changed is determined by Eq. (7b) for noiseless and lossless propagation. Hence, Q d km (L) can be replaced by Q d km (z m )e −4j 2 km (z m )l m , where l m = L − z m . Eq. (12) can be rewritten as Figure 1. Noise addition model. Fibre is divided into "M" segments. In each segment, the noise is modelled by approximating it with a noiseless propagation and point noise added at the end of the propagation. The perturbations of eigenvalues and spectra between segments are assumed to be independent. The total perturbation at the end of the fibre is approximated as a linear addition of the perturbations in each segment plus a compensations of the spectral phase based on the distances from each segment to the end of the fibre.  (13) is inconvenient to computer. In order to calculate the covariance matrices on the RHS of the equation, one has to save all the instances of noisy k and Q d k in all the segments and update them every time L changes. To simplify the model, the following approximation is applied to write the covariance matrices in the RHS of Eq. (13), which is evaluate at L, in terms of the covariance matrices evaluated at each segment ( L = z m ).
By separating the real and imaginary part of k as k = Re( k ) + jIm( k ) , we find Assume, the noiseless value for km is k0 , then we can write km as Replace km in Eq. (14) using Eq. (15) and assume km is sufficiently small, we get the following equations for the magnitude of the spectrum, for the phase, no approximation is needed, hence We , and the equations above can be rewritten as, Inserting Eqs. (18) into (13) and applying the covariance linear combination rules, each of the covariance matrices on the RHS of Eq. (13) can now be written as (subscript m is ignored for notation simplicity): where ⊙ is the Hadamard product such that: [A ⊙ B] ij = A ij B ij , and The significance of Eq. (19) is that instead of saving all the instances of noisy k and Q d k , and re-evaluate them for every L, one only need to evaluate and save the three matrices on the RHS of Eq. (19) once for every segment. These matrices are independent of propagation length L. Inserting Eq. (19) into Eq. (13), we finally get www.nature.com/scientificreports/ Since the signal propagates to each segment noiselessly, we can simply replace the propagation with adding a scaling factor and an extra phase to Q d k according to Eqs. (8). To use the model, one needs to first calculate the covariance matrices for a few different equally spaced within the range of [0,2π ). The optimum number of points can be determined by gradually increasing the number of points and checking for convergence. One also need to keep in mind that km in all the segments need to satisfy the condition 8Re( k0 )Im(� km )l m ≪ 1 or 8Im( k0 )Re(� km )l m ≪ 1 in order for Eq. (16) to be valid. In this report, ten points is used. Once the covariance matrices for different are obtained, map all the segments to corresponding and scale the terms in the matrices associated with |Q d k | by the factor e 4Im( 2 k )z m , such as where x stands for any one of the parameters of the covariance matrices. Finally, use Eq. (24) to predict the covariance matrices at desired propagation distance. With Eq. (24), one can predict the covariance matrix at any propagation distance by using a set of precalculated covariance matrices at the beginning of a fibre with different phases in Q d k without the need of doing statistical pulse propagation calculations at all.
Take a pulse with 1,2 = (6j, 5j) and Q d 1,2 = (−6j, −2j) as an example and consider a propagation distance L 2π corresponding to 2π phase difference ( L 2π = 2π/ 4Re( 2 1 − 2 2 ) ≈ 0.143 in this example), the corresponding correlation coefficient ρ a,b for all the top right half of the covariance matrix (symmetric with respect to the main diagonal axis) from L = 0 to 4L 2π are plotted in Fig. 2. The values predicted using the model are plotted as solid curves and ones simulated through nonlinear pulse propagation are plotted as dashed curves. A Splitstep Fourier method is applied for the pulse propagation simulation 20 and the rest of the simulation details can be found in the section "Numerical settings". As one can see, the model can produce reasonably good predictions of the covariance matrix at any propagation distance. Note that the length is a normalised unitless length. It can be scale to physical units according to proper fibre and signal parameters 8 . For instance, for a SMF-28 fibre, β 2 = −20 ps 2 /km , γ = 2 −1 km −1 at 1.55 μμm and a pulse q(t) = Asech(t/T0) with 100 ps pulse duration ( T 0 = 100 ps ), the L 2π length corresponds to approximately 143km. Note that the propagation is assumed to be lossless, i.e. a Raman amplified is used to compensate the loss. The model should not be applied to the cases with point amplification where erbium doped fibre amplifiers are used due to the continuous growth of the continuous spectrum.
Another example, Fig. 3 shows a case where the real part of the second eigenvalue is non-zero ( 1,2 = (6j, 1 + 5j) ). In this case, the signal separates apart into two pulses as the it propagates through the fibre. The results in Fig. 3 show good approximated predictions can still be made using the model. Parameter space. So far, we have described a model which allows us to find the propagation of correlations among eigenvalues, spectral magnitudes and phases at any propagation distance, by using pre-calculated covariance matrices at the beginning of the fibre. Those covariance matrices have eight degrees of freedom (two for each eigenvalue and spectrum). However, it is important to note the following properties of NFT 8 , due to which the number of degrees of freedom can be reduced.
Suppose q 1 (t) is obtained from q(t) through one of the above four transformations. It turns out that the covariance matrix for q 1 (t) can also be obtained directly from the covariance matrix for q(t) as well. To illustrate the idea, suppose that q 1 (t) = q(t − t o ) . Consider the ensemble of random signals q(t) + n(t) where n(t) is the additive white Gaussian noise added to q(t). Due to that n(t) is stochastic, the discrete eigenvalues and the spectral magnitudes of q(t) + n(t) (denoted by ( Q d k , k , k = 1, 2) ) will be stochastic as well. Using the transformation properties, the discrete eigenvalues and the spectral magnitudes of q(t − t o ) + n(t − t o ) will be distributed as (e −2j k t 0 Q d k , k , k = 1, 2) . As the white noise is "time-shift" invariant in the sense that stochastic properties of n(t) and n(t − t o ) are equivalent. The discrete eigenvalues and the spectral magnitudes of q(t − t o ) + n(t) will also be distributed as (e −2j k t 0 Q d k , k , k = 1, 2) . For that reason, we can determine the covariance matrix for q(t − t o ) + n(t) from that for q(t) + n(t).
Following the discussion above, it can be concluded that, when one needs to map out the covariance matrices of the whole parameter space while designing a two-eigenvalue system, only four degrees of freedom need to be considered. They are 1 , |Q d 1 |/|Q d 2 | and the NFT phase difference = ∠Q d 1 − ∠Q d 2 ( 2 as the fourth degree if the frequency is different). Once the covariance matrices along these degrees of freedom have been obtained, the covariance matrices along all other degrees of freedom can then be computed through the transformation mentioned above.
Numerical settings. This work relies heavily on numerical simulations. The time signal of the soliton pulses for the given Q d k and k are calculated using the Darboux transformation method 21 with a time window of [−10, 10] and a resolution of 2 14 points across the time window. The time window surround the pulse is chosen such that pulse magnitudes at the edges of the time windows are smaller than -70 dB with respect to the pulse peak. www.nature.com/scientificreports/ A split-step Fourier method 20 is used for pulse propagation simulation. A step size of 10 −5 is chosen for the pulse propagation simulation results in Figs. 2 and 3. In these results, White Gaussian noise where added into the simulation as stochastic noise. In Figs. 2 and 3, there are 40 points along the x-axis, each point corresponding to a segment of length 0.1L 2π or 1428 split steps.
To achieve high numerical accuracy, we have implemented a bi-redirection NFT algorithm 22 using a fourthorder Runge-Kutta method. With 2 14 sampling points, we can obtain accuracy more than 10 significant figures for eigenvalues and 8 significant figures for spectral magnitudes.
To obtain the covariance matrices, 5000 runs were involved in the Monte-Carlo simulation. This applies to both the pulse propagation simulation and the covariance matrix calculation at the beginning of the fibre. The Noise strength of κ = 0.01 is used in all simulations.

Conclusion
In this work, we have developed a model to estimate the covariance matrix of all 8 degrees of freedom of secondorder soliton pulses at any propagation distances based on the parameters at the beginning of the fibre. The model assumes ideal Raman amplification is provided to counteract the effect of loss in the fibre. We also discussed how one can map out the covariance matrix of the whole parameter space by using only 4 degrees of freedom. The model can be used to help design communication systems such as predicting minimal distances between clusters in the constellation diagrams as well as optimising the decision boundaries for detection. Furthermore, since NFT is the same as inverse scattering transform, the model may provide useful insights into other nonlinear physics research areas such as studying ocean rouge waves 23 .