Semianalytical estimation of the four-wave mixing noise based on extreme value theory

Abstract: Four-wave mixing (FWM) is one the limiting factors for existing and future wavelength division multiplexed optical networks. A semianalytical method based on Monte Carlo and Extreme Value theory is proposed and applied to study the influence of the FWM noise on the performance of WDM systems. The statistical behavior of the FWM noise is investigated while the Bit-Error rate is calculated for various combinations of the design parameters and for both single and multiple span WDM systems. The semianalytical method is also compared to the Multicanonical Monte Carlo (MCMC) method showing the same efficiency and accuracy with the former providing however the advantage of deriving closed-form approximations for the cumulative distribution functions of the photocurrents in the mark and the space state and the BER.


Introduction
It is widely accepted that multimedia applications and broadband services will play a significant role in the operation of future and already installed Wavelength Division Multiplexing (WDM) Optical Networks [1] due to the extremely huge volume of data that should be transmitted.These phenomena will be further enhanced as optical fibers get closer to the customers' premises, an architecture known as Fiber To The Home (FTTH) [2,3].The increased bit-rate of access networks will enable the extended use of capacity starving applications leading to the aggregation of high data volumes in transport and backbone networks.As a consequence, optical backbone networks must have the demanded characteristics (e.g.capacity) to support this vast data traffic in order to provide the quality of service required for these upgraded applications.Although, WDM optical networks concentrate several advantages such as low losses and high capacity, there are numerous linear and nonlinear effects degrading their performance and thus should be accurately evaluated.
Nonlinear effects such as self-and cross-phase modulation as well as four-wave mixing (FWM) seem to be more severe than linear ones since they cannot be easily compensated for by optical amplifiers and chromatic dispersion compensators.Both cross-phase modulation and FWM cause interference between channels of different wavelengths, resulting in an upper power limit for each WDM channel.However, as the channel spacing and/or the chromatic dispersion of a WDM network decrease, FWM greatly influences its performance.Using split-step Fourier simulation, it was shown that FWM-induced distortion dominates over both SPM and XPM effects in a WDM link with low dispersion fibers [4].Hence, a study of such systems is useful since low dispersion fibers (G.655 fibers) can provide many advantages (e.g., reduced need for dispersion compensation modules, etc.) and are used in numerous already installed links especially in transport and backbone networks.
Several techniques have been used for the evaluation of the Bit Error Rate under the influence of FWM [4][5][6][7].However, these techniques suffer from several drawbacks.The Split Step Fourier method [4] cannot be used for the evaluation of the BER and is limited to the estimation of the Q-factor under the assumption of Gaussian distributed noises.In [5], Monte Carlo simulations were performed, however the accuracy of the obtained results is low especially in the range of low BER values due to time constraints and the use of a symmetrical exponential approximation.The accuracy was slightly improved at [6] at the cost of increased time consumption.The Multicanonical Monte Carlo (MCMC) method [7,8] can greatly reduce the time required for the simulations but it cannot give a closed form formula for the probability density functions (pdfs) of the photocurrents that are used in the estimation of the BER.
In this paper, a hybrid method is proposed for the derivation of semi-analytical expressions for both the BER and the cumulative distribution functions (cdfs) of the photocurrents using Extreme Values Theory (EVT) [9,10] enabling the quick and accurate estimation of the FWM noise impact.

