Signal-noise interaction in nonlinear optical fibers: a hydrodynamic approach

We present a new perturbative approach to the study of signal-noise interactions in nonlinear optical fibers. The approach is based on the hydrodynamic formulation of the nonlinear Schr\"odinger equation that governs the propagation of light in the fiber. Our method is discussed in general and is developed in more details for some special cases, namely the small-dispersion regime, the continuous-wave (CW) signal and the solitonic pulse. The accuracy of the approach is numerically tested in the CW case.


Introduction
In fiber-optic communication systems, optical amplifiers are usually employed to periodically compensate for fiber loss. As a side effect of amplification, amplified spontaneous emission (ASE) noise is added to the optical signal at each amplification stage. Both signal and noise propagate through the optical fiber according to the nonlinear Schrödinger equation (NLSE), which accounts for attenuation, dispersion, and Kerr nonlinearity. In the linear regime (i.e., at low power), the ASE noise is not affected by propagation through the fiber and can be modeled as additive white Gaussian noise (AWGN) at the end of the optical link. On the other hand, when nonlinear effects are not negligible, signal and noise interact during propagation. This interaction manifests, for instance, in the parametric gain and modulation instability effects [1] and in the emergence of nonlinear phase noise [2]. As a result, noise becomes, in general, colored, non-stationary, and non-Gaussian. The knowledge of the statistics of the noisy signal at the output of the link is extremely important to design efficient optical systems, evaluate their performance, and establish the ultimate channel capacity. For instance, optimum detection (in a maximum likelihood or maximum a-posteriori probability sense) is based on the knowledge of the conditional distribution of the output signal given the input [3]. Unfortunately, this distribution is, in general, unknown. Exact statistical models in the absence of fiber dispersion have been studied in [4][5][6]. In the presence of dispersion, approximate models based on perturbation methods have been investigated in [7][8][9]. Some alternative (non-regular) perturbation methods have been also proposed to describe the noise as a non-additive perturbation of the noise-free solution. For instance, the combined regular-logarithmic perturbation (CRLP) [10] provides an accurate description of the non-Gaussian distribution through a non-linear combination of Gaussian variables, accounting for both parametric gain and nonlinear phase noise. Moreover, in [3] it is shown that a polar-Gaussian metric may provide a better performance than a Cartesian-Gaussian metric for maximum likelihood sequence detection. This suggests that a perturbation approach applied to a polar representation (e.g., amplitude and phase) of the signal may be a good alternative to obtain a simple but effective statistical model.
A similar problem arises in the different context of high-power fiber lasers, in which an incoherent or partially-coherent CW light is subject to spectral broadening during propagation through an optical fiber [11][12][13][14]. Also in this context, the important role of the interaction between fiber dispersion and nonlinearity and the suitability of a polar representation of the signal have been widely recognized.
In this paper we develop a novel approach to the problems outlined above, which is based on the hydrodynamic formulation of the NLSE. Such formulation, originally introduced by E. Madelung in a quantum mechanical framework [15] and subsequently applied also to optical fibers (e.g, in [16][17][18][19]), consists of a system of compressible, isothermal, Euler-like equations with a typical "quantum pressure" term, also known as Bohm potential. In such system the hydrodynamic variables have the physical meaning of power and angular frequency (playing the role of the fluid density and velocity, respectively). When the dispersion is small compared to the nonlinearity, the quantum pressure term can be formally neglected and we are left with a much simpler hydrodynamic system which is similar in form to the so called "shallow-water equations". Such a semiclassical approximation has been employed, for instance, in the study of self-similar propagation of high-power pulses in optical fiber amplifiers [20-23]. It must be noticed, however, that in this way the dynamics loses its dispersive character, which may lead to formation of shocks or even to ill-posedness of the mathematical problem [24]. Both the full and the approximated hydrodynamic equations are then used to study the propagation of an input signal affected by the ASE noise. Assuming the mean power 2σ 2 of the injected noise to be small with respect to the typical power of the signal, we perform a perturbation expansion at first-order in the small parameter σ . The leading order describes the fully deterministic propagation of the signal while the first-order obeys a linear system with stochastic input. Assuming the input ASE to be described by a band-limited white gaussian process, it turns out that, at first order in σ , the input power and frequency are still described by a joint gaussian process. Hence, the linearity of the stochastic system implies that the statistics remain gaussian after the propagation, which simplifies the computation of the output statistics.
The paper is organized as follows. The NLSE is recast in hydrodynamic form in Section 2, where also the small-dispersion regime is considered. In Sec. 3, the statistics of the hydrodynamic variables are computed up to terms of order σ 2 . The perturbative expansion is introduced in general in Sec. 4 and is applied to the particular regimes in the subsequent sections 5 (smalldispersion regime), 6 (CW signal) and 7 (solitonic pulse). In Sec. 8 the perturbative method is used to get equations that govern the evolution of the expected value of the solution. Finally, Sec. 9 is devoted to numerical experiments, in the CW case, that are aimed to test the accuracy of the method.

