Polarization Drift Channel Model for Coherent Fibre-Optic Systems

A theoretical framework is introduced to model the dynamical changes of the state of polarization during transmission in coherent fibre-optic systems. The model generalizes the one-dimensional phase noise random walk to higher dimensions, accounting for random polarization drifts, emulating a random walk on the Poincaré sphere, which has been successfully verified using experimental data. The model is described in the Jones, Stokes and real four-dimensional formalisms, and the mapping between them is derived. Such a model will be increasingly important in simulating and optimizing future systems, where polarization-multiplexed transmission and sophisticated digital signal processing will be natural parts. The proposed polarization drift model is the first of its kind as prior work either models polarization drift as a deterministic process or focuses on polarization-mode dispersion in systems where the state of polarization does not affect the receiver performance. We expect the model to be useful in a wide-range of photonics applications where stochastic polarization fluctuation is an issue.


Introduction
Enabled by digital signal processing, coherent detection allows spectrally efficient communication based on quadrature amplitude modulation formats, which carry information in both the intensity and the phase of the optical field, in both polarizations.Polarization-multiplexed quadrature phase-shift keying has recently been introduced for 100 Gb/s transmission per channel and it is expected that in the near future, higher-order modulation formats will become a necessity to reach even higher data rates.However, the demultiplexing of the polarization-multiplexed channels in the receiver requires knowledge of the state of polarization (SOP), which is drifting with time.This implies that the SOP must be tracked by a dynamic equalizer. 1,2 s the SOP changes with time in a random fashion, it is qualitatively different from the chromatic dispersion, which can be compensated for using a static equalizer.A deterministic or static behaviour would be straightforward to resolve, but with a nondeterministic SOP drift, the equalizer must be continuously adjusted.
In order to fully understand the impact of an impairment on the performance of a transmission system, a channel model is required, which should describe the behaviour of the channel as accurately as possible.Based on the statistical information that such models reveal, insights into how to treat the impairments optimally in order to maximize performance can be obtained and used as a result.On the other hand, a channel model that does not describe the fibre accurately may hinder the achievement of optimal performance.Therefore, it is important that the channel model matches the stochastic nature of the fibre closely.
Very few results on modelling of random polarization drifts are present in the literature.In the context of equalizer development, several models have indeed been suggested, but typically by either using a constant randomly chosen SOP [3][4][5][6] or by generating cyclic/quasi-cyclic deterministic changes, [7][8][9][10] which are usually nonuniform in their coverage of the possible SOPs.There is, on the other hand, a wealth of statistical models that describe differential group delay (DGD) and polarization mode dispersion dating back to the eighties. 11However, these results were typically developed to be applicable in systems using intensity modulation or single-polarization (differential) phase-shift keying formats, which are not affected by the SOP drift.Thus, instead of modelling the time evolution of the SOP, differential equations describing the SOP change with frequency and fibre length were typically given. 12,13 here are also some direct measurements of SOP changes, e.g., by Soliman et al. 14 However, in their measurement, fast SOP changes were induced in a dispersion-compensating module under laboratory conditions, without considering the SOP drift of an entire fibre link.It can be concluded that stochastic SOP drift has so far not been given much attention.
This paper suggests a model for random SOP drifts in the time domain by generalizing the one-dimensional (1D) phase noise random walk to a higher dimension.The model is based on a succession of random Jones matrices, where each matrix is parameterized by three random variables, chosen from a zero-mean Gaussian distribution with a variance set by a polarization arXiv:1507.00953v3[physics.optics]4 Nov 2015 linewidth parameter.The latter determines the speed of the drift and depends on the system details.The polarization drift has a random walk behaviour, where each step is independent of the previous steps and equally likely in all directions.The model is given in the Jones, Stokes formalisms and in the more general real 4D formalism. 15,16  argued above, the suggested model serves a different purpose than the models of polarization mode dispersion existing in the literature.In the latter case, the fibre is viewed either as a concatenation of a small number of segments with relatively high DGD, leading to the hinge model, 13,17,18 or the limiting case with segments of infinitesimal lengths, leading to a Maxwellian distribution of the DGD. 19The most common assumption is then to model the hinges as independent random SOP rotators. 13,18  the anisotropic hinge model, the SOP variation at the hinges is modelled as generalized waveplates parameterized by one random angle per hinge. 20The hybrid hinge model is a further generalized way to model the SOP changes at the hinges, 21 but in none of these publications is the SOP time evolution discussed.
The suggested model can be used in simulations for a wide range of photonics applications, where stochastic polarization fluctuation is an issue that needs to be modelled.For example, a fibre-optic communication link can be simulated, independently of the modulated data as well as of other considered impairments, which can be useful to, e.g., characterize receivers' performance.Moreover, it can reveal statistical knowledge of the received samples affected by polarization rotations, based on which the existing tracking algorithms can be optimally tuned or new, more powerful algorithms can be designed.High-precision transfer and remote synchronization of microwave and/or radio frequency signals 22 is another application that could benefit from a better understanding of how it is affected by polarization drifts, which is currently the limiting factor towards a better performance.Other applications such as fibre-optic sensors 23 have been developed for use in a broad range of applications, fibre-optic gyroscopes 24 and quantum key distribution 25 are strongly affected by polarization fluctuations and may benefit from a better understanding of the transmission medium.

