## 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.

© 2013 Optical Society of America

## 1. 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–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.

## 2. Four Wave Mixing (FWM)

#### 2.1 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*and

_{q}*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*. 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

_{r}*P*of the FWM product is given by:

_{pqr}*P*(

_{i}*i*=

*p*,

*q*,

*r*) represents the input peak power at the frequencies

*f*=

_{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 (

_{i}*P*=

_{i}*P*for

_{in}*i*= 1,2,…,

*N*). In (1),

_{ch}*γ*is the nonlinear coefficient of the fiber,

*α*is the fiber loss coefficient,

*L*is the total fiber length, ${L}_{eff}=\left(1-{e}^{-\alpha L}\right)/\alpha $is the effective length of the fiber,

*d*is the degeneracy factor (

_{pqr}*d*= 3 when

_{pqr}*p*=

*q*,

*d*= 6 when

_{pqr}*p*≠

*q*) and

*η*is the mixing efficiency given by:

*β*represents the phase mismatch and may be expressed in terms of the channel frequencies

*f*. For typical values of fiber’s chromatic dispersion coefficient

_{i}*D*, d

*D*/dλ (λ is the wavelength of the signal) and channel spacing Δ

*f*, phase mismatch can be approximated as follows:

*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:

*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*) is the input peak power of the central channel

_{in}*z*. In (3), (4) variables

*δ*,

*I*

_{m}_{,}

*are given by [6]:*

_{s}*θ*=

_{pqr}*θ*+

_{p}*θ*-

_{q}*θ*is the phase of the FWM noise generated from a channel combination (

_{r}*p*,

*q*,

*r*) with

*p*,

*q*,

*r*ϵ [1 . . .

*N*] and

_{ch}*r*=

*p*+

*q*-

*z*. Random phases

*θ*due to channel

_{i}*i*phase noise can be assumed mutually independent and uniformly distributed in [0, 2

*π*]. Furthermore,

*B*= 0 or

_{i}*B*= 1 is the bit value of channel

_{i}*i*.

#### 2.2 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*in the mark and the space state respectively by transforming the FWM efficiency of (2a) as follows [12]:

_{s}*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*and

_{m}*I*from systems parameters (see Appendix A). Their applicability assumes the FWM efficiency given by Eq. (2a).

_{s}## 3. 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*or minima

_{i}*M*’

*= min*

_{n}*X*(extremes),

_{i}*i*= 1, 2,…,

*n*as n→∞ is investigated by the EVT. The cdfs

*F*(

_{n}*x*) of

*M*and

_{n}*F*'(

_{n}*x*) of

*M*’

*can be written as a function of*

_{n}*F*(

*x*):

*X*is unknown and

_{i}*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

*μ*and scale parameters

_{n}*σ*, depending on

_{n}*n*, such that:

*ξ*is the shape parameter. For

*ξ*= 0, the cdf of the maxima

*F*is limited to the Gumbel distribution or simply the extreme value distribution:

_{n}*α*and

_{n}*u*depend on

_{n}*n*and are related to

*μ*and

_{n}*σ*. Moreover, in the cases

_{n}*ξ*> 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}## 4. 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*of (13) or the parameters

_{n}*a*’ and

_{n}*u*’ of (14), the sorted set of maxima (minima)

_{n}*x*

_{1},

*x*

_{2},…,

*x*such that

_{N}*x*

_{1}≤

*x*

_{2}… ≤

*x*, are plotted against the linearized probability

_{N}*y*= -ln[-ln[Λ(

_{i}*x*)]] or