Hydrodynamic formulation of the NLSE
We assume, as usual, that the electromagnetic field propagates along the fiber according to the NLSE [25] where the subscripts denote derivation with respect to z and t. Here, z > 0 is a longitudinal coordinate of the fiber axis, t ∈ R is the time coordinate corresponding to a time frame moving with the signal group velocity, β 2 ∈ R is the group velocity dispersion parameter and γ ≥ 0 and α ≥ 0 are, respectively, the nonlinear refraction and attenuation parameters. The variable u = u(z,t) is the envelope of the modulation of the electromagnetic field along the fiber. The dispersion parameter β 2 can be either positive (normal dispersion), or negative (anomalous dispersion), giving the NLSE, Eq. (1), a defocusing or focusing character, respectively. We apply to Eq. (1) the transformation where the real quantities ρ ≥ 0, φ and v, which are functions of z and t, are, respectively, power, phase, and angular frequency of the signal u. Substituting Eq. (2) into Eq. (1), and writing separately the equations for the real and imaginary part, we obtain the "hydrodynamic" formulation of the NLSE: Apart from the exchange of space and time coordinates, the presence of the dissipative term −αρ and, of course, different physical constants, system (3) is the hydrodynamic form of Schrödinger equation (for the nonlinear potential V = γρ) that was introduced by E. Madelung in 1926 [15]. They have the form of an Euler system of conservation laws for an isothermal, compressible gas with a "quantum pressure" term β 2 Q, where By considering a reference time T 0 and power P 0 (determined by the duration and mean power of a pulse injected in the fiber), it is known that the dimensionless parameter  (1)). By using dimensionless variables, it is possible to show that the error in neglecting the term β 2 Q in Eq. (3) is of order ε 2 . Then, for ε ≪ 1 we can formally write the approximated system For rigorous results about the limit ε → 0 of system (3) in the defocusing case, one can refer to [24].

Remark 2.1
Condition ε ≪ 1 can be interpreted as a condition on the typical frequency which is generally not satisfied in typical optical fiber communication systems.Thus, while the approximated system (6) can be used for specific applications (e.g., high-power fiber lasers or optical amplifiers), we will employ the complete system (3) to obtain some results of more general validity, suitable also for optical communications. (6) with α = 0 and β 2 > 0 is a hyperbolic system of conservation laws (indeed, it is formally equivalent to the well known equations describing shallow waters). For β 2 < 0 such system is not hyperbolic, which implies that the input-value problem is not well posed in general and leads to unphysical results [26]. This can be clearly seen in the simulations reported in Sec. 9, where the small-dispersion approximation applied to a fiber with anomalous dispersion predicts a non-physical, exponential growth.

Remark 2.2 Equation
System (3), or its small-dispersion approximation, Eq. (6), will be supplemented with input conditions (i.e. conditions at z = 0) of stochastic nature, which will be described in next section.

Stochastic input data
We assume that an optical amplifier injects in the fiber, at z = 0, an input signal of the form where u 0 (t) is the deterministic part of the input signal, and σ Z(t) is the stochastic part, that is the ASE noise. We assume that Z(t) = Z 1 (t) + iZ 2 (t) is a gaussian, band-limited, complex white noise [27], i.e. a stochastic process such that: i) for every fixed t ∈ R, Z 1 (t) and Z 2 (t) have normal distribution with mean 0 and variance 1 (so that the variance of the noise is 2σ 2 ), and are independent; ii) the autocorrelation function is (where · denotes the expectation), corresponding to a power spectral density where ω c > 0 is the cutoff frequency.
Note that 2σ 2 is the mean power of the noise and in the following we will assume

Statistics of the input power and phase
The input condition, Eq. (8), has now to be translated into an input condition for the hydrodynamic system, Eq. (6). This will be done at first order in σ (i.e. up to terms of order σ 2 , that are small compared to the typical pulse power P 0 ). As far as the power is concerned, we can write where, according to Eq. (2), the deterministic part of the input has been decomposed as The stochastic perturbation of the input power, is therefore a gaussian process with 0 mean and (by Eq. (9)) autocorrelation As far as the phase is concerned, we have where the perturbation is a gaussian process with 0 mean and autocorrelation

Statistics of the input frequency
We can now compute the statistics of the input frequency By standard results on differentiated processes [27] we have that V (t) is a gaussian process with 0 mean and autocorrelation function Then, by using Eq. (18), we obtain where We finally compute the cross-correlation between R and V . First of all, from Eqs. (9), (14) and (17), we immediately obtain Then, using