Present Phase and Polarization Drift Models
The fibre propagation through a linear medium is often described by a complex 2 × 2 Jones matrix, which, neglecting the nonlinear phenomena, relates the received optical field to the input.Another approach is to use the Stokes-Mueller formalism, where the evolution of the Stokes vectors is modelled by a Mueller 3 × 3 unitary matrix.The latter has the advantage that the Stokes vectors are experimentally observable quantities and can be easily visualized as points on a 3D sphere, called the Poincaré sphere.A further approach to describe the SOP rotations exists in the 4D Euclidean space, 15 where the wave propagation can be described by a 4 × 4 real unitary matrix. 16We will focus most of our discussion on the Jones formalism and connect it to the Stokes and real 4D formalisms later on.
An optical signal has two quadratures in two polarizations and can be described by a Jones vector as a function of time t and the propagation distance z as E(z,t) = (E x (z,t), E y (z,t)) T , where each element combines the real and imaginary parts of the electrical field in the x and y field components and (•) T denotes transposition.The kth discrete-time input into a transmission medium can be written as u k = E(0, kT ), where T is the sampling time and it is commonly related to the symbol interval.The received discrete signal, at distance L, is r k = E(L, kT ).
The propagation of the optical field can be described by a 2 × 2 complex-valued Jones matrix J k , which relates the received optical field r k ∈ C 2 , in the presence of optical amplifier noise and laser phase noise, to the input u k ∈ C 2 as where i = √ −1, φ k models the carrier phase noise and n k ∈ C 2 denotes the additive noise, which is represented by two independent complex circular zero-mean Gaussian random variables.Assuming that polarization-dependent losses are negligible, the channel matrix J k can be modelled as a unitary matrix, which preserves the input power during propagation.In this work, we assume that the chromatic dispersion has been compensated for and polarization mode dispersion is negligible.The transformation J k can be described using the matrix function J(α α α k ), which is expressed using the matrix exponential [26, p. 165]  parameterized by three variables α α α = (α 1 , α 2 , α 3 ) according to 27,28 where σ σ σ = (σ σ σ 1 ,σ σ σ 2 ,σ σ σ 3 ) is a tensor of the Pauli spin matrices 28 This notation of σ σ σ i complies with the definition of the Stokes vector, and it is different from the notation introduced by Frigo. 27he operation α α α • σ σ σ should be interpreted as a scaled linear combination of the three matrices σ σ σ 1 ,σ σ σ 2 ,σ σ σ 3 , and I 2 is the 2 × 2 identity matrix.In general, a 2 × 2 complex unitary matrix has four degrees of freedom (DOFs), but in this case we explicitly factored out the phase noise.Including the identity matrix in σ σ σ would make it possible to account for all four DOFs.The vector α α α can be expressed as a product of two independent random variables α α α = θ a, i.e., its length θ = α α α in the interval [0, π) and the unit vector a = (a 1 , a 2 , a 3 ), which represents its direction on the unit sphere.Since J k is unitary, the inverse can be found by the conjugate transpose operation or by negating α α α k , since The same principle holds for the phase noise (e iφ k ) −1 = (e iφ k ) H = e −iφ k .
Modern coherent systems require a local oscillator at the receiver in order to get access to both phase and amplitude of the electrical field.The local oscillator serves as a reference but it is not synchronized with the transmitter laser, resulting in phase noise, which creates the need for carrier synchronization.The phase noise is modelled by the angle φ , while α 1 , α 2 and α 3 model random fluctuations of the SOP caused by fibre birefringence and coupling.Both phase noise and SOP drift are dynamical processes that change randomly over time.The SOP drift time can vary from microseconds up to seconds, depending on the link type.It is usually much longer than the phase drift time, which is in the microsecond range for modern coherent systems. 29he update rule of the phase noise follows a Wiener process [30][31][32] where φ k is the innovation of the phase noise.An alternative form of equation ( 4) is which we will generalize later to account for the SOP drift.The innovation where σ 2 ν = 2π∆νT and ∆ν is the sum of the linewidths of the transmitter and local oscillator lasers.The accumulated phase noise at time k is the summation of k Gaussian random variables • φ k and the initial phase φ 0 , which becomes a Gaussian-distributed random variable with mean φ 0 and variance kσ 2 ν .The function e iφ k is periodic with period 2π, which means that the phase angle φ k can be limited to the interval [0, 2π) by applying the modulo 2π operation.In this case, the probability density function (pdf) of φ k becomes a wrapped (around the unit circle) Gaussian distribution.It can be straightforwardly verified that the phase drift model has the following properties: 1.The innovation 2. The innovation is symmetric, in the sense that the probabilities for phase changes φ and −φ are equally likely.
3. The most likely next phase state is the current state.5.As k increases, the distribution of e iφ k will approach the uniform distribution on the unit circle.The convergence rate towards this distribution increases with the ∆νT product.

