Jointcharacteristic Function of the First and SecondorderPolarizationmodedispersion Vectors in Linearly Birefringent Optical Fibers
 Author: Lee Jae Seung
 Organization: Lee Jae Seung
 Publish: Journal of the Optical Society of Korea Volume 14, Issue3, p228~234, 25 Sep 2010

ABSTRACT
This paper presents the joint characteristic function of the first and secondorder polarizationmodedispersion (PMD) vectors in installed optical fibers that are almost linearly birefringent. The joint characteristic function is a Fourier transform of the joint probability density function of these PMD vectors. We regard the random fiber birefringence components as white Gaussian processes and use a FokkerPlanck method. In the limit of a large transmission distance, our joint characteristic function agrees with the previous joint characteristic function obtained for highly birefringent fibers. However, their differences can be noticeable for practical transmission distances.

KEYWORD
Polarization mode dispersion , Optical fiber , Optical fiber transmission , Optical communication

I. INTRODUCTION
Random fluctuations of fiber birefringences cause polarization mode dispersion (PMD) effects detrimental to high bitrate channels. The PMD effects can be described using a PMD vector whose magnitude is equal to the differential group delay between two principal states [1], [2]. When there are many optical channels or the channel bit rate is high, the angular frequency derivative of the PMD vector becomes also important [3]. Often, the conventional PMD vector is called the firstorder PMD vector while its angular frequency derivative is called the secondorder PMD vector.
In [4], Foschini and Shepp derived a joint characteristic function for two white Gaussian vector processes in a compact functional form using a sinecosine Fourier expansion representation. The joint characteristic function is a Fourier transform of the joint probability density function (pdf) of two vector processes. Shortly after [4], it was shown that, when the fiber transmission distance is large, the joint characteristic function can be used for the first and secondorder PMD vectors in highly birefringent fibers [5]. This joint characteristic function has been used extensively to describe various distribution characteristics of PMD vectors [68]. The birefringence vector model used in [5] assumes a fixed linear birefringence with much smaller linear and circular random birefringences in equal amounts. Interestingly, the joint characteristic function of [5] agrees with the result of [9] that uses the model of [5] without the fixed linear birefringence component. However, installed optical fibers are almost linearly birefringent [10]. Thus it is worthwhile to justify the joint characteristic function of [5] for installed fibers.
Regarding the birefringence components as white Gaussian processes, a FokkerPlanck equation [11] was derived in [12] to investigate the pdf behavior of the firstorder PMD vector in linearly birefringent fibers. The FokkerPlanck method was used also in [13] to find a more correct joint characteristic function for the birefringence vector model of [9]. The secondorder angular frequency derivative of the fiber birefringence vector was included in [13] but was neglected in [5] and [9]. The birefringence vector distribution of [9] and [13], having linear and circular random birefringences in equal amounts, has a spherical symmetry in Poincare space. This symmetry helps to find the joint characteristic functions of [9] and [13]that are dependent only on the magnitude of two PMD vectors and the angle in between in their Fourier domain. However, the spherical symmetry does not hold for installed optical fibers and, to our knowledge, the joint characteristic function for installed optical fibers is still unknown.
In this paper, we find the joint characteristic function for installed optical fibers. We use the FokkerPlanck method and include the secondorder angular frequency derivative of the fiber birefringence vector. Our final result obtained after a complicated procedure has a compact functional form similar to the previous joint characteristic functions [5, 9, 13].
At first, we will present details of the FokkerPlanck method of [12] along with a new shortcut procedure for finding the pdf of the firstorder PMD vector. Then we extend the approach to find a FokkerPlanck equation (in Fourier domain) for the first and secondorder PMD vectors. The asymptotic solution of our FokkerPlanck equation yields the joint characteristic function.
II. PDF FOR THE FIRSTORDER PMD VECTOR
The dynamical PMD equation,
describes the evolution of the firstorder PMD vector, τ =(τ 1, τ 2, τ 3) along the fiber [1]. L is the fiber transmission distance.
β is the birefringence vector andβ Ω=∂β /∂Ω, where Ω is the angular frequency. For linearly birefringent fibers, we setwhere
is a zeromean white Gaussian process having a correlation property,
denotes the ensemble average and δij is the Kronecker delta. This modeling neglects the spatial size of the correlation length compared with the large parameter
L [9], [12], [13]. The proportionality factora is a function of Ω and we will denoteda /d Ω asa Ω. Similarly,d β /d Ω is simplified toβ Ω=a Ωetc.
We note that the vector equation (1) is composed of three Langevin equations,
where
N (=3) is the number of variables. The matrix formed byg _{ij} isThen we find the FokkerPlanck equation as follows [11]:
where
D _{i} andD _{ij} are drift and diffusion coefficients, respectively, defined asFrom these relations, we have
D _{1}=a^{2}τ_{1},D _{2} =a^{2}τ_{2},D _{3}=a^{2}τ _{3}, andP =P (τ,L ) is the pdf for τ . Since the input signal has experienced no PMD degradations, the initial condition is a Dirac delta function,P (τ ,0) =δ (τ).The Fokker Planck equation for the firstorder PMD vector is given by
where
As is shown in Fig. 1, we use a spherical coordinate using the following transform relations [14]:
The positive τ 3axis is the polar axis. θ and φ are polar and azimuth angles, respectively. Since our problem is invariant to the rotation about the τ 3axis, we may set ∂/∂φ=0. Then the differential operator
K can be simplified greatly asThis equation implies that (5) is a diffusion equation along τ 1, τ 2, and θ directions. Therefore, in the asymptotic region, where the transmission distance
L is large,P =P (τ ,L ) becomes smooth and widely spread in τ space.We expand
P (τ ,L ) in terms of Legendre polynomials asThe nth order Legendre polynomial,
P _{n}(cosθ), is an eigenmode of thea ^{2} term such thatSince
P_{n} (cosθ) is not the eigenmode ofa _{ω}^{2} terms of (10), there are couplings between these modes during the transient region. The initial condition, P(τ , 0)=δ(τ ), belongs to theP _{0}(cosθ) mode andA _{n}(τ , 0) forn ≠0. AsL increases from zero, the magnitude ofA _{n}(τ ,L ) (n ≠0) increases owing to the coupling with then =0 mode. Actually, these couplings exist only between even modes owing to the even symmetry of our problem with respect to the θ coordinate. Note that thea ^{2} term in (10) is dependent only on the angular variable θ while thea _{ω}^{2} terms are not. Thus, in the asymptotic region, whereP =P (τ ,L ) becomes smooth, thea ^{2} term becomes dominant compared with thea _{ω}^{2} terms forn ≠0 modes. Ifa _{ω}^{2} terms are neglected,A _{n}(τ ,L ) decays in proportion to exp{n (n +1)a ^{2}L }. It implies that, in the asymptotic region, we haveP (τ ,L )=A _{0}(τ ,L ) ultimately. This can be verified also using a multiplescale method [15].To find
A _{0}(τ ,L ) in the asymptotic region, we set ？/？θ= ？/？φ=0 for the spherical coordinate representations of ？^{2}/？ τ _{1} ^{2} and ？^{2}/？τ _{2} ^{2} terms inK , which givesNext, we remove the terms in
K that coupleA _{0}(τ ) with otherA _{n}(τ ,L ) (n ≠0) components. This can be achieved by applying (1/ 2)f ^{π}_{0}d θ sin θonK, which leaves only theP _{0}(cosθ) component inK . The result isThus, in the asymptotic region,
P (τ ,L ) becomes a Gaussian, expwith a variance σ^{2}(
L )=4a _{ω}^{2}L /3. This also implies that the magnitude of the first order PMD vector's pdf converges to a Maxwellian [2].We have used the same initial condition,
P (τ, 0)=δ(τ), to solve (13). This can be verified by taking a Laplace transform of (5) to findwhere
is the Laplace transform of
with respect to
L . The initial condition term, δ(τ), also belongs to theP _{0}(cosθ) mode and remains the same in the asymptotic region. Thus the initial condition,P (τ, 0)=δ(τ), can be used for (13).III. JOINT CHARACTERISTIC FUNCTION FOR THE FIRSTAND SECONDORDER PMD VECTORS
From (1), we find a differential equation for the secondorder PMD vector, τω=？τ/？ω=(τ_{ωl}, τ_{ω2}, τ_{ω3}), as follows:
where
β _{ωω}=？^{2}β /？ω^{2}=a_{ωω}(γ_{1}, γ_{1}, 0). We note that (1) and (14) can be regarded as six Langevin equations, ？ ？ =where
N =6.x _{i}=τi fori = 1, 2, 3 andx _{i}=τ_{ωi} fori = 4, 5, 6. Nowg _{ij} set forms a 6×6 matrix that will be used to find new drift and diffusion coefficients in the FokkerPlanck equation, (3). Its solution is the joint pdf of τ and τ_{ω} which will be denoted asP =P (τ, τ_{ω},L ).We will solve the FokkerPlanck equation in the Fourier domain by introducing a joint characteristic function,
where
k =(k _{1},k _{2},k _{3}) andk _{ω}=(k _{ω1},k _{ω2},k _{ω3}). Our definition of the join]t characteristic function is the complex conjugate of the conventional one. The differential equation for the joint characteristic function is given aswhere the differential operator
M is given byIn (18), we have introduced angular momentum operators defined as
L _{k} =jk ×∇_{k} = (L_{k1}, L_{k2}, L_{k3}), L_{kω}= jk ×∇_{kω}=(L _{kω1},L _{kω 2},L _{kω 3}), andL _{tot}=L _{k}+L _{kω}=(L _{tot1},L _{tot2},L _{tot3}) [16]. ∇_{k} and ∇_{kω} are del operators ink andk _{ω} domains, respectively.L _{tot}=L _{tot1}+L _{tot2} andL _{tot}=L _{tot1}L _{tot2} are raising and lowering operators, respectively, for the total angular momentum operator,L _{tot}. Since the input signal has experienced no PMD degradations, the initial condition isQ (k ,k _{ω}, 0) = 1 orP (τ ,τ_{ω}, 0)=δ(τ)δ(τ_{ω}).As is shown in the Appendix A, we can find
E {τ^{2}}=4a ^{2}_{ω}L from (16) using the relation,E {τ^{2}}=∇^{2}_{k}Q (k, k_{ω},L )'_{k=kω=0} . In a similar way,E {τ^{2}_{ω}} is found as 16a ^{4}_{ω}L ^{2} /3+4a ^{2} _{ωω}L (1a ^{4}_{ω} /9a ^{2} _{ωω}a ^{2})+2a ^{4}_{ω} {1exp (6a ^{2}L )}/ 27a ^{4}. AsL becomes large, the ratioE {τ^{2}_{ω}}/(E {τ^{2}})^{2} converges to 1/3 [5]. In this asymptotic region,Q(k, k _{ω},L ) becomes localized to the origink=k _{ω}=0.Since
E {τ^{2}_{ω}} >>E {τ^{2}},Q (k, k, L ) shrinks faster ink _{ω} domain than ink domain as L increases.The
L ^{2} _{tot} operator has eigenvalueL _{tot}(L _{tot} +1), whereL _{tot} is a nonnegative integer. We may expandQ(k, k_{ω}, L) using the eigenmodes ofL ^{2}_{tot}. TheL ^{2}_{tot} operator is composed of angular variables only and becomes dominant in the asymptotic region. The other terms inM become less effective. This is because, as is mentioned above,Q(k, k_{ω}, L) shrinks to the origin k=k_{ω}=0 and the shrink speed is faster ink _{ω} domain. Thus eigenmodes withL _{tot} ≠ 0 become negligible in the asymptotic region and only the eigenmode withL _{tot} = 0 survives. Note that this procedure is similar to the foregoing firstorder PMD vector analysis. TheL _{tot} = 0 condition also implies that the only possible eigenvalue ofL _{tot3} is zero since L _{tot}≤L _{tot3} ≤L _{tot}. Thus we haveL _{tot±}Q =L _{tot3}Q = 0 orFrom (24), we find
kLQ _{kω} =k _{ω}L _{k}Q = 0. Using these relations, it can be shown thatD andF terms in the operatorM disappear, i.e.,DQ =FQ = 0 . As a result, there are no differential terms fork _{ω} in the operatorM andk _{ω} becomes just a parameter in the asymptotic region. The symmetry in (24) implies that the asymptotic solution should be dependent only on the magnitudes ofk andk _{ω} vectors and the angle in between. This can be proved simply by rotating the coordinate such that one ofk ork _{ω} vector becomes a polar axis.To find the asymptotic solution satisfying this condition, let's assume that the
k _{ω} vector is located in thek vector domain as is shown in Fig. 2. Since our problem is invariant under the rotation about thek _{3}axis, we have setk _{ω}2. We rotate the (k_{1}, k_{2}, k_{3} ) coordinate about thek _{2} axis into the (k_{1}', k_{2}', k_{3}' ) coordinate, where thek _{3}' axis coincides with thek _{ω} vector's direction.The corresponding rotation matrix is
where θ_{ωk} is the polar angle of
k _{ω} before the rotation. With (25), we apply the transform relationsand
to
M , whereR _{ji} is the matrix element of the inverse rotation matrixR (θ_{kω}). Then the remaining operators,B, C , andE , are transformed intoIn the transformed (
k_{1}', k_{2}', k_{3}' ) coordinate, the asymptotic solution is not dependent on the azimuth angle φ′. Thus we set ？/？φ'=k_{1}'？/k_{2}'k_{2}'？/？k_{1}=0 inM . In addition, we applyon
M to remove θ_{kω} and φ′ dependences inM . This operation removes the couplings between the L _{tot} = 0 mode and the other modes with L _{tot} ≠ 0. Then cos^{2} θ_{ωk}, sin^{2} θ_{ωk}, and cosθ_{ωk}sinθ_{ωk}, factors become 1/3, 2/3, and 0, respectively. Also, cos^{2}φ', sin^{2}φ', and cosφ'sinφ' factors become 1/2, 1/2, and 0, respectively. Besides, we find ？/？k _{1}'^{2} = ？/？k _{2}'^{2}.Finally, we obtain the following differential equation for the asymptotic solution:
where we have omitted the prime notation. As before, the initial condition remains the same,
Q(k, k_{ω}, 0) = 1. We decomposeQ(k, k_{ω}, L) aswhere
q satisfiesEquation (31) can be solved using the separationofvariables method as
where
H _{2m}(^{.}) is the Hermite polynomial. Using an expansion formula for a Gaussian function,we find
Thus, for an arbitrary direction of
k _{ω}, the joint characteristic function becomeswhere σ^{2} =
E {τ_{1}^{2}}=E {τ_{2}^{2}}=E {τ_{3}^{2}}=4a ^{2}_{ω}L /3 as in the previous section.k _{//} is the component ofk parallel tok _{ω} andk ^{2}_{1}=k ^{2}k ^{2}_{//}.In the limit of large
L , the cosh (k _{ω}σ^{2}) term in (35) implies that the effective range ofk _{ω} decreases asO (L ^{1}) since σ^{2} ∝L . Also, the exponential terms in (35) imply that the effective range ofk _{i}(i =1, 2, 3) decreases asO (L ^{1/2}). Thus we may neglect thea _{ωω} term in (35) asL increases indefinitely. In this limit, (35) becomes the same as the joint characteristic function of [5] with reevaluated appropriate to the linearly birefringent fiber. This feature also holds when the fiber has both linear and circular birefringences in equal amounts [9], [13]. Thea _{ωω} term can be meaningful for practical transmission distances. For example, the secondorder PMD vector components can be Gaussianlike distributed instead of a soliton shape [2]. This point has been shown in [13] and presented in Appendix B.With the
condition, the firstorder PMD vector distribution has an azimuthal symmetry for all L in [12]. However, the azimuthal symmetry is not evident in our analysis of jointcharacteristic function where the number of variables is doubled from 3 to 6. In fact, our final asymptotic solution (35) has spherical symmetry. The coordinate orientation does not affect our result. This kind of spherical symmetry for the asymptotic jointcharacteristic function is also found in [5], [9], and [13].
It is also interesting that (35) is quite similar to the result of [13]. The only difference is the change of
L into 2L /3 in our case. In [13], the birefringence vector was assumed asis also a zeromean white Gaussian process having a correlation property,
Note that the number of random birefringence components is larger than our case. Thus, for given
a ,a _{ω}, anda _{ωω} values, the joint characteristic function reaches to its asymptotic form faster than our case as the transmission distance increases. In this respect, the scaling factor 2/3 can be regarded as the ratio of the birefringence component numbers.IV. CONCLUSION
We have found a joint characteristic function for the first and secondorder PMD vectors in installed fibers. As the fiber transmission distance increases indefinitely, our joint characteristic function resembles that of [5] obtained for highly birefringent fibers. This agreement may not hold for practical transmission distances owing to the secondorder angular frequency derivative of the fiber birefringence vector included in our analysis. As an example, we have shown that the secondorder PMD vector components can be Gaussianlike distributed instead of a soliton shape. If the fiber length is scaled by a factor of 2/3, our result obtained for linearly birefringent fibers has the same functional form with the previous joint characteristic function of [13] obtained for fibers having a spherically symmetric birefringence vector distribution in Poincare space.

[FIG. 1.] Spherical coordinate for τ space.

[FIG. 2.] Rotation of k vector space coordinates about the k2axis.