Perturbative treatment of the stochastic problem
Let us rewrite here the stochastic input-value problem given by Eqs. (3), (12) and (19): This problem will be treated perturbatively with respect to the small parameter σ ; then we expand the unknowns as where we assume ρ (1) ∼ v (1) ∼ σ . Substituting into Eq. (23) yields, at leading order, with input conditions (which, of course, means that at leading order we are left with the deterministic problem) and, at first order, where Q (1) is the first-order approximation of Q (see Eq. (4)), given by (linear in ρ (1) ). Note that Eq. (24) is a nonlinear deterministic problem (equivalent to the NLSE) while, once Eq. (24) has been solved, Eq. (25) is a linear stochastic problem.

Remark 4.1
The stochastic problem introduced in this section refers to the case of signal and ASE noise injected at the input of a single span of fiber. The analysis can be easily extended to a multi-span optical system by following the same approach employed in [10]. In this case, at each amplification stage, an independent perturbation of the order of σ (due to the ASE noise generated by the inline amplifier), with statistics as in Eqs. (15), (20) and (22), is added to the (already noisy) Madelung variables ρ and v, which are then propagated through the following span of fiber according to Eq. (23a).

The small-dispersion approximation
Let us consider the small-dispersion approximation (Q ≈ 0) of the input-value problems (24) and (25). The leading order equation, Eq. (24a), becomes that is the "shallow water" system (6). As already mentioned in Remark 2.2, according to the standard theory of first-order systems [26,28], we have that the system is hyperbolic if and only if β 2 > 0. In this case the characteristic velocities are and the characteristic curves are defined by the equations Moreover, we obtain the equations for the Riemann variables 2 ρ (0) ± β 2 γ v (0) . For a = 0, the Riemann variables are constants along the characteristic lines, in which case they are Riemann invariants for the problem. The condition α = 0 can always be assumed on a short fiber stretch, provided that the nonlinear coefficient γ is substituted with an effective coefficient γ eff [10,25]. Assume now that (ρ (0) , v (0) ) is the solution to problem (27), and consider the small-dispersion approximation of the first-order equation (25a): This is a linear, non-autonomous, problem for the unknowns (ρ (1) , v (1) ). For β 2 > 0, its characteristic curves are fixed, and are the same as those of problem (27), i.e. they are still given by Eq. (29). The characteristic equations read as follows: where (1) are the Riemann variables for the first-order problem.

The continuous-wave case
The simplest case in which problems (24) and (25) can be explicitly solved is that of constant deterministic input data ρ 0 and v 0 . We shall also assume α = 0. Under these conditions, problem (24) admits the constant solution Such solution is physically meaningful and corresponds to a continuous-wave (CW) solution of the NLSE (1) (with α = 0) [10] u(z,t) = e −i β 2 v 2 0 2 +γρ 0 z e iv 0 t √ ρ 0 (34) (up to a constant phase). The first-order equation, Eq. (25a), becomes At this point we can notice that we can get rid of the terms containing v 0 by means of the change of variable t → t − β 2 v 0 z (indeed, v 0 just corresponds to a global frequency shift), which leads to In order to make a comparison with the CRLP method [10], we return to the the phase variable φ t = v for which we easily obtain with the stochastic input conditions Note that an arbitrary function of z was set to 0 in Eq. (37a); this can be justified by looking at the Madelung equation for the phase, of whom Eq. (37a) is the linearization. The comparison of Eqs. (37) with the CRLP model is not straightforward since, in the latter, a three-variable perturbation expansion is used to represents the optical signal. However, as shown in [10], the CRLP expansion can be reduced to a simpler amplitude and phase representation by setting the quadrature component equal to zero. This is a reasonable approximation when nonlinearity dominates over dispersion. In this case, the equation governing the propagation of the CRLP amplitude and phase variables is analogous to Eq. (37a), but for obvious differences (a factor of two dividing and multiplying the right-hand side of the first and second equation, respectively) due to the fact that Eq. (37a) refers to a power and phase expansion. The Fourier transform of Eq. (37a) is and the associated transfer matrix T (z, ω), such that can be explicitly computed and reads as follows: This expression holds for β 2 > 0, or β 2 < 0 and ω 2 > 4γρ 0 |β 2 | . The expression of T (z, ω) for β 2 < 0 and ω 2 < 4γρ 0 |β 2 | is easily obtained by using cos(ix) = cosh(x) and sin(ix) = i sinh(x) (x real). In the small-dispersion approximation, Eq. (37a) drastically simplifies into and the corresponding transfer matrix, for β 2 > 0, becomes while the expression for β 2 < 0 can be obtained as explained above.