The outcome of two consecutive steps e i
The time evolution of the pdf of a fixed point corrupted by phase noise is exemplified in Fig. 1.
The initial phase difference between the two free running lasers has equal probability for every value, therefore it is common to model φ 0 as a random variable uniformly distributed in the interval [0, 2π), i.e., it has equal probability for every possible state.
The autocorrelation function (ACF) quantifies the level of correlation between two samples of a random process taken at different time instances by taking the expected value E[•] of the product of the samples.The ACF of the phase noise with lT time separation in response to a constant input u is 33 Phase noise pdf evolution.The pdf of e iφ k for φ 0 = π/4 and σ 2 ν = 0.0025 is shown.In (a), k = 1 corresponds to a single innovation and illustrates the second and third properties, i.e., the pdf is symmetric around the current state (the vertical line) and the peak of the pdf is at the current state.In (b), k = 5 and the pdf spreads over the circle.In (c), k = 8000 and the pdf approaches the uniform pdf, which supports the last property.
where the operations |•| and • denote absolute value and Euclidean norm, respectively.Here we neglected the SOP drift given by J k , in order to isolate the effects of the phase noise.We will compute the ACF of the SOP drift below.
In the literature, polarization rotations are generally modelled by the Jones matrix J k in equation ( 1) or subsets of it, obtained by considering only one or two of the three DOFs α 1 , α 2 , α 3 .Contrary to the phase noise, which is a random walk with respect to time, the matrix J k is in previous literature usually kept constant, [3][4][5][6]34 or it follows a deterministic cyclic/quasi-cyclic rotation pattern. The atter is obtained by modelling the parameters of J k as frequency components ωkT .For example, α 3 = ωkT 8, 9, 35 or α 2 , α 3 varied at different frequencies.2,7,10

Proposed Polarization Drift Model
We propose to model the polarization drift as a sequence of random matrices, which exploit all three DOFs of J k .The model simulates a random walk on the Poincaré sphere and we describe it in the Jones, Stokes and real 4D formalisms.In general, 4D unitary matrices have six DOFs, spanning a richer space than the Jones (four DOFs) or Mueller (three DOFs) matrices can.Out of the six DOFs, only four are physically realizable for propagating photons, and the remaining two are impossible to describe in the Jones or Stokes space. 16The Jones formalism can describe any physically possible phenomenon in the optical fibre making it sufficient for wave propagation.The 4D representation is preferred in some digital communication scenarios since the performance of a constellation corrupted by additive noise can be directly quantified in this space in contrast to the Stokes formalism.Even though the extra two DOFs do not model lightwave propagation, they can be used to account for transmitter and/or receiver imperfections, which cannot be done using Jones or Mueller matrices.However, the extra two DOFs are out of the scope of this paper and we will focus on the remaining four.