Single span systems
A WDM system with equally spaced N ch channels and amplitude-shift keying modulation, which is the most frequently used modulation scheme, is considered.All signals are assumed to be copolarized and synchronized, which represents a worst-case situation.Only the performance of the central channel has to be investigated since, in this case, the energy conservation requirement is satisfied by the largest number of frequency combinations.
FWM introduce intensity fluctuations that depend on the intensity of neighboring channels, resulting into interchannel interference.In detail, when three waves of frequencies f p , f q and f r (p ≠ q, r) interact due to the third-order electric susceptibility, a product wave is generated at frequency f pqr = f p + f q -f r .In case of systems with in-phase and equally spaced channels, most of the generated components will be located at the initial frequencies leading to induced distortion.The output power P pqr of the FWM product is given by: where P i (i = p, q, r) represents the input peak power at the frequencies f i = ω i /2π in the mark state.In a WDM system it can be assumed that all the peak powers at the mark state are equal (P i = P in for i = 1,2,…,N ch ).In (1), γ is the nonlinear coefficient of the fiber, α is the fiber loss coefficient, L is the total fiber length, ( ) is the effective length of the fiber, d pqr is the degeneracy factor (d pqr = 3 when p = q, d pqr = 6 when p≠q) and η is the mixing efficiency given by: ( ) ( ) In (2b), Δβ represents the phase mismatch and may be expressed in terms of the channel frequencies f i .For typical values of fiber's chromatic dispersion coefficient D, dD/dλ (λ is the wavelength of the signal) and channel spacing Δf, phase mismatch can be approximated as follows: ( )( ) where c is the speed of light in vacuum.
In practical applications, it can be assumed that Δβ>>α which generally holds for D ≥ 2ps/nm/Km and channel spacing Δf ≥ 10GHz.For large L one can also use the fact that exp(-αL)<<1.In this case, the expressions of the FWM noise photocurrents are: for the mark (the signal bit of the central channel is 1) and the space (the signal bit of the central channel is 0) state respectively.Note that k is the receiver's responsivity, E (m) and E (s)  are the amplitudes of the optical fields in the mark and the space state, respectively [11] and P z ( = P in ) is the input peak power of the central channel z.In (3), (4) variables δ, I m,s are given by [6]: ( ) for the mark and the space state respectively.Note that θ pqr = θ p + θ q -θ r is the phase of the FWM noise generated from a channel combination (p, q, r) with p, q, r ϵ [1 . ..N ch ] and r = p + q -z.Random phases θ i due to channel i phase noise can be assumed mutually independent and uniformly distributed in [0, 2π].Furthermore, B i = 0 or B i = 1 is the bit value of channel i.

Multispan systems
The analysis carried out in the previous subsection assumes a single-span system.However, practical systems consist of multiple spans.Moreover, if optical amplifiers are used in multispan WDM systems to compensate for the optical losses, then additional FWM noise products are generated in each span.Thus expansion of the above analysis is required and is of great importance.An easy way to take into account multiple spans is to approximate the pdf of the total FWM noise by the convolution of the individual (in each span) FWM noise pdfs.This approximation is based on the assumption that the phases of the channels at the beginning of each span, the phases of the products in different spans as well as FWM noise products in each span are independent [5].Also, the net dispersion in each span causes a walk-off of neighboring channels by at least one bit period.These considerations presuppose that the dispersion of each fiber span differs slightly and that the lengths of the fibers used in the span are different.
It is thus obvious that these assumptions may not be valid in all cases.If for example, the dispersion is completely compensated for in each span; then the phases of the channels cannot be considered independent.Hence, a more general approach should be followed.
In such cases and assuming equal spans, one can compute the photocurrents S m and S s in the mark and the space state respectively by transforming the FWM efficiency of (2a) as follows [12]: where M is the total span number and Δβ is the phase mismatch factor of each product.This essentially results in multiplying each term of the FWM products sum by sin(MΔβL/2)/sin(ΔβL/2).Doing so, however, destroys the independence of the auxiliary variables I m and I s from systems parameters (see Appendix A).Their applicability assumes the FWM efficiency given by Eq. (2a).

Basics of extreme value theory
In this section, the theoretical background of the EVT will be presented.To begin with, suppose a set of n independent and identically distributed random variables Χ 1 , Χ 2 ,…,X n , with a pdf f(x) and a cdf F(x).The statistical behavior of the maxima M n = max X i or minima M' n = min X i (extremes), i = 1, 2,…,n as n→∞ is investigated by the EVT.The cdfs F n (x) of M n and F n '(x) of M' n can be written as a function of F(x): Equations ( 9a) and (9b) are not really useful in practical cases since the initial cdf of X i is unknown and F n (x) depends on n.However, in the case of maxima, as n→∞ and according to Fisher-Tippett theorem [10], there are location parameters μ n and scale parameters σ n , depending on n, such that: where: ( ) ( ) It is straightforward to show that the cdf of the maxima is: ( ) where ξ is the shape parameter.For ξ = 0, the cdf of the maxima F n is limited to the Gumbel distribution or simply the extreme value distribution: It should be noted that α n and u n depend on n and are related to μ n and σ n .Moreover, in the cases ξ > 0 and ξ < 0, F n (x) follow Fréchet or Weibull distributions respectively.In the same manner, one can easily derive the following expression for the cdf F n '(x) of the minima M' n :