PSD matrix for the continuous wave
The transfer matrix T (z, ω) allows to obtain the power spectral density (PSD) matrix from its input value G(0, ω), through where T * z (ω) denotes the adjoint matrix (which is just the transpose, in this case). The input PSD matrix can be easily obtained from Eqs. (15), (18) and (21), with u 0 = √ ρ 0 , φ 0 = 0 and v 0 = 0, and reads as follows: Of course, in the "white noise limit" ω c → ∞ the rect function tends to 1.

Numerical results
The accuracy of the proposed approach is tested in the CW case, for which exact analytical solutions for the PSD and pdf of the received signal have been derived. We consider the propagation of a CW signal affected by band-limited AWGN noise through 50 km of lossless fiber. As a first test, we consider a non-zero dispersion-shifted fiber with GVD parameter β 2 = 5.1 ps 2 /km (normal dispersion) and nonlinear coefficient γ = 1.27 W −1 km −1 , a signal power of 10 mW, and a noise bandwidth of 50 GHz. The PSD matrix of the output power and phase variables is evaluated analytically from Eq.   the input one, Eq. (47) (which is diagonal and frequency-independent within the noise bandwidth), testifies the presence of a significant signal-noise interaction. At 20 dB of SNR, there is a perfect agreement between the full model and numerical simulations, while the approximated model, as expected, is accurate only at low frequency. On the other hand, at 10 dB of SNR, a small discrepancy between numerical simulations and the full model can be observed. This is due to the neglected higher-order terms in σ which, at 10 dB of SNR, become relevant. In order to investigate the impact of fiber dispersion on the accuracy of the model, the same analysis is performed also for a standard fiber with a GVD parameter β 2 = −21.7 ps 2 /km (anomalous dispersion), leaving all the other parameters unchanged. Figures 2(a) and 2(b) show the corresponding results for an SNR of 20 and 10 dB, respectively.
In this case, the anomalous dispersion of the standard fiber is responsible for a significant amplification of the PSD of the power variable at some frequencies due to the presence of hyperbolic functions in Eq. (40). This corresponds to an overall increase of the noise power due   to signal-noise interaction. Moreover, the higher value of the modulus of the GVD parameter confines signal-noise interaction to a narrower bandwidth. Also in this case, the accuracy of the full model is confirmed, though the discrepancy at 10 dB of SNR is slightly more relevant than in the previous case because of the aforementioned amplification of the noise power. However, the approximated model is very inaccurate because of the higher dispersion. As a final test, we investigate the ability of the model to reproduce the non-Gaussian distribution of the output signal. We consider the same CW power and the same kind of fibers considered in Figs. 1 and 2, an SNR of 20 dB, and reduce the noise bandwidth to 10 GHz (trying to keep all the noise within a bandwidth in which it is mostly affected by signal-noise interaction, as shown by the PSDs in the previous figures). Figures 3(a) and 3(b) show the contour plots of the joint pdf of the real and imaginary parts of the (normalized and derotated) optical signal u(z,t)e iγρ 0 z / √ ρ 0 at the output of the fiber at an arbitrary time t (as, for a CW input signal, the process is stationary) for the non-zero dispersion-shifted fiber and the standard fiber, respectively. The solid curves are obtained according to the full first-order model, evaluating the PSD matrix of the power and phase variable as discussed before. Given the input conditions and the linearity of the propagation equations, the output variables are jointly Gaussian and their covariance matrix is evaluated by integrating the PSD matrix. The joint distribution of the real and imaginary parts is finally obtained through a change of variables. The dotted curves are ob-tained by using the split-step Fourier method for the propagation of the noisy signal through the fiber, and applying the multicanonical Monte Carlo algorithm (see [29] and references therein) to estimate the pdf with uniform accuracy down to low values. In both cases, there is a good agreement between theory and simulations in the modal region of the pdfs, while some discrepancy appears in the tails. This is, again, due to the neglected higher-order terms, which are more relevant for higher values of the perturbation variables (i.e., in the tails of the pdf). Also in this case, the discrepancy is higher in the standard fiber because of the noise amplification that is possible with anomalous dispersion.

Conclusions
We have introduced a perturbative approach, based on the hydrodynamic formulation of the NLSE, to the nonlinear propagation of signal and noise in an optical fiber. In the smalldispersion approximation, the proposed method corresponds to the semiclassical approximation of the Madelung equations in quantum mechanics. However, we have seen that such approximation is well-behaved only in the defocusing case (normal dispersion), for which the resulting first-order PDE system has a hyperbolic character. On the other hand, by removing this approximation, the method gets a more general validity and its accuracy is greatly improved.
The developed formalism is rather general and, in principle, can be applied to input signals of any form. In this paper we have considered particular input data, namely the CW signal and the solitonic pulse, for which simple explicit equations have been derived. In the CW case, in particular, we have obtained explicit expressions for the PSD and pdf of the received signal that have been compared with numerical simulations (performed by employing the split-step Fourier method). The numerical results confirm the soundness of the proposed approach.