As the bit rate and the transmission distance of optical fiber communication systems continue to increase, the signal distortion caused by polarization mode dispersion (PMD) in optical fiber channels is becoming a major limitation of system performance improvement. Optical compensators are effective to reduce the system degrading effects caused by PMD.¹ Recently, we proposed a new higher-order optical PMD compensator in Ref. 2 where the outage probability (the probability of the PMD-induced bit error rate (BER) exceeding the allowed value ${10}^{-12}$ ) is used to evaluate the efficiency of our proposed compensator. Using standard Monte Carlo (MC) simulation to produce the optical fiber channel samples,² the outage probability of the compensated systems is estimated down to less than ${10}^{-3}$ . In practical system design, one demands that the PMD-induced outage probability is typically ${10}^{-6}$ or less. Estimating the outage probability of the compensated systems down to extremely low values can give the system designer more meaningful information. Using MC as the simulation tool cannot meet this requirement because an extremely large number of system configuration samples must be produced to obtain a reliable estimate.

Markov chain Monte Carlo (MCMC) simulation can provide a fast and accurate estimate of outage probability down to extremely low values. MCMC is based on a series of Markov chain fiber channel samples and the simulated annealing algorithm introduced by Metropolis ^{3, 4, 5} Secondini have applied this approach to evaluate the outage probability in the first-order compensated system.⁶ In this letter, the application of MCMC to accurately estimating the outage probability of the optical communication system compensated by higher-order optical PMD compensators is described. Markov chain fiber channel samples can drastically reduce the number of simulations required to calculate the outage probability of the compensated systems down to extremely low values. Using MCMC, we estimate the outage probability of the optical system compensated by our proposed higher-order optical PMD compensator down to ${10}^{-7}$ . Our proposed compensator always outperforms the referenced higher-order optical compensation scheme. The accuracy of MCMC is affirmed, too.

The principle states of polarization (PSPs) and the differential group delay (DGD) in the fiber channel that determine the PMD effects in the optical communication system are described by the PMD vector.¹ The frequency-dependent PMD vector can be defined in terms of a Taylor expansion in angular frequency $\omega$ . Second-order PMD is decided by the change of DGD with wavelength (PCD) and the rotation of PSPs with frequency (PSP depolarization).¹ PCD does not contribute significantly to the system degrading caused by second-order PMD.^{1, 7} The frequency-dependent eigenvector of the transmission matrix of the PMD-dominated fiber channel tends to process at a nearly constant PSP depolarization rate within a certain limited bandwidth. Compensating the PSP depolarization within certain limited bandwidth around the carrier frequency can mitigate the higher-order PMD degrading effects effectively.^{2, 8} Based on this theory, we proposed the new higher-order optical PMD compensator, which requires a simpler control algorithm than that for the typical higher-order optical PMD compensators and can always have better performance than first-order optical compensators.² The transmission matrix of the proposed higher-order optical compensator is described in detail in Ref. 2.

The PMD-induced outage probability $P_{out}$ of the system can be described as:

$X={\log }_{10}\left(E\right)$ , where $E$ is the BER of the optical system. $I\left(X\right)$ is the indicator function. $I\left(X\right)=u(X+12)$ , where $u(\bullet )$ is the unit step function. $p\left(X\right)$ is the probability density function (pdf) of $X$ . In order to estimate the PMD-induced $P_{out}$ of the system, a series of system configuration samples describing the PMD effects in the fiber channel must be produced. We use a concatenation of $M$ birefringent fibers with random rotation relative to each other as the optical fiber channel model.⁹ In Ref. 2, the fiber channel samples are produced by MC. The length of each birefringent segment $l_j$ is randomly generated from a Gaussian distribution and the rotation angle ${\alpha }_j$ is a uniformly distributed random variable in $[0,2\pi ]$ (Ref. 9). We also consider the temperature fluctuations along the fiber. The phase angle change ${\varphi }_j$ is added to the transmission signal in each birefringent segment,⁹ and it is uniformly distributed in $[0,2\pi ]$ . So the fiber channel model is decided by $(\mathbf{L},\mathbf{\alpha },\mathbf{\phi })=(l_1,l_2,\dots ,l_M,{\alpha }_1,{\alpha }_2,\dots ,{\alpha }_M,{\varphi }_1,{\varphi }_2,\dots ,{\varphi }_M)$ , with the random vectors $\mathbf{L}=(l_1,l_2,\dots ,l_M),\mathbf{\alpha }=({\alpha }_1,{\mathbf{\alpha }}_2,\dots ,{\alpha }_M)$ , and $\mathbf{\phi }=({\varphi }_1,{\varphi }_2,\dots ,{\varphi }_M)$ . $X=f(\mathbf{L},\mathbf{\alpha },\mathbf{\phi })$ and the produced samples $X_1,X_2,\dots ,X_N$ used to estimate $P_{out}$ are independent of each other.