Jones Space Description
Similarly to the phase noise update equation (5), we model the time evolution of J k as where J( • α α α k ) is a random innovation matrix defined as equation ( 2).This mathematical formulation is not new to the polarization community, as it is commonly used to describe the polarization evolution in space (each Jones matrix represents a waveplate).However, here it models the temporal evolution of the SOP, where each innovation matrix corresponds to a time increment.The parameters of the innovation J( • α α α k ) are random and drawn independently from a zero-mean Gaussian distribution at each time instance k where σ 2 p = 2π∆p T , and we refer to ∆p as the polarization linewidth, which quantifies the speed of the SOP drift, analogous to the linewidth describing the phase noise, cf.equation (6).Drawing α α α from a zero-mean Gaussian distribution results in special cases of θ and a, where the former becomes a Maxwell-Boltzmann distributed random variable, and the vector a is uniformly distributed over the 3D unit sphere, implying that the marginal distribution of each a i is uniform in [−1, 1]. 36t is important to note that, contrary to phase noise, the equivalent vector α α α k parameterizing J k in equation ( 8) does not follow a Wiener process, i.e., α α α k = • α α α k +α α α k−1 , because in general J(α α α 1 )J(α α α 2 ) = J(α α α 1 +α α α 2 ).Equality occurs for (anti-)parallel α α α 1 and α α α 2 , and holds approximately when α α α 1 1 and α α α 2 1.In a prestudy for this work, 37 we incorrectly used a Wiener process model for the polarization drift.We will return to that model later and discuss its shortcomings.

Stokes Space Description
The evolution of the SOP is often analysed in the Stokes space, where the Jones vectors are expressed as real 3D Stokes vectors and can be easily visualized on the Poincaré sphere.The transmitted Jones vector u k can be expressed as a Stokes vector [38, eq.(2.5.26)] where the ith component of s u k is given by u H k σ σ σ i u k and (•) * denotes conjugation.The equivalent Stokes propagation model of equation ( 1) can be written It is important to note that only s r k and s u k can be obtained by applying equation ( 10) to r k and u k , respectively.The noise component s n k consists of three terms where the first two represent the signal-noise interaction and the last one the noise source.It should be noted that there is no time averaging in equation (10) and it represents an instantaneous mapping to the Stokes space.Thus, the noise terms are polarized and rapidly varying.
The matrix M k is a 3 × 3 Mueller matrix, corresponding to the Jones matrix J k , and the polarization transformation introduced by it can be seen as a rotation of the Poincaré sphere.It has three parameters, the same as J k , and can be written as M k = M(α α α k ), where the function M(•) can be expressed using the matrix exponential 28 M(α α α) = exp(2K(α α α)) = exp(2θ K(a)) where θ = α α α , a = α α α/θ and K(a) denotes The inverse can be obtained by The transformation in equation ( 13) can be viewed in the axis-angle rotation description as a rotation around the unit vector a by an angle 2θ .
Note that any operation that applies a phase rotation on both polarizations with the same amount, such as phase noise or frequency offsets, cannot be modelled in the Stokes space.This can be seen in equation (11), where the phase noise does not alter the transmitted Stokes vector directly, but only contributes to the additive noise in equation (12), which would not exist in the absence of n k .
Analogously to equation ( 8), the time evolution of M k is where M( • α α α k ) is the innovation matrix parameterized by the random vector • α α α k defined in equation ( 9).Fig. 2 shows an example of an SOP drift as it evolves with time.The line represents the evolution of the vector M k (0, 0, 1) T for k = 1, . . ., 3000.

4D Real Space Description
[41] In this case, the transmitted/received sample and the noise term in equation ( 1) can all be represented as a real four-component vector (Re(E x ), Im(E x ), Re(E y ), Im(E y )) T .The transformations induced by the channel are modelled by 4 × 4 real unitary matrices.The e iφ k J k transformation in equation ( 1) can be combined into and the function R(•) can be described using the matrix exponential of a linear combination of four basis matrices 16 where ρ ρ ρ = (ρ ρ ρ 1 ,ρ ρ ρ 2 ,ρ ρ ρ 3 ) and λ λ λ 1 are four constant basis matrices [16, eqs.( 20)-( 23)] The inverse can be obtained as The update of R k in equation ( 16) can be expressed analogously to equations ( 8) and (15) as where the phase innovation