_{i}*y*’ = ln[-ln[1-Λ(

_{i}*x*)]] respectively, where Λ(

_{i}*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 $pdf=B{e}^{\pm b{I}_{m,s}}=\mp bA{e}^{\pm b{I}_{m,s}}$leading to an exponential cdf $cdf=1-A{e}^{\pm b{I}_{m,s}}$. Furthermore, one can easily show that exponential cdfs always satisfy von-Mises convergence criterion, $\underset{x\to \infty}{\mathrm{lim}}\frac{d}{dx}\left(\frac{1-cdf}{pdf}\right)=0$, which determines the distributions *F*(*x*) 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*’. Hence, parameters

_{n}*a*and

_{n}*u*or

_{n}*a*’ and

_{n}*u*’ are obtained through linear fit such that

_{n}*y*=

_{i}*a*(

_{n}*x*-

_{i}*u*) or

_{n}*y*’ =

_{i}*a*’ (

_{n}*x*-

_{i}*u*’). To compute the BER, one can use:

_{n}*S*and

_{m}*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

_{s}*Q*is the decision level. The terms in the last equality of (15) can be easily estimated using the photocurrents’ cdfs which in turn is related to

*F*’ and

_{n}*F*through Eqs. (9a) and (9b).

_{n}Assuming that the maxima of the photocurrents at the space state and the minima of the photocurrent at the mark state follow the distributions ${F}_{n,\text{\hspace{0.17em}}G}$ and ${F}_{n,\text{\hspace{0.17em}}G}^{\text{'}}$ (Eqs. (13) and (14)), as mentioned before, and using (9a) and (9b), the cdfs of the photocurrents are given by:

*Q*is chosen such as to minimize the BER.

## 5. 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.

#### 5.1 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*= 7.83 × 10

_{n}^{−5},

*a*’ = 2.99 × 10

_{n}^{4}and

*u*’ = 6.69 × 10

_{n}^{−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].

#### 5.2 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*= 32,

_{ch}*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}.

#### 5.3 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).

Figure 5 depicts 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*= 4dBm, the BER equals to 3 × 10

_{in}^{−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.

## 6. 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.

## Appendix A

The optical fields amplitudes, *E*^{(}^{m}^{)} and *E*^{(}^{s}^{)}, at a given channel *z* in the mark and the space states, respectively, are given by:

*P*and

_{z}*θ*are the input peak power and the phase in the mark state, respectively, of given channel

_{z}*z*, while

*P*and

_{pqr}*θ*=

_{pqr}*θ*+

_{p}*θ*-

_{q}*θ*are the peak power and the phase, respectively, of the FWM noise generated from a channel combination (

_{r}*p*,

*q*,

*r*).

The photocurrent at the detector is written as:

*P*of the FWM product is:

_{pqr}*η*, is given by:

*D*, d

*D*/d

*λ*(

*λ*is the wavelength of the signal) and channel spacing Δ

*f*, phase mismatch can be approximated as follows:

*c*is the speed of light in vacuum.

For *D* ≥ 2ps/nm/Km and channel spacing Δ*f* ≥ 10GHz, one can assume that Δβ>>α while for large *L*, exp(-αL)<<1. In this case, the photocurrent in the mark state can be simplified to:

*P*=

_{p}*P*=

_{q}*P*=

_{r}*P*, the following equation can be easily derived:

_{in}*α*,

*L*,

*γ*,

*λ*,

*D*,

*P*and Δ

_{in}*f*. By doing this, system parameters can be removed outside summation and packed into the variable

*δ*.

In the case of a multispan system, the efficiency of the four-wave mixing is given by:

^{2}(

*M*Δ

*βL*/2)/ sin

^{2}(Δ

*βL*/2) is forced to participate in the summation (due to the presence of

*p*,

*q*and

*r*in Δ

*β*), trapping system parameters. This prohibits the isolation of system parameters from the rest terms of the photocurrent.

In the space state, it is straightforward to show that:

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.

## References and links

**1. **B. Mukherjee, *Optical WDM Networks* (Springer, 2006).

**2. **C.-H. Lee, S.-M. Lee, K.-M. Choi, J.-H. Moon, S.-G. Mun, K.-T. Jeong, J. H. Kim, and B. Kim, “WDM-PON experiences in Korea [Invited],” J. Opt. Netw. **6**(5), 451–464 (2007). [CrossRef]

**3. **C. H. Lee, W. V. Sorin, and B. Y. Kim, “Fiber to the Home Using a PON Infrastructure,” Lightwave Technology, Journalism **24**, 4568–4583 (2006).

**4. **I. Neokosmidis, T. Kamalakis, A. Chipouras, and T. Sphicopoulos, “New techniques for the suppression of the four-wave mixing-induced distortion in nonzero dispersion fiber WDM systems,” J. Lightwave Technol. **23**, 1137–1144 (2005).

**5. **M. Eiselt, “Limits on WDM Systems Due to Four-Wave Mixing: A Statistical Approach,” J. Lightwave Technol. **17**(11), 2261–2267 (1999). [CrossRef]

**6. **I. Neokosmidis, T. Kamalakis, A. Chipouras, and T. Sphicopoulos, “Evaluation by Monte Carlo Simulations of the Power Limits and Bit-Error Rate Degradation in Wavelength-Division Multiplexing Networks Caused by Four-Wave Mixing,” Appl. Opt. **43**(26), 5023–5032 (2004). [CrossRef] [PubMed]

**7. **I. Neokosmidis, T. Kamalakis, A. Chipouras, and T. Sphicopoulos, “Estimation of the four-wave mixing noise probability-density function by the multicanonical Monte Carlo method,” Opt. Lett. **30**(1), 11–13 (2005). [CrossRef] [PubMed]

**8. **R. Holzlöhner and C. R. Menyuk, “Use of multicanonical Monte Carlo simulations to obtain accurate bit error rates in optical communications systems,” Opt. Lett. **28**(20), 1894–1896 (2003). [CrossRef] [PubMed]

**9. **E. Castillo, A. S. Hadi, N. Balakrishman, and J. M. Sarabia, *Extreme value and related models with applications in engineering and science* (John Wiley and Sons, 2005).

**10. **S. Savory, F. Payne, and A. Hadjifotiou, “Estimating Outages Due to Polarization Mode Dispersion Using Extreme Value Statistics,” J. Lightwave Technol. **24**(11), 3907–3913 (2006). [CrossRef]

**11. **K. Inoue, K. Nakanishi, K. Oda, and H. Toba, “Crosstalk and power penalty due to fiber four-wave mixing in multichannel transmissions,” Lightwave Technology, Journalism **12**, 1423–1439 (1994).

**12. **K. Inoue, “Phase-mismatching characteristic of four-wave mixing in fiber lines with multistage optical amplifiers,” Opt. Lett. **17**(11), 801–803 (1992). [CrossRef] [PubMed]

**13. **G. H. Einarsson, *Principles of Lightwave Communications* (John Wiley & Sons, 1996).