Now we use MCMC to produce the PMD-dominated optical fiber channel samples. The main idea of the MCMC method is to rely on dependent Markov chain sequences, which have a limiting distribution corresponding to a distribution of interest.⁴ A simulated annealing algorithm is used to produce the required Markov chain samples $X_i$ ’s. In order to cause large DGD events to occur more frequently, the Markov chain $\left\{X_i\right\}$ is produced according to the biased pdf $p^{*}\left(X\right)$ instead of $p\left(X\right)$ with $p^{*}\left(X\right)\propto e^{X∕T}$ (Refs. 3, 4) $i=1,2,\dots ,N$ , where $N$ is the number of samples produced. Increasing the number of outage events can reduce the number of simulations required. We can choose $p^{*}\left(X\right)=Cp\left(X\right)e^{X∕T}$ (Ref. 6). $C$ is a constant, and $T$ is the temperature determining the probability enhancement of outage events. The samples $X_i$ ’s are decided by the three random vector samples ${\mathbf{L}}_i$ ’s, ${\mathbf{\alpha }}_i$ ’s, and ${\mathbf{\phi }}_i$ ’s. That is, $X_i=f({\mathbf{L}}_i,{\mathbf{\alpha }}_i,{\mathbf{\phi }}_i)$ with ${\mathbf{L}}_i=(l_{i1},l_{i2},\dots ,l_{iM}),{\mathbf{\alpha }}_i=({\alpha }_{i1},{\alpha }_{i2},\dots ,{\alpha }_{iM})$ , and ${\mathbf{\phi }}_i=({\phi }_{i1},{\phi }_{i2},\dots ,{\phi }_{iM})$ . Starting from ${\mathbf{L}}_1,{\mathbf{\alpha }}_1$ , and ${\mathbf{\phi }}_1(X_1$ is determined), generate ${\stackrel{\mathbf{̃}}{\mathbf{L}}}_{i+1},{\tilde{\mathbf{\alpha }}}_{i+1}$ , and ${\tilde{\mathbf{\phi }}}_{i+1}$ from ${\stackrel{\mathbf{̃}}{\mathbf{L}}}_{i+1}={\mathbf{L}}_i+\mathrm{\Delta }{\mathbf{L}}_i,{\tilde{\mathbf{\alpha }}}_{i+1}={\mathbf{\alpha }}_i+\mathrm{\Delta }{\mathbf{\alpha }}_i$ , and ${\tilde{\mathbf{\phi }}}_{i+1}={\mathbf{\phi }}_i+\mathrm{\Delta }{\mathbf{\phi }}_i$ , where $\mathrm{\Delta }{\mathbf{L}}_i=(\mathrm{\Delta }{l}_{i1},\mathrm{\Delta }{l}_{i2},\dots ,\mathrm{\Delta }{l}_{iM}),\mathrm{\Delta }{\mathbf{\alpha }}_i=(\mathrm{\Delta }{\alpha }_{i1},\mathrm{\Delta }{\alpha }_{i2},\dots ,\mathrm{\Delta }{\alpha }_{iM})$ , and $\mathrm{\Delta }{\mathbf{\phi }}_i=(\mathrm{\Delta }{\varphi }_{i1},\mathrm{\Delta }{\varphi }_{i2},\dots ,\mathrm{\Delta }{\varphi }_{iM})$ are the random perturbations with Gaussian distribution independent of ${\mathbf{L}}_i,{\mathbf{\alpha }}_i$ , and ${\mathbf{\phi }}_i$ .

The corresponding sample ${\tilde{X}}_{i+1}=f({\stackrel{\mathbf{̃}}{\mathbf{L}}}_{i+1},{\tilde{\mathbf{\alpha }}}_{i+1},{\tilde{\mathbf{\phi }}}_{i+1})$ can be calculated from the produced fiber channel sample. By simulated annealing algorithm,^{3, 4} the fiber channel samples ( ${\mathbf{L}}_i$ ’s, ${\mathbf{\alpha }}_i$ ’s, and ${\mathbf{\phi }}_i$ ’s) and the $X_i$ ’s (BER samples) are accepted as follows:

with $\rho =\min \left[\exp \right({\tilde{X}}_{i+1}-X_i∕T),1]$ . In our simulations, $\mathrm{\Delta }{l}_{ij},\mathrm{\Delta }{\alpha }_{ij}$ , and $\mathrm{\Delta }{\varphi }_{ij}$ are all random variables with zero mean values. The variance of $\mathrm{\Delta }{l}_{ij}$ is $(0.3\overline{l}*{\sigma }_X{)}^2∕M$ . The average length $\overline{l}$ of the birefringent fiber element decides the mean DGD of the system when the PMD coefficient of the fiber channel is determined. The variances of $\mathrm{\Delta }{\alpha }_{ij}$ and $\mathrm{\Delta }{\varphi }_{ij}$ are both $(\pi *{\sigma }_X{)}^2∕M$ . Here ${\sigma }_X$ means the average step moving in the $X$ space.

The Markov chain $\left\{X_i\right\}$ produced using the method above has a limiting pdf of $p^{*}\left(X\right)$ as $N$ is increasing.^{3, 4} From Eq. 1, $P_{out}$ can be estimated⁶ using the produced samples $X_i$ ’s by ${\widehat{P}}_{out}={\sum }_{i=1}^{N}I\left(X_i\right)e^{-X_i∕T}∕{\sum }_{i=1}^N{e}^{-X_i∕T}$ . The parameter $N$ used in simulations is controlled by the relative accuracy ${\sigma }_{p_{out}}∕{\widehat{P}}_{out}$ (Ref. 6). ${\sigma }_{p_{out}}$ is the standard variance of the outage probability and is estimated^{5, 6} by the produced sequence $\left\{X_i\right\}$ . The generation of samples is stopped when the desired accuracy is reached. Step ${\sigma }_X$ and temperature $T$ are important parameters to determine the convergence and accuracy of the result curves. As suggested in Refs. 4, 6, the acceptance rate defined as the number of accepted samples divided by $N$ and the outage rate defined as the number of outage samples divided by $N$ are used to check the validity of the values of ${\sigma }_X$ and $T$ . In each simulation, the acceptance rate^{4, 6} should be around 0.23, and the outage rate⁶ should be around 0.5.

In the simulations, we launch $40\mathrm{Gb}∕\mathrm{s}{2}^7-1\mathrm{PRBS}$ sequences with raised cosine NRZ pulses into the same optical communication system model used in Ref. 2, while the optical fiber channel samples and the BER samples are generated by MCMC as suggested earlier in this letter. The parameter $M$ is 50 in our simulations.

In Fig. 1, the curves of ${\widehat{P}}_{out}$ and the relative accuracy versus the number of fiber channel samples used in simulations are given when the optical communication system is compensated by our proposed higher-order optical PMD compensator² and the mean DGD of the channel normalized to the bit time is 0.3. A fast convergence can be observed in the figure. About $30000$ fiber channel samples are used when the relative accuracy is reduced to 0.1. If MC is used to estimate the outage probability, about ${10}^7$ fiber channel samples should be produced to guarantee the same accuracy.⁶ In fact, MC can estimate $P_{out}$ only down to ${10}^{-4}$ . The superiority of MCMC is obvious compared with MC. All the outage probability values calculated by MCMC in our simulations are determined when the corresponding relative accuracies are approximately reduced to 0.1.

We use MCMC to prove the efficiency of our proposed higher-order optical PMD compensator further also by comparison with the higher-order compensation scheme of Shtaif ^{2, 10} In Fig. 2, the curves of the outage probability versus the mean DGD of the fiber channel normalized to the bit time are given when the input power of the system is $-22.6\mathrm{dBm}$ . Estimated by MCMC, the outage probability values of the system compensated by our proposed higher-order optical compensator are shown as circles, and the outage probability values of the system compensated by the referenced higher-order optical compensator are shown as squares. In order to prove the accuracy of MCMC, some outage probability values of the compensated system estimated using the system configuration samples² produced by MC are also shown in Fig. 2. (Stars: compensated by our proposed optical compensator; crosses: compensated by the referenced optical compensator.) We can see a good agreement between the two simulation tools. MCMC can estimate the outage probability down to ${10}^{-7}$ . In Fig. 2, we can also see that our proposed higher-order compensation scheme always has better performance. As analyzed in Ref. 2, in order to give a fair and creditable comparison, the values of the variable parameters to control the higher-order optical compensators are not optimized. The optimized controlling parameters can ensure the higher-order optical compensators to have better performances.

In this letter, we describe how to use MCMC to accurately estimate the outage probability of the optical communication system compensated by higher-order optical PMD compensators. The efficiency of our proposed higher-order optical PMD compensator is proved further using MCMC as the simulation tool to estimate the outage probability down to extremely low values $\left({10}^{-7}\right)$ . The accuracy and superiority of MCMC are also affirmed.