Using EVT in the BER calculation
In this section, the proposed framework for the estimation of system performance in terms of the BER will be described.In fact, EVT will be used to quickly find semi-analytical expressions for both the photocurrents' cdfs in the mark and the space state and the BER.
According to the proposed method, a set of p = n × N values of the photocurrent S (m) (or S (s) ) are produced using Monte Carlo simulations on Eqs. ( 3)- (7).The p photocurrent values are segmented into N groups of n samples and each group forms the sequence of the independent and identically distributed (iid) random variables X 1 , X 2 ,…, X n .The maximum (minimum) value of each group is then calculated for the space (mark) state.The obtained maxima (minima) are sorted in the ascending order as analyzed in [10].In order to estimate the parameters a n and u n of ( 13) or the parameters a n ' and u n ' of (14), the sorted set of maxima (minima) ] respectively, where Λ(x i ) = i / (1 + N) is the empirical distribution of the maxima (minima).
In [6], the pdfs of the variable Is and Im (left part) and thus the pdfs of the photocurrents were shown to be of the form leading to a Gumbel distribution for the extremes.This can also be confirmed by the example application provided in section II.D of [10] describing how an exponential distribution function leads to a Gumbel extreme distribution.Using these facts, in the following the Gumbel distribution is used as the cdf of the maxima F n or minima F n '.Hence, parameters a n and u n or a n ' and u n ' are obtained through linear fit such that y i = a n (x i -u n ) or y i ' = a n ' (x iu n ').To compute the BER, one can use: where ( ) ( ) , m s S S f f are the pdfs of the photocurrents, S m and S s , in the mark and the space state respectively.The integrals are the probabilities that an error will occur in the space and the mark states, respectively and Q is the decision level.The terms in the last equality of (15) can be easily estimated using the photocurrents' cdfs which in turn related to F n ' and F n through Eqs.(9a) and (9b).
Assuming that the maxima of the photocurrents at the space state and the minima of the photocurrent at the mark state follow the distributions , n G F and ' , n G F (Eqs. ( 13) and ( 14)), as mentioned before, and using (9a) and (9b), the cdfs of the photocurrents are given by: Hence, Eq. ( 15) can be transformed as follows: where Q is chosen such as to minimize the BER.

Results and discussion
We applied the proposed hybrid Extreme value / MC method for the estimation of the impact of the FWM noise on the performance of WDM systems.Several values for both the iterations and the samples were examined, however there was no good convergence between the Gumbel distribution and the photocurrents distribution.In the following, N = 1000 iterations were used, each involving the generation of 100 samples [10] leading to a good convergence as shown in the next subsection.In the remainder of the paper the parameters used in the calculations are as follows: λ = 1.55μm, c = 3 × 10 8 m/s, L = 80Km, α = 0.2dB/Km, γ = 2.4 (Km × W) −1 and k = 1.28 A/W.

Comparison with multicanonical Monte Carlo method
In order to assess WDM systems performance using the proposed hybrid method, one should initially investigate its accuracy.Towards this direction, the Extreme value / MC method is compared to the multicanonical MC (MCMC) method.
The results obtained by the hybrid method for the central channel (z = 8) of a 16-channel WDM system are plotted as points in Fig. 1 for various values of the chromatic dispersion and the channel spacing.Shown as lines are the BER obtained by the MCMC method.From Fig. 1, it can be deduced that the results of the proposed hybrid method are almost identical to those of the MCMC method proving its high accuracy.It is interesting to note that although both methods can be used to accurately estimate even very low BER values (i.e., 10 −15 ) with almost same computation time requirements, the proposed method has a big advantage; It provides a semi-analytical expressions for the cdfs in the mark (Eqs.( 14) and ( 16)) and the space state (Eqs.( 13) and ( 17)) respectively as well as for the BER.For example, the normalized cdfs of a 16 channels single span WDM system with D = 5ps/km/nm, Δf = 25GHz and P in = 9dBm obtained by the hybrid method are illustrated in Fig. 2. In this case the parameters of Eqs. ( 13) and ( 14) are found a n = 2.12 × 10 4 , u n = 7.83 × 10 −5 , a n ' = 2.99 × 10 4 and u n ' = 6.69 × 10 −5 .This is important since it enables further processing of the obtained results without the need of repeating simulations.For example, one can approximately but easily estimate the statistics of the photocurrents in a WDM multispan system (for various values of spans) using the semi-analytical expressions obtained once for the single span system along with the convolutional approach described in [5].
Furthermore, the advantage of semi-analytical expressions can be used in order to incorporate other noises, such as thermal noise, influencing actual systems performance.To include these noises, one can use a technique based on the moment generating function (MGF) of the decision variable whereas the error probability can be estimated by use of the saddle-point approximation through the MGF as described in [13].