Polarization Drift Model Properties
Due to the similarities between the phase noise and the SOP drift, we will use the properties of the phase noise previously presented as guidelines to validate the proposed channel model.These properties can be reformulated as follows: 1.The innovation 2. The innovation is isotropic, in the sense that all possible orientations of the changes of the Stokes vector corresponding to a movement given by one innovation are equally likely.9) are T = 2.2 h (set by the measurements) and ∆p = 60 nHz (obtained by fitting the dash-dotted ACF line in Fig. 4).In (a) and (d), k = 2 innovation steps are plotted, whereas k = 8 in (b), (e) and k = 16 in (c), (f).Gaussian-like isotropic distributions can be noted in all cases, simulations and measurements, leading to a good (visual) agreement.The spread over the sphere increases with k and the pdf will become uniform if we let k grow large enough.Unfortunately, our measured data do not cover a long enough time period such that uniformity is achieved.
3. The most likely next SOP is the current state.

The outcome of two consecutive steps M(
5. As k increases, the pdf of the product M k s u , for a constant s u , will approach the uniform distribution on the Poincaré sphere.The convergence rate towards this distribution increases with the ∆p T product.
In the following, we will investigate these properties and discuss their validity.

Independent Innovations
This can be easily concluded from the updating rule in equation ( 8), (15) or (19), where the parameters of the innovation do not depend on neither the previous innovations nor the state.

Isotropic Innovation
We will use the following theorem to show that the innovation M( • α α α) is isotropic.Theorem 1.Let a random unit vector a ∈ R 3 be uniformly distributed over the 3D sphere, γ be a random angle with an arbitrary pdf and x ∈ R 3 an arbitrary unit vector.The pdf of the vector y = M(γa)x, where M(•) is defined in equation (13), is invariant to rotations around x, i.e., M(β x)y has the same pdf regardless of β .
In other words, the theorem states that the pdf of the random vector y = M(γa)x does not change if it is rotated by any angle around x.This can be true only if y is isotropic (centred and equally likely in all directions) around x.The technical details of the proof are presented in Supplementary Section I.
Note that Theorem 1 holds for our proposed Stokes innovation matrix M( • α α α) since • α α α = θ a and a is uniformly distributed over the 3D sphere, hence the vector M( • α α α)s u is isotropic around s u .The evolution of M k s u can be seen as an isotropic random walk on the Poincaré sphere starting at M 0 s u and taking random steps, of size proportional to σ p , equally likely in all directions.Curiously, the isotropic property is given only by the fact that the rotation axis a is uniformly distributed over the sphere, while the angle θ may have any pdf.In fact, the pdf of the angle θ determines the shape of the pdf of M( • α α α)s u , which we will discuss later.In contrast, our preliminary SOP drift model 37 does not fulfil the isotropicity condition.The update method of the channel matrix was done by modelling α α α k as Wiener processes, which does not fulfil Theorem 1.

Pdf of a Random Step
Unfortunately, we were not able to find a closed form expression of the pdf of the point M( • α α α)s u for a fixed s u and a random M( • α α α).Using approximations, valid for σ 2 p 1, which is the case for most practical scenarios, the pdf of M( • α α α)s u can be approximated by a bivariate Gaussian pdf centred at s u of variance σ 2 p I 2 on the plane normal to s u .The peak of the pdf is located at s u and we can conclude that the next most likely SOP is the current one.The details of the derivations are provided in Supplementary Section II.In the same section, under the same assumption of small σ 2 , we show that the outcome of two consecutive random innovations can be achieved by a single random innovation if the variance of the random variables • α α is doubled, which fulfils property Uniformity point M u an random on Poincaré therefore the of k the of sphere approach meaning all will equally This intuitive taking steps all with preferred will to uniform This was by et 42

measurements a cable in converge of Poincaré
We compare the model with experimental results in 3, where the evolution of the of a Stokes vector affected by polarization drift after numbers of innovations starting from a point is shown.The experimental data was obtained measuring the Jones matrices of a km long buried fibre from 1505 nm to 1565 nm in steps of 0.1 nm for 36 days at every ∼ 2.2 h.The technicalities of the measurement setup and postprocessing have been published elsewhere. 43The histograms corresponding to the measurements were captured from all the Stokes vectors obtained by applying equation (10) to the product of the matrix J t (ω)J −1 t−k (ω), i.e., the measured matrix corresponding to k innovation steps, with a constant Jones vector, where J t (ω) denotes the measured Jones matrix at time t and wavelength ω.

Initial State
The initial channel matrix M 0 should be chosen such that all the SOPs on the Poincaré sphere are uniformly distributed, i.e., equally likely, after the initial step M 0 s u , for any s u .Analogously, the initial phase φ 0 in equation ( 4) should be chosen from a uniform distribution in the interval [0, 2π).In order to generate such a matrix M 0 , the axis a must be uniformly distributed over the unit sphere and the distribution of the angle θ ∈ [0, π/2) must be (1 − cos θ )/π. 44The generation of the angle θ is not straightforward, and therefore we present an alternative to generate the axis and the angle simultaneously. 41The vector α α α 0 = θ a is formed from the unit vector (cos θ , a 1 sin θ , a 2 sin θ , a 3 sin θ ) T = g/ g , where g ∼ N (0, I 4 ), which will satisfy the conditions for both axis a and angle θ .This approach of generating α α α 0 of the initial channel matrix can be used regardless of the considered method to characterize the polarization evolution, i.e., Jones , Stokes or real 4D formalism.

The SOP Autocorrelation Functions
The ACF is commonly used to quantify how quickly, on average, the channel changes in time, frequency, etc. [45][46][47] The ACF of the SOP drift separated by a time difference of lT in response to a constant input u can be expressed as where r k = J k u is the received sample.The details of this derivation can be found in Supplementary Section III.Using the approximation (1 − |l|σ ACF of phase noise eq. ( 7) ACF of SOP drift in the Jones/4d space eq.( 20) ACF of SOP drift in the Stokes space eq.( 22) experimental data Solid/dashed lines refer to the analytic expressions, whereas the triangles are extracted from measurements.We observe excellent agreement between the experiment and theory, except in the tails of the experimental ACF.This inconsistency can be caused by the lack of accuracy in that region.
The ACF expressions in equations ( 20) and ( 21) hold for the Jones and real 4D space descriptions.Using an analogous derivation as for equation (20), it can be shown that the ACF of the SOP drift in the Stokes space is where s r k = M k s u is the Stokes received sample of a constant input s u affected by polarization drift.Using the following approximations: exp(−2σ p ≈ − p σ p 0 (1 8π|l|T |l| exp(−8π|l|T ∆p) for σ 2 p 1, can be approximated as Both ACFs of the polarization drift depend only on the time separation l but not on the absolute time k.Comparing equation ( 21) with (23), it is interesting to note that even if they describe the same physical phenomenon, on average, a small movement in the Jones/real 4D space will result in a movement that is 8/3 larger in the Stokes space.This result was also previously observed in similar setups analysing polarization-mode dispersion, where the autocorrelation, with respect to frequency, was derived. 43,45 he ACF of the phase noise equation ( 7) has the slowest decay rate, being a factor of three slower than equation (21), and a factor of eight slower than equation (23).Fig. 4 shows the ACFs of the phase noise and the SOP drift, and the later is compared with experimental results.Here the analytic equations ( 20) and ( 22) match the experiment for ∆p = 80 nHz and ∆p = 60 nHz, respectively, and T = 2.2 h.We attribute the discrepancy between the two values to the nonunitary of the measured Jones matrices, which, through the nonlinear operation in equation (10), cause the inconsistency.

Linewidth Parameters
The choice of ∆ν and ∆p is important when a system is simulated considering phase noise and/or SOP drifts.Both parameters reflect physical properties of the system.The quality of the (transmitter and receiver) lasers can be quantified by the ∆ν parameter, which represents the sum of the linewidths of the deployed lasers.Distributed feedback lasers are commonly employed in transmission systems due to their cost efficiency.The linewidth of such lasers varies from 100 kHz to 10 MHz. 30 The polarization linewidth parameter depends on the installation details.Measurements of polarization fluctuations have been reported varying from weeks (this work) to seconds, 48 milliseconds 49,50 or microseconds under mechanical perturbations. 51he polarization linewidth could be quantified through measurements of the ACF, either in the Jones or Stokes space, as in Fig. 4.

Summary
This section is intended to summarize the previous sections and provide an easy implementable form of the proposed model without requiring knowledge about the details of the derivations.
First, the parameters of the channel must be selected: • The sampling time T , often equal to the symbol interval.
• The laser linewidth ∆ν, where the contributions of the transmitter and receiver lasers are taken into account.
The model is then implemented as follows: 0. Set the initial state of the channel: • φ 0 ∼ U {0, 2π}.
2. Apply the model for every transmitted sample u k , which results in the received sample r k : • r k = e iφ k J k u k +n k , where n k denotes the additive noise represented by two independent complex circular zero-mean Gaussian random variables.
Repeat the last two steps for every sample u k .Alternatively, the Jones formalism used above can be replaced by the Stokes formalism using equation ( 13) instead of equation ( 2) and neglecting the phase noise e iφ ; or 4D formalism by combining e iφ k and J k into R k given by equation (17).In the former case, u k and r k are 3D Stokes vectors defined as equation (10), whereas in the latter case, u k , r k and n k are 4D real vectors.
In the summary above, we have included phase and additive noise for completeness.Without loss of generality, the model can be applied with/without phase and additive noise; also other impairments can be considered.

Conclusions
We have proposed a channel model to emulate random polarization fluctuations based on a sequence of random matrices.The model is presented in the three formalisms, Jones, Stokes and real 4D, and it is parameterized by a single variable, called the polarization linewidth.
The model has an isotropic behaviour, which has been successfully verified using experimental data, where every step on the Poincaré sphere is equally likely in all directions emulating an isotropic random walk and can be easily coupled with any other impairments to form a complete channel model.
Compared to the existing literature, the fundamental advantages of the proposed model are randomness and statistical uniformity.Such a model is relevant for a wide range of fibre-optical applications where stochastic polarization fluctuations are an issue.It can potentially lead to improved signal processing that accounts optimally for this impairment and more realistic simulations can be carried out in order to accurately quantify system performance.

II 3D Distribution Analysis
In this section, we derive the approximate pdf of the point s r = M( • α α α)s u for a fixed s u and random M( • α α α).Exact expressions are very difficult to obtain, therefore we make use of the approximations sin 2θ ≈ 2θ and cos 2θ ≈ 1, valid for σ 2 T and it can be noted that s r then has a bivariate Gaussian distribution on the plane normal to s u and the peak of the distribution centred at s u .
Using equation (S8) and by removing high order terms, such as α i α j for any i, j, the multiplication of two matrices M(α α α) can be approximated as From equation (S9) we can conclude that two consecutive small innovations can be replaced by a single innovation, i.

III Autocorrelation
In this section, we derive the ACF of the SOP drift in equation (20).The derivation uses the Jones description of the model but the result is valid for the 4D description as well.
At first we will calculate the expectation of the innovation matrix from equation ( 2)

• φ 1
e i • φ 2 can be emulated by a single step e i • φ t by doubling the variance of • φ k .

3 fFigure 3 .
Figure 3.The histograms of M k s u for different steps k and a fixed s u = [1, 1, 1] T / √ 3 obtained from the model (top row) and from measurements (bottom row) are shown.The highest density is represented by dark red and the lowest by dark blue, the outer part of the density.The parameters of the simulated drift in equation (9) are T = 2.2 h (set by the measurements) and ∆p = 60 nHz (obtained by fitting the dash-dotted ACF line in Fig.4).In (a) and (d), k = 2 innovation steps are plotted, whereas k = 8 in (b), (e) and k = 16 in (c), (f).Gaussian-like isotropic distributions can be noted in all cases, simulations and measurements, leading to a good (visual) agreement.The spread over the sphere increases with k and the pdf will become uniform if we let k grow large enough.Unfortunately, our measured data do not cover a long enough time period such that uniformity is achieved.

Figure 4 .
Figure 4. ACF comparison.The normalized ACF of the phase noise and SOP drift is plotted versus normalized time.Solid/dashed lines refer to the analytic expressions, whereas the triangles are extracted from measurements.We observe excellent agreement between the experiment and theory, except in the tails of the experimental ACF.This inconsistency can be caused by the lack of accuracy in that region.