Single span systems
In this subsection, the proposed method is used in the case of a single span WDM system.Although single span systems can be considered as a simplification of practical systems, they can be proved useful and powerful for the investigation of FWM noise statistical characteristics providing a first glance in FWM impact on systems performance.
Figures 3(a)-3(c) depict the BER as a function of the input peak power P in for various values of the number of channels, the chromatic dispersion and the channel spacing using (18).As expected, when the dispersion or the channel spacing is increased lower values of the BER are obtained.Similar behavior is observed as the number of channels is reduced.
For example, for N ch = 16 channels, channel spacing Δf = 25 GHz and input peak power P in = 4 dBm, an increase in chromatic dispersion coefficient D from 2 to 5ps/nm/km causes a reduction of the BER from 3 × 10 −3 to 10 −7 .
Such plots are important for the design of practical systems.They are useful in determining the maximum input power allowed in a WDM link, given its characteristics, i.e., channel spacing Δf, fiber dispersion coefficient D, and the total number of channels Ν ch .For example, for N ch = 32, D = 2ps/nm/km, and Δf = 50 GHz, the input peak power of each channel must not exceed 5dBm if the BER is not to exceed the threshold of 10 −9 .

Multispan systems
In order to investigate the more practical case of multispan systems, a series of hybrid Extreme Value / MC simulations is performed using this time the FWM efficiency of (8).
In Fig. 4, the results for the cdf of the FWM photocurrent in a multispan system are presented for M = 2, 4, 8 spans, N ch = 32 channels, chromatic dispersion coefficient D = 2ps/km/nm, channel spacing Δf = 50GHz and input peak power 6 dBm.Although, Fig. 4 cannot be used in order to extract useful information, it can provide a first insight on the impact of multiple spans on the performance of WDM systems.
As shown in Fig. 4, an increase in the number of spans causes a consequent decrease of both cdfs (mark and space) slopes.This can be an indication of the detrimental influence of the number of spans on systems performance.
To confirm the negative impact of number of spans as well as to derive detailed results regarding FWM induced distortion in a multispan WDM system, the BER is estimated as a function of the input peak power (Fig. 5).5, an increase in the number of spans results in an increase of the BER.For example, in case of P in = 4dBm, the BER equals to 3 × 10 −11 , 10 −8 and 8 × 10 −6 for M = 2, 4 and 8 spans respectively.In other words, an increase of the number of spans from 2 to 4, 4 to 8 and 2 to 8 spans, imposes a power penalty of almost 1, 1.5 and 2.5dB respectively in the case of BER = 10 −9 .This can be attributed to the fact that then additional FWM noise products are generated in each span while the already existing ones are further enhanced.

Conclusion
In this paper, a method based on MC and extreme value theory has been proposed for the semianalytical estimation of the impact of FWM noise in WDM optical systems.The presented method exhibits the same accuracy and computational time requirements with the MCMC while providing close-form BER formulas.This method can also be used to study the statistical characteristics of the FWM noise and provide semi-analytical expressions for the BER and the cdfs of the photocurrents in both single and multiple span systems.
Since P pqr is included in the summations of Eq. (A12), the same conclusion for the multispan system can be derived in the same manner.
one can easily show that exponential cdfs always satisfy vonthe distributions F(x)

Fig. 1 .
Fig.1.Comparison between the proposed hybrid method and the MCMC method for Nch = 16 channels.

Fig. 3 .
Fig. 3. Bit Error Rate as a function of the input peak power P in for a) N ch = 8, b) N ch = 16 and c) N ch = 32 channels.D 1 = 2ps/km/nm and D 2 = 5ps/km/nm.

Fig. 5 .
Fig. 5. Bit Error Rate as a function of the input peak power Pin of multispan WDM systems with N ch = 32 channels, chromatic dispersion coefficient D = 2ps/km/nm, channel spacing Δf = 50GHz and various values of the number of spans.

Figure 5
Figure5depicts the BER as a function of the transmitter's power for various values of the number of spans.The multispan WDM systems under investigation are assumed to have N ch = 32 channels, dispersion coefficient D = 2ps/km/nm and channel spacing Δf = 50GHz.As shown in Fig.5, an increase in the number of spans results in an increase of the BER.For example, in case of P in = 4dBm, the BER equals to 3 × 10 −11 , 10 −8 and 8 × 10 −6 for M = 2, 4 and 8 spans respectively.In other words, an increase of the number of spans from 2 to 4, 4 to