Borrowing methodology from quantum mechanics, we propose and develop a density-matrix formalism for modal coupling and dispersion in mode-division multiplexing communications systems. The central concept in our formalism is the density matrix, from which all observable information of an optical field can be handily accessed. In the formalism, we derive fundamental evolution equations and concatenation rules for the key elements that characterize essential modal properties, and construct a statistical model ready for the numerical analysis of stochastic light propagation in randomly perturbed fibers. Unlike the Stokes-vector formalism that requires J2 − 1 auxiliary Gell-Mann matrices, the density-matrix formalism can be directly formulated for arbitrary modal‐space dimension J. Based on the density-matrix formalism, the statistical modal properties of a 4-mode fiber under random perturbation are numerically investigated, which raises an interesting possibility of optimizing the modal dispersion by manipulation of the random perturbation.
1. Introduction
Modal coupling and dispersion in fibers is a fundamental problem of mode-division multiplexing (MDM) technology [1–4] in modern optical communications. Even a “single-mode” fiber actually supports two orthogonally-polarized modes that might be mutually coupled and differ in group delays under perturbation of the fiber, and the highly successful theory to analyze polarization-mode dispersion (PMD) in such $2$-dimensional $(J=2)$ modal space is the Stokes formalism [5–8] . As current MDM technology exploits more orthogonal optical modes for signal transmission, study on coupling and dispersion of multiple ($J>2$) modes becomes increasingly crucial to high-capacity and high-performance communications networks [9,10]. For example, in the OAM-based multiplexing [11–15] that serves as a promising candidate of MDM, optical modes (HE$_{\nu p}$, EH$_{\nu p}$, $\nu =1,2,\ldots ,p=1,2,\ldots$) of well-defined orbital angular momentum (OAM) [16–18] are employed to encode and transmit signals. As the fiber is perturbed, the even and odd modes of the same HE$_{\nu p}$ or EH$_{\nu p}$ could substantially couple to each other, and HE$_{\nu +1,p}$ and EH$_{\nu -1,p}$ modes of the same $\nu$ might also be significantly coupled in weakly-guided fibers. In such cases, one needs to treat the coupling and dispersion of at least $J=4$ orthogonal modes, and the Stokes formalism has been extended to modal spaces of arbitrary dimension to investigate modal coupling and dispersion in fibers supporting multiple modes [19–22] .
At heart of the Stokes formalism lies the Stokes vector that describes the state of electromagnetic fields. The essential necessity of the Stokes vector is for the description of incoherent (or mixed) field states, which cannot be described by the Jones vector. Moreover, in the PMD theory, the Stokes space corresponding to the $2$-dimensional complex modal space (Jones space) is a $3$-dimensional real vector space, such that the Stokes vector can be illustratively viewed as a geometric vector in the Poincaré Sphere, which gives the Stokes formalism additional powerful advantage in the analysis of PMD. However, for a modal space of higher dimension ($J>2$), the Stokes vector (of dimension $J^{2}-1>3$) can no longer be treated geometrically, and its appliction requires $J^{2}-1$ auxiliary Gell-Mann matrices [19,20]. Therefore, it should be appealing to devise a theoretical tool that can properly describe general (i.e., both pure and mixed) field states and is simple to apply in modal spaces of any dimension $J$.
In this paper, borrowing methodology from quantum mechanics, we introduce such a theoretical tool named density matrix operator $\hat {\rho }$ for the description of general field states, and develop a formalism that can be conveniently applied to analyze modal coupling and dispersion of multiple-mode fields. Although notations similar to the density matrix operator for pure states occasionally appeared in literature, they were mostly invoked as an auxiliary (without explicit physical meaning) in certain intermediate steps of derivations [6,19,20]. Contrastingly, the density matrix operator in its most general form plays a central role in our proposed formalism. From the most fundamental propagation rules of electromagnetic fields in fiber, we construct the differential equations for the evolution of $\hat {\rho }$ over propagation distance $z$ and frequency $\omega$. The evolution equations closely resemble the Liouville equation of the density matrix operator in quantum mechanics, with $z$ or $\omega$ playing the role of time $t$, and the modal-propagation (MP) operator $\hat {Z}$ or group-delay (GD) operator $\hat {\Omega }$ playing the role of Hamiltonian. We further link $\hat {Z}$, $\hat {\Omega }$, and $\hat {\Omega }_{\omega }\equiv \frac {\partial }{\partial \omega }\hat {\Omega }$ to essential modal properties of the fiber, and derive the evolution equations of these operators. Particularly useful for numerical analysis, we establish the concatenation rules that determine $\hat {\Omega }$, $\hat {\Omega }_{\omega }$ of a concatenated fiber-segment assembly from $\hat {\Omega }$, $\hat {\Omega } _{\omega }$ of the individual segments. Combining the density-matrix framework and the stochastic-process theory, we formulate a statistical model for the study of multiple-mode coupling and dispersion in randomly-perturbed fibers, which could be critical to the performance of MDM communications networks. Finally, as a typical example to illustrate the application of our proposed formalism, we present detailed numerical simulations and analysis for the statistical properties of coupling and dispersion in a $4$-dimensional modal space with realistic fiber parameters. We emphasize that our work mainly concerns a general theoretical formalism (beyond a specific model) to treat modal coupling and dispersion in optical communications, and the essential distinctions between our theoretical formalism and the conventional Stokes-vector formalism [23,24] are the fundamental elements (particularly, density matrix versus Stokes vector) characterizing field states and modal properties.
The paper is organized as follows. In Sec. 2, after briefly reviewing the Stokes formalism and discussing its application, we introduce the density matrix operator $\hat {\rho }$ to fully describe general field states in a modal space of arbitrary dimension $J$; in Sec. 3, we derive the differential equations for the evolution of $\hat {\rho }$ and other key elements (e.g., $\hat {\Omega }$ and $\hat {\Omega }_{\omega }$) in our proposed formalism; in Sec. 4, we establish the concatenation rules for $\hat {\Omega }$ and $\hat {\Omega }_{\omega }$; in Sec. 5, we construct the statistical model for multiple-mode coupling and dispersion in randomly-perturbed fibers; in Sec. 6, we numerically simulate and analyze the statistical modal properties of a fiber that supports $4$ orthogonal modes; a summary of our main work is given in Sec. 7.
2. Density matrix operator for description of general field states
A monochromatic electromagnetic field of frequency $\omega$ propagating (along the $z$-direction) in a fiber can be written as $\vec {E}\left ( x,y,z\right ) e^{i\omega t}$, with $\vec {E}$ being the electric field amplitude. Throughout this paper, unless explicitly stated otherwise, we assume all fields to be monochromatic and omit the time-harmonic factor $e^{i\omega t}$. In close analogy to quantum mechanics, the field $\vec {E} \left ( x,y,z\right )$ at a longitudinal position $z$ is described by a ket vector $\left \vert \psi \left ( z\right ) \right \rangle$, which can presumably be expanded by $J$ orthogonal and normalized modes, formally called a “basis” and typically (but not necessarily) chosen to be the “free modes”, i.e., the ideal-fiber eigenmodes $\left \{ \left \vert j\right \rangle _{f},j=1,2,\ldots ,J\right \}$ of effective refractive indices (ERIs) $n_{j}^{\left ( 0\right ) }$. If the fiber is perturbed at position $z$, the local eigenmodes will change to $\left \{ \left \vert j;z\right \rangle _{e},j=1,2,\ldots ,J\right \}$ of ERIs $n_{j}\left ( z\right )$, connected to the free modes via a unitary $J\times J$ matrix $\bar {\Gamma }\left ( z\right )$
$$ \left( \left\vert 1;z\right\rangle _{e},\left\vert 2;z\right\rangle _{e},\ldots,\left\vert J;z\right\rangle _{e}\right) ^{\mathrm{T}}=\bar{\Gamma} \left( z\right) \left( \left\vert 1\right\rangle _{f},\left\vert 2\right\rangle _{f},\ldots,\left\vert J\right\rangle _{f}\right) ^{\mathrm{T}}, \bar{\Gamma}^{\dagger }\left( z\right) \bar{\Gamma}\left( z\right) =\bar{I}, $$
In any selected basis $\left \{ \left \vert 1\right \rangle ,\left \vert 2\right \rangle ,\ldots ,\left \vert J\right \rangle \right \}$, a ket $\left \vert \psi \right \rangle$ can be represented by an explicit column vector, i.e., the Jones vector $\psi ,$
(1)$$\psi =\left( \psi _{1},\psi _{2},\ldots,\psi _{J}\right) ^{\mathrm{T}},\;\psi _{j}=\left\langle j\right. \left\vert \psi \right\rangle ,$$
with $\left \langle j\right . \left \vert \psi \right \rangle$ being the inner product of $\left \vert j\right \rangle$ and $\left \vert \psi \right \rangle$. Obviously, given $\psi$ in the basis, $\left \vert \psi \right \rangle$ is uniquely determined, and vice versa. Therefore, the ket $\left \vert \psi \right \rangle$ and the Jones vector $\psi$ are normally referenced interchangeably, and the terms “modal space” (with emphasis on physical fields) and “Jones space” (with emphasis on mathematical descriptions) both refer to the $J$-dimensional vector space for $\left \vert \psi \right \rangle$ or $\psi$. Propagation of the electromagnetic field is characterized by the evolution of state $\left \vert \psi \left ( z\right ) \right \rangle$ over $z$, and, in absence of loss and dissipation, the evolution goes unitarily and preserves the state modulus $\left \langle \psi \left ( z\right ) \right . \left \vert \psi \left ( z\right ) \right \rangle$. Unless otherwise indicated, ket vectors (and correspondingly the modes) in this paper are always understood to be normalized, i.e., $\left \langle \psi \right . \left \vert \psi \right \rangle$ $=1$.2.1 Review of Stokes vectors
Before proposing and developing our theoretical formalism, we take a brief review of the Jones space and Stokes space in the polarization mode dispersion (PMD) theory, where the electromagnetic fields are expanded by two orthogonally-polarized modes (e.g., the $x$- and $y$-polarization states $\left \vert x\right \rangle$ and $\left \vert y\right \rangle$) and the modal space is of dimension $J=2$. Most of the materials in this subsection (Sec. 2.1) are already standard knowledge in the PMD theory [6,8], and they are presented here mainly for two purposes: (1) to make a concise comparison of the Stokes and Jones formalisms, so as to explore potential advantages of the Jones formalism if provided with a theoretical element (later introduced as the “density matrix operator” $\hat {\rho }$, a main point of this paper) that can describe both pure and mixed field states; (2) to discuss several key Jones-space concepts that will be generalized to multiple-mode fields as we develop our “density-matrix formalism” of modal coupling and dispersion.
For a coherent superposition of states $\left \vert \psi ^{\left ( 1\right ) }\right \rangle ,\left \vert \psi ^{\left ( 2\right ) }\right \rangle ,\ldots$, the electromagnetic field can be properly described by $\left \vert \psi \right \rangle =\alpha _{1}\left \vert \psi ^{\left ( 1\right ) }\right \rangle +\alpha _{2}\left \vert \psi ^{\left ( 2\right ) }\right \rangle +\cdots$, where $\alpha _{1},\alpha _{2},..$ are complex coefficients. However, for an incoherent mixture of states $\left \vert \psi ^{\left ( 1\right ) }\right \rangle ,\left \vert \psi ^{\left ( 2\right ) }\right \rangle \cdots$ respectively with noise variance $p_{1},p_{2},\ldots$ (analogous to the “probability” for mixed states in quantum mechanics, although $p_{1}+p_{2}+\cdots$ is not necessarily unity here), it is incorrect to describe the field by superposing $\left \vert \psi ^{\left ( 1\right ) }\right \rangle ,\left \vert \psi ^{\left ( 2\right ) }\right \rangle \cdots$ as above, since the ket vector $\left \vert \psi \right \rangle$ can only characterize pure states. In the PMD theory, the general field state (either pure or mixed) is represented by the Stokes vector $\vec {S}$,
(2)$$\vec{S}=\sum_{q=1,2,..}p_{q}\vec{S}^{\left( q\right) },$$
where $\vec {S}^{\left ( q\right ) }$ $\equiv \left ( S_{1}^{\left ( q\right ) },S_{2}^{\left ( q\right ) },S_{3}^{\left ( q\right ) }\right ) ^{\mathrm {T}}$ are the Stokes vector corresponding to the individual pure states $\left \vert \psi ^{\left ( q\right ) }\right \rangle ,$ (3)$$S_{1}^{\left( q\right) }=\psi _{x}^{\left( q\right) }\psi _{x}^{\left( q\right) \ast }-\psi _{y}^{\left( q\right) }\psi _{y}^{\left( q\right) \ast },\;S_{2}^{\left( q\right) }=2\textrm{Re}\psi _{x}^{\left( q\right) }\psi _{y}^{\left( q\right) \ast },\;S_{3}^{\left( q\right) }=-2\textrm{Im}\psi _{x}^{\left( q\right) }\psi _{y}^{\left( q\right) \ast },$$
with $\psi _{x}^{\left ( q\right ) }=\left \langle x\right . \left \vert \psi ^{\left ( q\right ) }\right \rangle ,\;\psi _{y}^{\left ( q\right ) }=\left \langle y\right . \left \vert \psi ^{\left ( q\right ) }\right \rangle$. In the current paper, we always assume that the total variance $p_{1}+p_{2}+\cdots =1$, and generalizations beyond this restriction are straightforward. Obviously, $\vec {S}$ in Eq. (2) represents a pure state if there exists only one term in the summation. For clarity, we refer to the $3$-dimensional real-vector space for $\vec {S}$ as the “Stokes space”, which corresponds to a $2$-dimensional complex-vector modal space or Jones space for $\left \vert \psi \right \rangle$. Evidently, the essential advantage of the Stokes vector over the Jones vector is for the description of mixed states. A less essential, but nonetheless a huge technical advantage of the Stokes vector is that it can be illustratively depicted as a geometrical vector in the Poincaré Sphere.In a lossless fiber, a pure state $\left \vert \psi \left ( z,\omega \right ) \right \rangle$ evolves over $z$ according to [6,8]
(4)$$i\frac{\partial }{\partial z}\left\vert \psi \left( z,\omega \right) \right\rangle =\hat{Z}\left( z,\omega \right) \left\vert \psi \left( z,\omega \right) \right\rangle ,$$
where we have explicitly allowed the frequency ($\omega$) dependence in the field state. $\hat {Z}\left ( z,\omega \right )$ is an Hermitian operator ($\hat {Z}^{\dagger }=\hat {Z})$ in the Jones space, and its eigenvectors and eigenvalues represent the local eigenmodes (LEs) $\left \vert l_{\pm };z,\omega \right \rangle$ and the corresponding propagation constants $n_{\pm }\left ( z,\omega \right ) \frac {\omega }{c}$, respectively [$n_{\pm }$ are the ERIs, with $n_{+}\geq n_{-}$, and $c$ is the speed of light in vacuum]. The Stokes-space counterpart of Eq. (4), valid for the general Stokes vector in Eq. (2), reads (5)$$\frac{\partial }{\partial z}\vec{S}\left( z,\omega \right) =\vec{\beta} \left( z,\omega \right) \times \vec{S}\left( z,\omega \right) ,$$
where (6)$$\;\vec{\beta}=\left( \beta _{1},\beta _{2},\beta _{3}\right) ^{\mathrm{T} },\;\beta _{1}=\bar{Z}_{xx}-\bar{Z}_{yy},\;\beta _{2}=2\textrm{Re}\left( \bar{Z }_{xy}\right) ,\;\beta _{3}=-2\textrm{Im}\left( \bar{Z}_{xy}\right) ,$$
is the “birefrigence vector” in the Stokes space, with $\bar {Z}_{jj^{\prime }}=\left \langle j\right \vert \hat {Z}\left \vert j^{\prime }\right \rangle ,$ ($j,j^{\prime }=x$ or $y$). The modulus $\left \vert \vec {\beta }\right \vert$ and unit vector $\frac {\vec { \beta }}{\left \vert \vec {\beta }\right \vert }$ of the birefrigence vector give the propagation-constant difference $\left ( n_{+}-n_{-}\right ) \frac {\omega }{c}$ and the LE of larger propagation constant $n_{+}\frac {\omega }{c}$ (or ERI $n_{+}$), respectively.The solution of Eq. (4) can be formally written as
(7)$$\left\vert \psi \left( z,\omega \right) \right\rangle =\hat{U}\left( z,z_{ \mathrm{I}},\omega \right) \left\vert \psi \left( z_{\mathrm{I}},\omega \right) \right\rangle ,$$
where $z_{\mathrm {I}}$ is the “initial” position, and the evolution operator $\hat {U}\left ( z,z_{\mathrm {I}},\omega \right )$ obeys (8)$$i\frac{\partial }{\partial z}\hat{U}\left( z,z_{\mathrm{I}},\omega \right) = \hat{Z}\left( z,\omega \right) \hat{U}\left( z,z_{\mathrm{I}},\omega \right) ,\;\hat{U}\left( z_{\mathrm{I}},z_{\mathrm{I}},\omega \right) =\hat{I},$$
with $\hat {I}$ being the unit operator. From Eq. (8) and the Hermiticity of $\hat {Z}\left ( z,\omega \right )$, it follows that $\hat {U} \left ( z,z_{\mathrm {I}},\omega \right )$ is unitary (9)$$\hat{U}^{\dagger }\left( z,z_{\mathrm{I}},\omega \right) \hat{U}\left( z,z_{ \mathrm{I}},\omega \right) =\hat{U}\left( z,z_{\mathrm{I}},\omega \right) \hat{U}^{\dagger }\left( z,z_{\mathrm{I}},\omega \right) =\hat{I}.$$
Eq. (4), or alternatively, Eqs. (7) and (8), govern the evolution of field state $\left \vert \psi \right \rangle$ over propagation position $z$. For the variation of $\left \vert \psi \right \rangle$ over frequency $\omega$, assuming that the state at initial position $z_{\mathrm {I}}$ does not depend on $\omega$, i.e., $\left \vert \psi \left ( z_{\mathrm {I}},\omega \right ) \right \rangle =\left \vert \psi \left ( z_{\mathrm {I}}\right ) \right \rangle$, we take the $\omega$-derivative of Eq. (7) and apply the unitarity of $\hat {U}$ to obtain (10)$$i\frac{\partial }{\partial \omega }\left\vert \psi \left( z,\omega \right) \right\rangle =\hat{\Omega}\left( z,z_{\mathrm{I}},\omega \right) \left\vert \psi \left( z,\omega \right) \right\rangle ,$$
where $\hat {\Omega }$ is the “group-delay operator” (11)$$\hat{\Omega}\left( z,z_{\mathrm{I}},\omega \right) =i\hat{U}_{\omega }\left( z,z_{\mathrm{I}},\omega \right) \hat{U}^{\dagger }\left( z,z_{\mathrm{I} },\omega \right) ,$$
with $\hat {U}_{\omega }\left ( z,z_{\mathrm {I}},\omega \right ) \equiv \frac { \partial }{\partial \omega }\hat {U}\left ( z,z_{\mathrm {I}},\omega \right )$. $\hat {\Omega }\left ( z,z_{\mathrm {I}},\omega \right )$ is an Hermitian Jones-space operator, whose eigenvectors $\left \vert \mathrm {p}_{\pm };z,z_{ \mathrm {I}},\omega \right \rangle$ and eigenvalues $\tau _{\mathrm {p}_{\pm }}\left ( z,z_{\mathrm {I}},\omega \right )$ [$\tau _{\mathrm {p}_{+}}\geq \tau _{\mathrm {p}_{-}}$] represent the principal states of the fiber segement from $z_{\mathrm {I}}$ to $z$ and the corresponding group delays, respectively. Equation (10) has the Stokes-space counterpart that is valid for the general Stokes vector in Eq. (2), (12)$$\frac{\partial }{\partial \omega }\vec{S}\left( z,\omega \right) =\vec{\tau} \left( z,z_{\mathrm{I}},\omega \right) \times \vec{S}\left( z,\omega \right) ,$$
where $\vec {\tau }$ is the “PMD vector” (13)$$\vec{\tau}=\left( \tau _{1},\tau _{2},\tau _{3}\right) ^{\mathrm{T}},\;\tau _{1}=\bar{\Omega}_{xx}-\bar{\Omega}_{yy},\;\tau _{2}=2\textrm{Re}\left( \bar{ \Omega}_{xy}\right) ,\;\tau _{3}=-2\textrm{Im}\left( \bar{\Omega}_{xy}\right) ,$$
with $\bar {\Omega }_{jj^{\prime }}=\left \langle j\right \vert \hat {\Omega } \left \vert j^{\prime }\right \rangle$ ($j,j^{\prime }=x$ or $y$). The modulus $\left \vert \vec {\tau }\right \vert$ and unit vector $\frac {\vec {\tau } }{\left \vert \vec {\tau }\right \vert }$ of the PMD vector give the group-delay difference $\left ( \tau _{\mathrm {p}_{+}}-\tau _{\mathrm {p}_{-}}\right )$ and the principal state of larger group delay $\tau _{\mathrm {p}_{+}}$, respectively.It is instructive to compare $\vec {\beta }$, $\vec {\tau }$ (Stokes-space vectors) and $\hat {Z}$, $\hat {\Omega }$ (Jones-space opertors) in the analysis of modal coupling and dispersion. In the PMD theory, i.e., for a modal space of dimension $J=2$, both $\vec {\beta }$ and $\hat {Z}$ yield the local eigenmodes and the difference of propagation constants, and similarly, both $\vec {\tau }$ and $\hat {\Omega }$ yield the principal states and the difference of group delays. Although $\hat {Z}$ and $\hat {\Omega }$ actually give the propagation constant for each local eigenmode and the group delay for each principal state, in most situations, only the differences really matter. Thus, for $J=2$, $\vec {\beta },\vec {\tau }$ and $\hat {Z}$, $\hat {\Omega }$ provide essentially the same modal information of the fiber. However, for modal spaces of dimension $J>2$, the application of $\vec {\beta },\vec {\tau }$, e.g., accessing the modal eigenvectors (local eigenmodes and principal states) and eigenvalues (propagation constants and group delays), needs the assistance of $J^{2}-1$ traceless and mutually orthogonal Gell-Mann matrices [19,20], while the extension of $\hat {Z}$, $\hat {\Omega }$ to any higher-dimensional ($J>2$) modal space is straightforward, and the modal eigenvectors and eigenvalues can be directly acquired from $\hat {Z}$, $\hat {\Omega }$ (see Sec. 3).
Still, the Stokes space is advantageous over the Jones space in one crucial aspect (even for multiple-mode fields), i.e., the Stokes vector $\vec {S}$ can describe mixed field states while the Jones vector $\left \vert \psi \right \rangle$ cannot, and, for $J=2,$ $\vec {S}$ can be treated as a geometric vector in the Poincaré Sphere. On the other hand, in higher-dimensional ($J>2$) modal spaces, the application of $\vec {S }$ also require $J^{2}-1$ auxiliary Gell-Mann matrices, and the geometric advantage of the Poincaré Sphere is largely lost since the dimension ($J^{2}-1$) of $\vec {S}$ is now larger than $3$. Therefore, it should be desirable to device some Jones-space tool that is simple in definition and straightforward in application to represent general (both pure and mixed) field states for modal spaces of arbitrary dimension. In Sec. 2.2, we will introduce such a theoretical tool $\hat {\rho }$ in the Jones space, named “density matrix operator”, as the central concept of our proposed formalism for modal coupling and dispersion.
2.2 Density matrix operator in Jones space
From now on, all modal spaces under consideration will be of arbitrary dimension $J$, unless explicitly specified otherwise. As already mentioned, in the Jones space one uses the ket vector $\left \vert \psi \right \rangle$ to characterize the state of a pure field. Alternatively, again in close analogy to quantum mechanics, we can introduce the “density matrix operator” $\hat {\rho } =\left \vert \psi \right \rangle \left \langle \psi \right \vert$ to describe the field state, since any physical information derivable from $\left \vert \psi \right \rangle$ can also be achieved from $\hat {\rho }$, and vice versa. Given a basis $\left \{ \left \vert 1\right \rangle ,\left \vert 2\right \rangle ,\ldots ,\left \vert J\right \rangle \right \}$, $\left \vert \psi \right \rangle$ is represented by the Jones (column) vector $\psi$ in Eq. (1 ). Accordingly, the density matrix operator $\hat {\rho }$ is represented by an $J\times J$ matrix $\bar {\rho }$ of entries $\bar {\rho }_{jj^{\prime }}\equiv \left \langle j\right \vert \hat {\rho }\left \vert j^{\prime }\right \rangle =\left \langle j\right . \left \vert \psi \right \rangle \left \langle \psi \right . \left \vert j^{\prime }\right \rangle =\psi _{j}\psi _{j^{\prime }}^{\ast }$, $ j,j^{\prime }=1,2,\ldots ,J$. If one further writes $\left \langle \psi \right \vert$ in the basis as a row vector $\left ( \left \langle \psi \right . \left \vert 1\right \rangle ,\left \langle \psi \right . \left \vert 2\right \rangle ,\ldots ,\left \langle \psi \right . \left \vert J\right \rangle \right ) =\psi ^{\dagger }$, then $\bar {\rho }$ can be written as $\bar {\rho }=\psi \psi ^{\dagger }$. Naturally, $\bar {\rho }$ is called “density matrix”, and, in absence of ambiguity, we will frequently refer to $\hat {\rho }$ and $\bar {\rho }$, as well as their names, interchangeably.
For an incoherent mixture of $\left \vert \psi ^{\left ( 1\right ) }\right \rangle ,\left \vert \psi ^{\left ( 2\right ) }\right \rangle ,\ldots$ (with noise variance $p_{1},p_{2},\ldots$, respectively), the field state is fully described by the density matrix operator in the general form,
(14)$$\hat{\rho}=\sum_{q=1,2,\ldots}p_{q}\hat{\rho}^{\left( q\right) },\;\hat{\rho} ^{\left( q\right) }=\left\vert \psi ^{\left( q\right) }\right\rangle \left\langle \psi ^{\left( q\right) }\right\vert ,$$
which would reduce to the pure-state form if only one term exists in the summation. In any chosen basis, the corresponding density matrix reads $\bar {\rho }=\sum _{q=1,2,\ldots }p_{q}\psi ^{\left ( q\right ) }\psi ^{\left ( q\right ) \dagger }$, with $\psi ^{\left ( q\right ) }$ being the Jones vector for $\left \vert \psi ^{\left ( q\right ) }\right \rangle$. Like in quantum mechanics, all observable information of the field state can be extracted from $\hat {\rho }$, which is an Hermitian operator of $J$ real and non-negative eigenvalues. From $\left \langle \psi ^{\left ( q\right ) }\right \vert \left . \psi ^{\left ( q\right ) }\right \rangle =1$ and $\sum _{q=1,2,\ldots }p_{q}=1$, it follows that Tr$\hat {\rho }=1$, where “Tr” stands for the trace of an operator or matrix. Also, one has $\frac{1}{N} \leq \operatorname{Tr} \hat{\rho}^{2} \leq 1$ (Tr$\hat {\rho }^{2}=1$ iff $\hat {\rho }$ is a pure state).In principle, a pure electromagnetic field $\left \vert \psi ^{\left ( q\right ) }\right \rangle$ can be decomposed (via proper apparatus) into any complete set of orthogonal modes $\left \{ \left \vert \phi ^{\left ( 1\right ) }\right \rangle ,\left \vert \phi ^{\left ( 2\right ) }\right \rangle ,\ldots ,\left \vert \phi ^{\left ( J\right ) }\right \rangle \right \}$, e.g., demultiplexing a signal field into its consituent modes. The portion of energy in mode $\left \vert \phi ^{\left ( j\right ) }\right \rangle$, i.e., ratio $P_{j}^{\left ( q\right ) }$ of the energy $E_{j}^{\left ( q\right ) }$ in $\left \vert \phi ^{\left ( j\right ) }\right \rangle$ to the total energy $E_{T}^{\left ( q\right ) }$ in $\left \vert \psi ^{\left ( q\right ) }\right \rangle$, can be made analogous to the measurement probability ($\sim P_{j}^{\left ( q\right ) }$) of finding a quantum system ($\sim \left \vert \psi ^{\left ( q\right ) }\right \rangle$) in an eigenstate ($\sim \left \vert \phi ^{\left ( j\right ) }\right \rangle$) of some observable $\hat {O}$ [25]$,$ and is thus given by
$$ P_{j}^{\left( q\right) }\equiv \frac{E_{j}^{\left( q\right) }}{E_{T}^{\left( q\right) }}=\left\vert \left\langle \phi ^{\left( j\right) }\right\vert \left. \psi ^{\left( q\right) }\right\rangle \right\vert ^{2}=\left\langle \phi ^{\left( j\right) }\right\vert \left. \psi ^{\left( q\right) }\right\rangle \left\langle \psi ^{\left( q\right) }\right\vert \left. \phi ^{\left( j\right) }\right\rangle . $$
$$ C_{jj^{\prime }}^{\left( q\right) }=\left( \left\langle \phi ^{\left( j\right) }\right\vert \left. \psi ^{\left( q\right) }\right\rangle \right) \left( \left\langle \phi ^{\left( j^{\prime }\right) }\right\vert \left. \psi ^{\left( q\right) }\right\rangle \right) ^{\ast }=\left\langle \phi ^{\left( j\right) }\right\vert \left. \psi ^{\left( q\right) }\right\rangle \left\langle \psi ^{\left( q\right) }\right\vert \left. \phi ^{\left( j^{\prime }\right) }\right\rangle . $$
(15)$$\begin{aligned} P_{j} &=\sum_{q}p_{q}P_{j}^{\left( q\right) }=\left\langle \phi ^{\left( j\right) }\right\vert \left( \sum_{q}p_{q}\left\vert \psi ^{\left( q\right) }\right\rangle \left\langle \psi ^{\left( q\right) }\right\vert \right) \left\vert \phi ^{\left( j\right) }\right\rangle =\left\langle \phi ^{\left( j\right) }\right\vert \hat{\rho}\left\vert \phi ^{\left( j\right) }\right\rangle ,\\ C_{jj^{\prime }} &=\sum_{q}p_{q}C_{jj^{\prime }}^{\left( q\right) }=\left\langle \phi ^{\left( j\right) }\right\vert \left( \sum_{q}p_{q}\left\vert \psi ^{\left( q\right) }\right\rangle \left\langle \psi ^{\left( q\right) }\right\vert \right) \left\vert \phi ^{\left( j^{\prime }\right) }\right\rangle =\left\langle \phi ^{\left( j\right) }\right\vert \hat{\rho}\left\vert \phi ^{\left( j^{\prime }\right) }\right\rangle . \end{aligned}$$
More generally, for the information of an observable quantity $A$ (e.g., the orbital angular momentum) with orthogonal eigenmodes $\left \{ \left \vert \phi _{A}^{\left ( j=1,2,\ldots ,J\right ) }\right \rangle \right \}$ and corresponding eigenvalues $a_{j=1,2,\ldots ,J}$, Eq. (15) gives the portion of energy in $\left \vert \phi _{A}^{\left ( j\right ) }\right \rangle$ as $\left \langle \phi _{A}^{\left ( j\right ) }\right \vert \hat {\rho } \left \vert \phi _{A}^{\left ( j\right ) }\right \rangle$, and thus the mean value of $A$ is (16)$$\left\langle A\right\rangle =\sum_{j=1}^{J}P_{j}a_{j}=\sum_{j=1}^{J}\left\langle \phi _{A}^{\left( j\right) }\right\vert \hat{\rho}\left\vert \phi _{A}^{\left( j\right) }\right\rangle a_{j}=\textrm{Tr}\left( \hat{\rho}\sum_{j=1}^{J}a_{j}\left\vert \phi _{A}^{\left( j\right) }\right\rangle \left\langle \phi _{A}^{\left( j\right) }\right\vert \right) \equiv \textrm{Tr}\left( \hat{\rho}\hat{A} \right) .$$
Here, by defining operator $\hat {A}=\sum _{j=1}^{J}a_{j}\left \vert \phi _{A}^{\left ( j\right ) }\right \rangle \left \langle \phi _{A}^{\left ( j\right ) }\right \vert$ to represent the observable quantity $A$, the formula for $\left \langle A\right \rangle$ has been cast into the same form as in quantum mechanics. In principle, all observable information of the field state can be achieved via Eqs. (15) and (16). The simple rules (15)–(16) to acquire physical information in the field state illustrate a main strength of the density-matrix formalism.We conclude the section with two important remarks. First, one should realize that the Stokes-vector formalism and the density-matrix formalism are eventually equivalent in the mathematical sense, and thus it is in principle possible to access via the Stokes vector the same physical information extracted from the density matrix. Indeed, some physical interpretation of the Stokes parameters was achieved by representing the multi-modal signal in high-order Poincaré spheres [26,27], and the Stokes formalism could be used for modal-dispersion/modal-dependent-loss measurement with direct-detection receivers [28]. Nonethelss, the density matrix appears to be capable of more directly treating or yielding some crucial information, e.g., the energy distribution $P_{j}$ and the cohernce $C_{jj^{\prime }}$ as in Eq. (15). Another important remark states that our proposed density matrix formalism borrows methodology from quantum mechanics, and similar transfer of methodology from quantum mechanics to optics can be found in literature. For example, the “coherency matrix” employed in Ref. [29] to help characterize the noise fields in multiple-mode-fiber optical communications channels is also mathematically equivalent to the density matrix in quantum mechanics, and it could be of great interest to unify the density-matrix or coherency-matrix approaches in different realms of optics and broaden the scope of this useful theoretical tool.
3. Evolution equations of density matrix operator and other key elements
In this section, we construct the fundamental equations that govern the evolution of density matrix operator $\hat {\rho }$ and other key elements that characterize essential modal properties. We emphasize that all the analysis and conclusions are valid for modal spaces of any dimension $J\geq 2$.
In a lossless fiber, a monochromatic pure field of frequency $\omega$ evolves over the propagation position $z$ according to the direct generalization of Eq. (4) from $J=2$ to arbitrary $J$. The Hermitian “modal propagation (MP)” operator $\hat {Z}\left ( z,\omega \right )$ can be written in terms of its eigenvectors $\left \vert l_{j};z,\omega \right \rangle$ and eigenvalues $\frac {\omega }{c}n_{j}\left ( z,\omega \right )$ ($j=1,2,\ldots ,J$)
$$ \hat{Z}\left( z,\omega \right) =\frac{\omega }{c}\sum_{j=1}^{J}n_{j}\left( z,\omega \right) \left\vert l_{j};z,\omega \right\rangle \left\langle l_{j};z,\omega \right\vert \equiv \frac{\omega }{c}\hat{N}\left( z,\omega \right) , $$
(17)$$i\frac{\partial }{\partial z}\hat{\rho}\left( z,\omega \right) =\hat{Z} \left( z,\omega \right) \hat{\rho}\left( z,\omega \right) -\hat{\rho}\left( z,\omega \right) \hat{Z}\left( z,\omega \right) \equiv \left[ \hat{Z}\left( z,\omega \right) ,\hat{\rho}\left( z,\omega \right) \right] ,$$
where we have introduced the commutator of two operators, $\left [ \hat {A}, \hat {B}\right ] \equiv \hat {A}\hat {B}-\hat {B}\hat {A}$. Furthermore, for a mixed state, e.g., the initial field at $z_{\mathrm {I}}$ being a mixture of pure states $\left \vert \psi ^{\left ( 1\right ) }\left ( z_{\mathrm {I}}\right ) \right \rangle ,\left \vert \psi ^{\left ( 2\right ) }\left ( z_{\mathrm {I} }\right ) \right \rangle ,\ldots$ with noise variance $p_{1},p_{2},\ldots$, since the evolution of each individual pure state $\hat {\rho }^{\left ( q=1,2,\ldots \right ) }\left ( z,\omega \right ) =\left \vert \psi ^{\left ( q\right ) }\left ( z,\omega \right ) \right \rangle \left \langle \psi ^{\left ( q\right ) }\left ( z,\omega \right ) \right \vert$ is governed by Eq. (17), the overall density matrix operator $\hat {\rho }\left ( z,\omega \right ) =\sum _{q=1,2,\ldots }p_{q}\left \vert \psi ^{\left ( q\right ) }\left ( z,\omega \right ) \right \rangle \left \langle \psi ^{\left ( q\right ) }\left ( z,\omega \right ) \right \vert$ also obeys the same equation. Therefore, the fundamental equation governing the propagation (in the $z$ -direction) of a general field state $\hat {\rho }$ is Eq. (17), which closely resembles the Liouville equation for the density matrix in quantum mechanics, with $\hat {Z}$ and $z$ playing the role of Hamiltonian and time, respectively. In a basis $\left \{ \left \vert 1\right \rangle ,\left \vert 2\right \rangle ,\ldots .,\left \vert J\right \rangle \right \}$, one has the matrix representation of Eq. (17), $i \frac {\partial }{\partial z}\bar {\rho }\left ( z,\omega \right ) =\left [ \bar {Z} \left ( z,\omega \right ) ,\bar {\rho }\left ( z,\omega \right ) \right ]$. For clarity, we note that, given the basis, each Jones-space equation of operators $\hat {A}$ ($\hat {A}=\hat {\rho },\hat {Z},\cdots$) and/or ket vectors $\left \vert \psi \right \rangle$ corresponds to a representation equation in exactly the same form, with $\hat {A}$ and $\left \vert \psi \right \rangle$ replaced by matrices $\bar {A}$ and Jones (column) vectors $\psi$, i.e., $\hat {A}\rightarrow \bar {A}$: $\bar {A}_{jj^{\prime }}=\left \langle j\right \vert \hat {A}\left \vert j^{\prime }\right \rangle$, $ \left \vert \psi \right \rangle \rightarrow \psi$: $\psi _{j}=\left \langle j\right \vert \left . \psi \right \rangle$, $ j,j^{\prime }=1,2,\ldots ,J$. So, in absence of ambiguity, we will regard all Jones-space operators and ket vectors as synonymous with the corresponding matrices and Jones vectors. Accordingly, all equations of Jones-space operators and/or ket vectors will be regarded as synonymous with the corresponding equations of matrices and/or Jones vectors.The direct generalization of Eqs. (7) – (9) to modal spaces of arbitrary dimension $J$ formally solves for $\left \vert \psi \left ( z,\omega \right ) \right \rangle$ via the evolution operator $\hat {U}.$ It follows from Eq. (7) that the formal solution of Eq. (17) can be expressed in terms of $\hat {U}$,
(18)$$\hat{\rho}\left( z,\omega \right) =\hat{U}\left( z,z^{\prime },\omega \right) \hat{\rho}\left( z^{\prime },\omega \right) \hat{U}^{\dagger }\left( z,z^{\prime },\omega \right) ,$$
where a more flexible notation $z^{\prime }$ (instead of $z_{\mathrm {I}}$) has been adopted to denote the initial position. Also, right-multiplying both sides of the first equation in (8) by $\hat {U}^{\dagger }\left ( z,z^{\prime },\omega \right )$ and applying Eq. (9) yields the relation of $\hat {Z}$ and $\hat {U}$, (19)$$\hat{Z}\left( z,\omega \right) =i\hat{U}_{z}\left( z,z^{\prime },\omega \right) \hat{U}^{\dagger }\left( z,z^{\prime },\omega \right) ,$$
with $\hat {U}_{z}\left ( z,z^{\prime },\omega \right ) \equiv \frac {\partial }{ \partial z}\hat {U}\left ( z,z^{\prime },\omega \right )$.Equation (17), or equivalently Eqs. (8) and (18 ), govern the propagation of monochromatic fields along a fiber, i.e., the variation of $\hat {\rho }$ over propagation position $z$. To analyze modal dispersion, it is necessary to determine the variation of $\hat {\rho }$ over frequency $\omega .$ To this end, we first note that Eqs. (10) and (11) can also be directly extended to arbitrary $J$-dimensional modal spaces. Again, we emphasize that Eq. (10) is valid only if the “initial” state $\left \vert \psi \left ( z^{\prime }\right ) \right \rangle$ does NOT depend on frequency $\omega$. $\hat {\Omega }\left ( z,z^{\prime },\omega \right )$ is an Hermitian operator that can be written in terms of its $J$ eigenvectors $\left \vert \mathrm {p} _{j=1,2,\ldots ,J};z,z^{\prime },\omega \right \rangle$ and eigenvalues $\tau _{ \mathrm {p}_{j=1,2,\ldots ,J}}\left ( z,z^{\prime },\omega \right )$,
$$ \hat{\Omega}\left( z,z^{\prime },\omega \right) =\sum_{j=1}^{J}\tau _{ \mathrm{p}_{j}}\left( z,\omega \right) \left\vert \mathrm{p}_{j};z,z^{\prime },\omega \right\rangle \left\langle \mathrm{p}_{j};z,z^{\prime },\omega \right\vert . $$
(20)$$i\frac{\partial }{\partial \omega }\hat{\rho}\left( z,\omega \right) =\left[ \hat{\Omega}\left( z,z^{\prime },\omega \right) ,\hat{\rho}\left( z,\omega \right) \right] ,$$
where, for consistency, the field state at the initial position $z^{\prime }$ has to be $\omega$-independent, $\hat {\rho }\left ( z^{\prime },\omega \right ) =\hat {\rho }\left ( z^{\prime }\right ) .$Unlike $\hat {Z}\left ( z,\omega \right )$ that is determined by local properties of the fiber at $z$, $\hat {\Omega }\left ( z,z^{\prime },\omega \right )$ depends on global properties of the fiber segment from $z^{\prime }$ to $z$. The differential equation governing the evolution of $\hat {\Omega }$ over $z$ can be achieved by taking the $z$-derivative of Eq. (11),
(21)$$i\frac{\partial }{\partial z}\hat{\Omega}\left( z,z^{\prime },\omega \right) =i\hat{Z}_{\omega }\left( z,\omega \right) +\left[ \hat{Z}\left( z,\omega \right) ,\hat{\Omega}\left( z,z^{\prime },\omega \right) \right] ,$$
with $\hat {Z}_{\omega }\equiv \frac {\partial }{\partial \omega }\hat {Z}$. Here, we have applied Eqs. (9) and (19), as well as (22)$$\begin{aligned} \hat{U}_{z}\left( z,z^{\prime },\omega \right) \hat{U}^{\dagger }\left( z,z^{\prime },\omega \right) &=-\hat{U}\left( z,z^{\prime },\omega \right) \hat{U}_{z}^{\dagger }\left( z,z^{\prime },\omega \right) , \end{aligned}$$
(23)$$\begin{aligned} \hat{U}_{\omega }\left( z,z^{\prime },\omega \right) \hat{U}^{\dagger }\left( z,z^{\prime },\omega \right) &=-\hat{U}\left( z,z^{\prime },\omega \right) \hat{U}_{\omega }^{\dagger }\left( z,z^{\prime },\omega \right) \end{aligned}$$
that follow from Eq. (9). Noting $\hat {U}\left ( z^{\prime },z^{\prime },\omega \right ) =\hat {I}$, one immediately gets the initial condition for $\hat {\Omega }\left ( z,z^{\prime },\omega \right )$ from Eq. (11), (24)$$\hat{\Omega}\left( z^{\prime },z^{\prime },\omega \right) =0.$$
Equations (17) and (21), or equivalently, Eqs. (8), (18) and (11), subject to proper initial conditions, form the foundation of our theoretical formalism for modal coupling and dispersion in multiple-mode fibers. Once the MP operator $\hat {Z}\left ( z,\omega \right )$ are known (e.g., calculated from the optical structure of the fiber), one can solve Eqs. (17) and (21) for $\hat {\rho }\left ( z,\omega \right )$ and $\hat {\Omega } \left ( z,z^{\prime },\omega \right )$ as functions of $z$ at any $\omega$. Alternatively, we can solve Eq. (8) for the evolution operator $\hat {U}\left ( z,z^{\prime },\omega \right )$ as function of $z,\omega$ and further obtain $\hat {\rho }\left ( z,\omega \right )$ and $\hat {\Omega }\left ( z,z^{\prime },\omega \right )$ via Eqs. (18) and (11 ). Subsequently, observable information on modal coupling and dispersion can be extracted from $\hat {\rho }\left ( z,\omega \right )$ and $\hat {\Omega } \left ( z,z^{\prime },\omega \right )$.
The group-delay (GD) operator $\hat {\Omega }\left ( z,z^{\prime },\omega \right )$ yields direct information on first-order modal dispersion (group delay) of the fiber. Second-order dispersion is attained from the first-order $\omega$-derivative of the GD operator, i.e., $\hat { \Omega }_{\omega }\equiv \frac {\partial }{\partial \omega }\hat {\Omega }$, and third-order dispersion from the second-order $\omega$-derivative, $\hat { \Omega }_{\omega \omega }\equiv \frac {\partial ^{2}}{\partial \omega ^{2}} \hat {\Omega }$, and so on. In principle, $\hat {\Omega }_{\omega }$, $\hat { \Omega }_{\omega \omega },\ldots$ at any particular frequency $\omega$ can be achieved by solving for $\hat {\Omega }$ in the neighborhood of $\omega$ and calculating the frequency-derivatives via definition. A practically more feasible way is to construct the evolution equations for $\hat { \Omega }_{\omega }$, $\hat {\Omega }_{\omega \omega }$, $\cdots$ by differentiating Eq. (21) to the desired order and solve these equations together with Eqs. (17) and (21). Here, since in most situations modal dispersion beyond second-order is of essentially no practical significance, we only write down the equation for $\hat {\Omega }_{\omega }$ [following immediately from the $\omega$-derivative of Eq. (21)],
(25)$$\begin{aligned} &i\frac{\partial }{\partial z}\hat{\Omega}_{\omega }\left( z,z^{\prime },\omega \right)\\ &=i\hat{Z}_{\omega \omega }\left( z,\omega \right) +\left[ \hat{Z}_{\omega }\left( z,\omega \right) ,\hat{\Omega}\left( z,z^{\prime },\omega \right) \right] +\left[ \hat{Z}\left( z,\omega \right) ,\hat{\Omega}_{\omega }\left( z,z^{\prime },\omega \right) \right] , \end{aligned}$$
with $\hat {Z}_{\omega \omega }\equiv \frac {\partial ^{2}}{\partial \omega ^{2}}\hat {Z}$. Equation (25) is subject to the initial condition (26)$$\hat{\Omega}_{\omega }\left( z^{\prime },z^{\prime },\omega \right) =0,$$
a direct consequence of Eq. (24).Once $\hat {\Omega }\left ( z,z^{\prime },\omega \right )$ and $\hat {\Omega } _{\omega }\left ( z,z^{\prime },\omega \right )$ at a frequency $\omega$ are available, the second-order modal dispersion information (at $\omega$) can be deduced as follows. For a sufficiently small increment $\delta \omega$ of frequency, one has
(27)$$\hat{\Omega}\left( \omega +\delta \omega \right) =\hat{\Omega}\left( \omega \right) +\hat{\Omega}_{\omega }\left( \omega \right) \delta \omega ,$$
where $z$ and $z^{\prime }$ have been formally omitted to simplify the notations. $\hat {\Omega }\left ( \omega \right )$ and $\hat {\Omega }\left ( \omega +\delta \omega \right )$ are then diagonalized to yield the eigenvalues (group delays) $\tau _{\mathrm {p}_{1}}\left ( \omega \right ) \geq \tau _{\mathrm {p}_{2}}\left ( \omega \right ) \geq \cdots .\geq \tau _{\mathrm {p} _{J}}\left ( \omega \right )$, $\tau _{\mathrm {p}_{1}}\left ( \omega +\delta \omega \right ) \geq \tau _{\mathrm {p}_{2}}\left ( \omega +\delta \omega \right ) \geq \cdots .\geq \tau _{\mathrm {p}_{J}}\left ( \omega +\delta \omega \right )$, and the corresponding eigenvectors (principle states) $\left \vert \mathrm {p}_{j=1,2,\ldots ,J};\omega \right \rangle$, $\left \vert \mathrm {p} _{j=1,2,\ldots ,J};\omega +\delta \omega \right \rangle$. The second-order modal dispersion includes two aspects: (1) mode-dependent chromatic dispersion (MDCD), i.e., the change of group delays with $\omega$, $$ \tau _{\mathrm{p}_{j},\omega }\left( \omega \right) \equiv \frac{\partial \tau _{\mathrm{p}_{j}}\left( \omega \right) }{\partial \omega }=\frac{\tau _{ \mathrm{p}_{j}}\left( \omega +\delta \omega \right) -\tau _{\mathrm{p} _{j}}\left( \omega \right) }{\delta \omega }, $$
$$ \hat{\rho}_{\mathrm{p}_{j},\omega }\left( \omega \right) \equiv \frac{ \partial }{\partial \omega }\left( \left\vert \mathrm{p}_{j};\omega \right\rangle \left\langle \mathrm{p}_{j};\omega \right\vert \right) =\frac{ \left\vert \mathrm{p}_{j},\omega +\delta \omega \right\rangle \left\langle \mathrm{p}_{j},\omega +\delta \omega \right\vert -\left\vert \mathrm{p} _{j},\omega \right\rangle \left\langle \mathrm{p}_{j},\omega \right\vert }{ \delta \omega }. $$
(28)$$\left\vert \vec{\tau}_{\omega \perp }\left( \omega \right) \right\vert =\left( \tau _{\mathrm{p}_{+}}\left( \omega \right) -\tau _{\mathrm{p} _{-}}\left( \omega \right) \right) \frac{\left\vert \hat{\rho}_{\mathrm{p} _{+},\omega }\left( \omega \right) \right\vert +\left\vert \hat{\rho}_{ \mathrm{p}_{-},\omega }\left( \omega \right) \right\vert }{2}.$$
A natural generalization of $\left \vert \vec {\tau }_{\omega \perp }\right \vert$, (29)$$T_{\mathrm{p}}\left( \omega \right) =\left( \tau _{\mathrm{p}_{1}}\left( \omega \right) -\tau _{\mathrm{p}_{J}}\left( \omega \right) \right) \frac{ \left\vert \hat{\rho}_{\mathrm{p}_{1},\omega }\left( \omega \right) \right\vert +\left\vert \hat{\rho}_{\mathrm{p}_{2},\omega }\left( \omega \right) \right\vert +\cdots+\left\vert \hat{\rho}_{\mathrm{p}_{J},\omega }\left( \omega \right) \right\vert }{J},$$
can be employed to measure the TPSC in modal spaces of arbitrary dimension $J$.4. Concatenation rules
Normally, to solve the fundamental equations in Sec. 3 for the evolution of $\hat {\Omega }$ and $\hat {\Omega }_{\omega }$, the fiber is divided into sufficiently small segments, with each segment $m$ possessing a constant (i.e., $z$-independent) $\hat {Z}\left ( z,\omega \right ) =\hat {Z}_{m}\left ( \omega \right )$, to convert the differential equations into difference equations. This approach can be regarded as a special implementation of the powerful concatenation rules that determine $\hat {\Omega }$ and $\hat {\Omega } _{\omega }$ for an assembly of concatenated fiber segments from those for the individual segments. While the concatenation rules in the PMD theory primarily concern the Stokes-space PMD vector $\vec {\tau }$ and its $\omega$-derivative $\vec {\tau }_{\omega }$, in our formalism, the concatenation rules will be formulated for the GD operator $\hat { \Omega }$ and its $\omega$-derivative $\hat {\Omega }_{\omega }$. We emphasize that the concatenation rules generally do not require the individual fiber segments be small or $\hat {Z}$ be constant within each segment.
According to Eq. (4), for two concatenated fiber segments from $z_{1}$ to $z_{2}$ and from $z_{2}$ to $z_{3}$, the “initial” states are transformed into the “final” states via
$$ \left\vert \psi \left( z_{2},\omega \right) \right\rangle =\hat{U}\left( z_{2},z_{1},\omega \right) \left\vert \psi \left( z_{1},\omega \right) \right\rangle ,\;\left\vert \psi \left( z_{3},\omega \right) \right\rangle = \hat{U}\left( z_{3},z_{2},\omega \right) \left\vert \psi \left( z_{2},\omega \right) \right\rangle , $$
(30)$$\left\vert \psi \left( z_{3},\omega \right) \right\rangle =\hat{U}\left( z_{3},z_{2},\omega \right) \hat{U}\left( z_{2},z_{1},\omega \right) \left\vert \psi \left( z_{1},\omega \right) \right\rangle .$$
The larger fiber segment from $z_{1}$ to $z_{3}$ has its own evolution operator $\hat {U}\left ( z_{3},z_{1},\omega \right )$ to connect the “initial” and “final” states (31)$$\left\vert \psi \left( z_{3},\omega \right) \right\rangle =\hat{U}\left( z_{3},z_{1},\omega \right) \left\vert \psi \left( z_{1},\omega \right) \right\rangle .$$
Since $\left \vert \psi \left ( z_{1},\omega \right ) \right \rangle$ can be any ket vector, comparing Eqs. (30) and (31) immediately leads to (32)$$\hat{U}\left( z_{3},z_{1},\omega \right) =\hat{U}\left( z_{3},z_{2},\omega \right) \hat{U}\left( z_{2},z_{1},\omega \right) .$$
Inserting Eq. (32) into Eq. (11) and taking into account Eq. (9), one acquires the basic (two-segment) concatenation rule for $\hat {\Omega },$ (33)$$\hat{\Omega}\left( z_{3},z_{1},\omega \right) =\hat{\Omega}\left( z_{3},z_{2},\omega \right) +\hat{U}\left( z_{3},z_{2},\omega \right) \hat{ \Omega}\left( z_{2},z_{1},\omega \right) \hat{U}^{\dagger }\left( z_{3},z_{2},\omega \right) ,$$
Furthermore, the basic concatenation rule for $\hat {\Omega }_{\omega }$ can be obtained by differentiating Eq. (33) with respective to $\omega$ and applying Eqs. (9) and (23), (34)$$\begin{aligned} \hat{\Omega}_{\omega }\left( z_{3},z_{1},\omega \right) &=\hat{\Omega} _{\omega }\left( z_{3},z_{2},\omega \right) +\hat{U}\left( z_{3},z_{2},\omega \right) \hat{\Omega}_{\omega }\left( z_{2},z_{1},\omega \right) \hat{U}^{\dagger }\left( z_{3},z_{2},\omega \right)\\ &\quad-i\left[ \hat{\Omega}\left( z_{3},z_{2},\omega \right) ,\hat{\Omega}\left( z_{3},z_{1},\omega \right) \right] . \end{aligned}$$
The concatenation rules for an assembly of $M$ concatenated segments, from $z_{0}$ to $z_{M}$ and divided at $z_{1},z_{2},\ldots ,z_{M-1}$, are reached by repeating the two-segment rules. For the evolution operator $\hat {U}$, repeated application of Eq. (32) yields
(35)$$\hat{U}\left( z_{M},z_{0},\omega \right) =\hat{U}\left( z_{M},z_{M-1},\omega \right) \hat{U}\left( z_{M-1},z_{M-2},\omega \right) \cdots\hat{U}\left( z_{1},z_{0},\omega \right) .$$
For the GD operator $\hat {\Omega }$, by repeatedly applying Eq. (33), one acquires the multiple-segment cancatenation rule (36)$$\begin{aligned} \hat{\Omega}\left( z_{M},z_{0},\omega \right) &=\sum_{m=1}^{M}\hat{W}_{m} \hat{\Omega}\left( z_{M-m+1},z_{M-m},\omega \right) \hat{W}_{m}^{\dagger },\\ \hat{W}_{m} &=\hat{U}\left( z_{M},z_{M-1},\omega \right) \hat{U}\left( z_{M-1},z_{M-2},\omega \right) \cdots\hat{U}\left( z_{M-m+2},z_{M-m+1},\omega \right) \end{aligned}$$
By use of Eq. (35), Eq. (36) can be formally simplified to (37)$$\hat{\Omega}\left( z_{M},z_{0},\omega \right) =\sum_{m=1}^{M}\hat{U}\left( z_{M},z_{M-m+1},\omega \right) \hat{\Omega}\left( z_{M-m+1},z_{M-m},\omega \right) \hat{U}^{\dagger }\left( z_{M},z_{M-m+1},\omega \right) ,$$
which further turns into an integral form in the continuous limit (see Appendix B) (38)$$\hat{\Omega}\left( z,z_{\mathrm{I}},\omega \right) =\int_{z_{\mathrm{I} }}^{z}dz^{\prime }\;\hat{U}\left( z,z^{\prime },\omega \right) \hat{Z} _{\omega }\left( z^{\prime },\omega \right) \hat{U}^{\dagger }\left( z,z^{\prime },\omega \right) .$$
Taking the $\omega$-derivative of Eqs. (37) and (38 ) yields the corresponding formulae for $\hat {\Omega }_{\omega }$, (39)$$\begin{aligned} \hat{\Omega}_{\omega }\left( z_{M},z_{0},\omega \right) &=\sum_{m=1}^{M} \hat{U}\left( z_{M},z_{M-m+1},\omega \right) \hat{\Omega}_{\omega }\left( z_{M-m+1},z_{M-m},\omega \right) \hat{U}^{\dagger }\left( z_{M},z_{M-m+1},\omega \right)\\ &\quad-i\sum_{m=1}^{M}\left[ \hat{\Omega}\left( z_{M},z_{M-m+1},\omega \right) , \hat{\Omega}\left( z_{M},z_{M-m},\omega \right) \right] , \end{aligned}$$
and (40)$$\begin{aligned} \hat{\Omega}_{\omega }\left( z,z_{\mathrm{I}},\omega \right) &=\int_{z_{ \mathrm{I}}}^{z}dz^{\prime }\;\hat{U}\left( z,z^{\prime },\omega \right) \hat{Z}_{\omega \omega }\left( z^{\prime },\omega \right) \hat{U}^{\dagger }\left( z,z^{\prime },\omega \right)\\ &\quad-i\int_{z_{\mathrm{I}}}^{z}dz^{\prime }\;\left[ \hat{\Omega}\left( z,z^{\prime },\omega \right) ,\hat{U}\left( z,z^{\prime },\omega \right) \hat{Z}_{\omega }\left( z^{\prime },\omega \right) \hat{U}^{\dagger }\left( z,z^{\prime },\omega \right) \right] , \end{aligned}$$
where Eqs. (9) and (23) have again been applied.5. Statistical analysis
In the preceding sections, we have shown that the density matrix operator $\hat {\rho }$, GD operator $\hat {\Omega }$ and its $\omega$-derivative $\hat {\Omega }_{\omega }$ evolve over $z$ according to Eqs. (17), (21) and (25). Alternatively, the evolution of $\hat {\Omega }$ and $\hat {\Omega }_{\omega }$ might be more conveniently achieved via the concatenation rules (33) and (34). For practical purposes (e.g., numerical simulations), we now use matrices in the free-eigenmode basis $\left \{ \left \vert j\right \rangle _{f},j=1,2,\ldots ,J\right \}$ to represent the operators. If the fiber is under random perturbation, the MP matrix $\bar {Z}$ randomly changes over $z,$
(41)$$\begin{aligned} \bar{Z}\left( z,\omega \right) &=\bar{Z}_{\mathrm{ID}}\left( z,\omega \right) +\Delta \bar{Z}\left( z,\omega \right) ,\\ \bar{Z}_{\mathrm{ID}}\left( z,\omega \right) &=\frac{\omega }{c}\bar{N}_{ \mathrm{ID}}\left( \omega \right) ,\;\Delta \bar{Z}\left( z,\omega \right) = \frac{\omega }{c}\Delta \bar{N}\left( z\right) , \end{aligned}$$
where $\bar {Z}_{\mathrm {ID}}\left ( z,\omega \right )$ and $\Delta \bar {Z} \left ( z,\omega \right )$ are respectively the ideal-fiber MP matrix and the random variation due to perturbation, with the normalized parts $\bar {N} _{\mathrm {ID}}\left ( \omega \right )$ independent of $z$ and $\Delta \bar {N} \left ( z\right )$ independent of $\omega$. In our statistical model, we employ the Langevin differential equation to treat $\Delta \bar {N}$, (42)$$\frac{d}{dz}\Delta \bar{N}\left( z\right) =-L_{C}^{-1}\left\{ \Delta \bar{N} \left( z\right) +\sqrt{L_{C}}\bar{g}\left( z\right) \right\} ,\;\Delta \bar{N }\left( z_{\mathrm{I}}\right) =0,$$
with $L_{C}$ as the correlation length. The Hermitian matrix $\bar {g}\left ( z\right )$ is a white-noise stochastic process satifying $\left \langle \bar {g}\left ( z\right ) \right \rangle =0$ and (43)$$\begin{aligned} &\left\langle \bar{g}_{jj}\left( z\right) \bar{g}_{jj}\left( z^{\prime }\right) \right\rangle =\gamma _{jj}^{2}\delta \left( z-z^{\prime }\right) ,\;\left\langle \textrm{Re}\left( \bar{g}_{jj^{\prime }}\left( z\right) \right) \textrm{Re}\left( \bar{g}_{jj^{\prime }}\left( z^{\prime }\right) \right) \right\rangle =\gamma _{jj^{\prime }\left( \textrm{Re}\right) }^{2}\delta \left( z-z^{\prime }\right) ,\\ \; &\left\langle \textrm{Im}\left( \bar{g}_{jj^{\prime }}\left( z\right) \right) \textrm{Im}\left( \bar{g}_{jj^{\prime }}\left( z^{\prime }\right) \right) \right\rangle =\gamma _{jj^{\prime }\left( \textrm{Im}\right) }^{2}\delta \left( z-z^{\prime }\right) ,\;\;j<j^{\prime }. \end{aligned}$$
Moreover, the entries on and above the main diagonal of matrix $\bar {g}$ (i.e., $\bar {g}_{j<j^{\prime }}$), as well as the real and imaginary parts of the same entry $\bar {g}_{j<j^{\prime }}$, are all statistically independent from each other. Meanwhile, we no longer regard the entries below the main diagonal as independent variables since $\bar {g} _{j^{\prime }j}=\bar {g}_{jj^{\prime }}^{\ast }$. Once $\bar {N}_{\mathrm {ID} }\left ( \omega \right )$ and the probability distributions of $\bar {g} _{j\geq j^{\prime }}$ are specified, the statistical properties of modal coupling and dispersion can be deduced via the evolution equations and/or concatenation rules. It is worth noting that the Langevin equation has also been employed in PMD statistics. For example, one of the models in Ref. [30] assumed that $\beta _{1}$ and $\beta _{2}$ of the birefrigence vector $\vec {\beta }$ are independent Langevin processes and $\beta _{3}=0,$ which is equivalent to the assumption that $\Delta \bar {N}_{11}-\Delta \bar {N}_{22}$ and Re$\left ( \Delta \bar {N}_{12}\right )$ are independent Langevin processes and Im$\left ( \Delta \bar {N}_{12}\right ) =0$. This is similar to Eq. (42) [with $\gamma _{jj^{\prime }\left ( \textrm {Im}\right ) }^{2}=0$] for $J=2$, with a significant difference that $\Delta \bar {N}_{11}$ and $\Delta \bar {N}_{22}$ vary independently according to Eq. (42).To simulate the stochastic evolution of $\bar {\rho }\left ( z,\omega \right )$, $\bar {\Omega }\left ( z,z_{\mathrm {I}},\omega \right )$ and $\bar {\Omega } _{\omega }\left ( z,z_{\mathrm {I}},\omega \right )$, we divide the fiber (from initial point $z_{\mathrm {I}}$ to final point $z_{\mathrm {F}}$) into segments of length $\delta z\ll L_{C},$ such that $\Delta \bar {N}$ can be treated as constant within each segment, i.e., $\Delta \bar {N}\left ( z\right ) =\Delta \bar {N}\left ( z_{m-1}^{+}\right )$ for the $m$th segment $z\in \left [ z_{m-1},z_{m}\right )$ ($z_{m}-z_{m-1}=\delta z$, $m=1,2,\ldots M$, $z_{0}=z_{\mathrm {I}},z_{M}=z_{\mathrm {F}}$, and $z_{m}^{+}$ means approaching $z_{m}$ from the right). In our current model, Re$\left ( \bar {g} _{j<j^{\prime }}\right )$, Im$\left ( \bar {g}_{j<j^{\prime }}\right )$ and $\bar {g}_{jj}$ each randomly take a value of $\frac {\gamma _{\cdots }}{ \sqrt {\delta z}}$ or $-\frac {\gamma _{\cdots }}{\sqrt {\delta z}}$ with equal probability $p=\frac {1}{2}$, where $\gamma _{\cdots }$ refer to the corresponding $\gamma$ in Eq. (43). Now we can acquire $\Delta \bar {N}\left ( z_{m}^{+}\right )$, $\bar {\rho }\left ( z_{m},\omega \right )$, $\bar {\Omega }\left ( z_{m},z_{\mathrm {I}},\omega \right )$ and $\bar {\Omega } _{\omega }\left ( z_{m},z_{\mathrm {I}},\omega \right )$ from $\Delta \bar {N} \left ( z_{m-1}^{+}\right )$, $\bar {\rho }\left ( z_{m-1},\omega \right )$, $\bar {\Omega }\left ( z_{m-1},z_{\mathrm {I}},\omega \right )$ and $\bar {\Omega } _{\omega }\left ( z_{m-1},z_{\mathrm {I}},\omega \right )$, as elaborated next.
By independently generating the random values for Re$\left ( \bar {g} _{j<j^{\prime }}\right )$, Im$\left ( \bar {g}_{j<j^{\prime }}\right )$ and $\bar {g}_{jj}$, $\Delta \bar {N}\left ( z_{m}^{+}\right )$ is obtained from $\Delta \bar {N}\left ( z_{m-1}^{+}\right )$ via Eq. (42),
$$ \Delta \bar{N}\left(z_{m}^{+}\right)=\Delta \bar{N}\left(z_{m-1}^{+}\right)-L_{C}^{-1}\left\{\Delta \bar{N}\left(z_{m-1}^{+}\right)+\sqrt{L_{c}} \bar{g}\right\} \delta z $$
$$ \bar{\rho}\left( z_{m},\omega \right) =\bar{U}\left( z_{m},z_{m-1},\omega \right) \bar{\rho}\left( z_{m-1},\omega \right) \bar{U}^{\dagger }\left( z_{m},z_{m-1},\omega \right) , $$
$$ \bar{U}\left( z_{m},z_{m-1},\omega \right) =e^{-i\bar{Z}\left( z_{m-1}^{+},\omega \right) \delta z}=e^{-i\frac{\omega }{c}\left\{ \bar{N}_{ \mathrm{ID}}\left( \omega \right) +\Delta \bar{N}\left( z_{m-1}^{+}\right) \right\} \delta z} $$
$$ \bar{\Omega}\left( z_{m},z_{\mathrm{I}},\omega \right) =\bar{\Omega}\left( z_{m},z_{m-1},\omega \right) +\bar{U}\left( z_{m},z_{m-1},\omega \right) \bar{\Omega}\left( z_{m-1},z_{\mathrm{I}},\omega \right) \bar{U}^{\dagger }\left( z_{m},z_{m-1},\omega \right) . $$
$$ \bar{\Omega}\left( z_{m},z_{m-1},\omega \right) =i\frac{\bar{U}\left( z_{m},z_{m-1},\omega +\delta \omega \right) -\bar{U}\left( z_{m},z_{m-1},\omega \right) }{\delta \omega }\bar{U}^{\dagger }\left( z_{m},z_{m-1},\omega \right) , $$
(44)$$\begin{aligned} &\bar{U}\left( z_{m},z_{m-1},\omega +\delta \omega \right) =\exp \left\{ -i\left( \bar{Z}\left( z_{m-1}^{+},\omega \right) +\bar{Z}_{\omega }\left( z_{m-1}^{+},\omega \right) \delta \omega +\frac{1}{2}\bar{Z}_{\omega \omega }\left( z_{m-1}^{+},\omega \right) \left( \delta \omega \right) ^{2}\right) \delta z\right\}\\ &=\exp \left\{ -i\frac{\omega }{c}\left( \mathcal{\bar{N}}_{0}+\mathcal{ \bar{N}}_{1}\delta \omega +\mathcal{\bar{N}}_{2}\left( \delta \omega \right) ^{2}\right) \delta z\right\} \end{aligned}$$
where $$ \mathcal{\bar{N}}_{0}=\bar{N}_{\mathrm{ID}}\left( \omega \right) +\Delta \bar{N}\left( z_{m-1}^{+}\right) ,\;\mathcal{\bar{N}}_{1}=\frac{1}{\omega } \left( \bar{N}_{\mathrm{ID}}^{\left( \mathrm{G}\right) }\left( \omega \right) +\Delta \bar{N}\left( z_{m-1}^{+}\right) \right) ,\;\mathcal{\bar{N}} _{2}=\frac{1}{2\omega }\bar{N}_{\mathrm{ID},\omega }^{\left( \mathrm{G} \right) }\left( \omega \right) , $$
$$\begin{aligned} \bar{\Omega}_{\omega }\left( z_{m},z_{\mathrm{I}},\omega \right) &=\bar{ \Omega}_{\omega }\left( z_{m},z_{m-1},\omega \right) +\bar{U}\left( z_{m},z_{m-1},\omega \right) \bar{\Omega}_{\omega }\left( z_{m-1},z_{\mathrm{ I}},\omega \right) \bar{U}^{\dagger }\left( z_{m},z_{m-1},\omega \right) \\ &-i\left[ \bar{\Omega}\left( z_{m},z_{m-1},\omega \right) ,\bar{\Omega} \left( z_{m},z_{\mathrm{I}},\omega \right) \right] , \end{aligned}$$
$$\begin{aligned} &\bar{\Omega}_{\omega }\left( z_{m},z_{m-1},\omega \right) =i\frac{\partial }{\partial \omega }\left\{ \bar{U}_{\omega }\left( z_{m},z_{m-1},\omega \right) \bar{U}^{\dagger }\left( z_{m},z_{m-1},\omega \right) \right\} \\ &\quad=i\bar{U}_{\omega \omega }\left( z_{m},z_{m-1},\omega \right) \bar{U} ^{\dagger }\left( z_{m},z_{m-1},\omega \right) +i\bar{U}_{\omega }\left( z_{m},z_{m-1},\omega \right) \bar{U}_{\omega }^{\dagger }\left( z_{m},z_{m-1},\omega \right) , \end{aligned}$$
(45)$$\begin{aligned} &\bar{U}_{\omega \omega }\left( z_{m},z_{m-1},\omega \right) \equiv \frac{ \partial ^{2}}{\partial \omega ^{2}}\bar{U}\left( z_{m},z_{m-1},\omega \right)\\ &=\frac{\bar{U}\left( z_{m},z_{m-1},\omega +2\delta \omega \right) -2\bar{U} \left( z_{m},z_{m-1},\omega +\delta \omega \right) +\bar{U}\left( z_{m},z_{m-1},\omega \right) }{\left( \delta \omega \right) ^{2}}. \end{aligned}$$
With the initial conditions $\Delta \bar {N}\left ( z_{\mathrm {I}}^{+}\right ) =0$, $\bar {\rho }\left ( z_{\mathrm {I}},\omega \right ) =\bar {\rho }_{0}$, $\bar { \Omega }\left ( z_{\mathrm {I}},z_{\mathrm {I}},\omega \right ) =0$, $\bar {\Omega } _{\omega }\left ( z_{\mathrm {I}},z_{\mathrm {I}},\omega \right ) =0$, and applying the equations in the previous paragraph for $m=1,2,\ldots M$ successively, we would eventually reach $\bar {\rho }\left ( z_{\mathrm {F} },\omega \right )$, $\bar {\Omega }\left ( z_{\mathrm {F}},z_{\mathrm {I}},\omega \right )$, $\bar {\Omega }_{\omega }\left ( z_{\mathrm {F}},z_{\mathrm {I} },\omega \right )$, and complete one implementation (sample) of the stochastic evolution process. Repeating the procedure for a large number ($Q$) of times, an ensemble of $Q$ samples are produced,
$$ \left\{ \bar{\rho}^{\left( q\right) }\left( z_{m=1,2,\ldots,M},\omega \right) , \bar{\Omega}^{\left( q\right) }\left( z_{m=1,2,\ldots ,M},z_{\mathrm{I}},\omega \right) ,\bar{\Omega}_{\omega }^{\left( q\right) }\left( z_{m=1,2,\ldots ,M},z_{ \mathrm{I}},\omega \right) ,q=1,2,\ldots ,Q\right\} . $$
(46)$$\bar{\rho}\left( z_{m},\omega \right) =\frac{1}{Q}\sum_{q=1}^{Q}\bar{\rho} ^{\left( q\right) }\left( z_{m},\omega \right) ,$$
which contains all observable information of the field at each $z_{m}$ during propagation. Particularly, the modal-coupling information, such as the energy distribution and coherence among a complete set of orthogonal modes, can be extracted from $\bar {\rho }\left ( z_{m=1,2,\ldots ,M},\omega \right )$ [see Eq. (15)]. For the modal dispersion, one needs to calculate relevant quantities (e.g., $\tau _{\mathrm {p}_{j}}$, $\tau _{ \mathrm {p}_{j},\omega }$, $\left \vert \hat {\rho }_{\mathrm {p}_{j},\omega }\right \vert$, and $T_{\mathrm {p}}$ in Sec. 3) from each $\bar {\Omega } ^{\left ( q\right ) }$ and $\bar {\Omega }_{\omega }^{\left ( q\right ) }$ in the ensemble via Eqs. (27) – (29), and then evaluate the statistical moments or probability density functions of the quantities over the ensemble.6. Numerical simulations
As a typical application of the density matrix formalism developed above, we now numerically simulate the statistical properties of coupling and dispersion in a $2$-dimensional ($2$D) and a $4$-dimensional ($4$D) modal space, respectively, with the $2$D case mainly serving the purpose of comparison to the PMD theory (see Appendix A for the correspondance of the density-matrix and Stokes-vector formalisms in the 2D modal space). Generally, for a modal space of dimension $J$, the ideal-fiber normalized MP matrix has a simple diagonal form in the free-eigenmode basis $\left \{ \left \vert j\right \rangle _{f},j=1,2,\ldots ,J\right \} ,$
$$ \bar{N}_{\mathrm{ID}}\left( \omega \right) =\textrm{diag} \left( n_{1}^{\left( 0\right) }\left( \omega \right) ,n_{2}^{\left( 0\right) }\left( \omega \right) ,\ldots,n_{J}^{\left( 0\right) }\left( \omega \right) \right) , $$
$$ \bar{N}_{\mathrm{ID}}^{\left( \mathrm{G}\right) }\left( \omega \right) = \textrm{diag}\left( n_{\mathrm{G,}1}^{\left( 0\right) }\left( \omega \right) ,n_{\mathrm{G,}2}^{\left( 0\right) }\left( \omega \right) ,\ldots,n_{\mathrm{G,} J}^{\left( 0\right) }\left( \omega \right) \right) , $$
We first conduct the simulations for the $2$-dimensional modal space, and the parameters are selected as: vacuum wavelength $\lambda =1.55$ $\mu m$ ($\omega =\frac {2\pi }{\lambda }$), degenerate ideal-mode ERIs $n_{1}^{\left ( 0\right ) }=n_{2}^{\left ( 0\right ) }=1.47777$ and EGRIs $n_{\mathrm {G} ,j}^{\left ( 0\right ) }=n_{j}^{\left ( 0\right ) }$ (neglecting ideal-mode chromatic dispersion), correlation length $L_{C}=50$ $m$, total fiber length $L=z_{\mathrm {F}}-z_{\mathrm {I}}=100L_{C}$, and “random kick” factors
$$ \gamma _{11}=\gamma _{22}=\gamma _{12\left( \textrm{Re}\right) }=10^{-5},\;\gamma _{12\left( \textrm{Im}\right) }=0, $$
(47)$$\sqrt{\left\langle \Delta \tau ^{2}\left( z,z_{\mathrm{I}}\right) \right\rangle }=\sqrt{\frac{3\pi }{8}}\left\langle \Delta \tau \left( z,z_{ \mathrm{I}}\right) \right\rangle =\sqrt{2\left\langle \Delta \tau _{C}^{2}\right\rangle }\left( e^{-\frac{z-z_{\mathrm{I}}}{L_{C}}}+\frac{z-z_{ \mathrm{I}}}{L_{C}}-1\right) ^{\frac{1}{2}},$$
with $\left \langle \Delta \tau _{C}^{2}\right \rangle$ being the value of $\left \langle \Delta \tau ^{2}\right \rangle$ at $z-z_{\mathrm {I}}\simeq 1.2L_{C}$, and (48)$$\begin{aligned} &\mathcal{P}_{\Delta \tau }\left( x\geq 0\right) =\frac{32x^{2}}{\pi ^{2}\left\langle \Delta \tau \right\rangle ^{3}}\exp \left( -\frac{4}{\pi } \left( \frac{x}{\left\langle \Delta \tau \right\rangle }\right) ^{2}\right) , \\\end{aligned}$$
$$\begin{aligned}\mathcal{P}_{\Delta \tau _{\omega }}\left( x\right) =\frac{2}{\left\langle \Delta \tau \right\rangle ^{2}}\textrm{sech}^{2}\left( \frac{4x}{\left\langle \Delta \tau \right\rangle ^{2}}\right) ,\end{aligned}$$
(49)$$\begin{aligned} &\mathcal{P}_{\left\vert \vec{\tau}_{\omega \bot }\right\vert }\left( x\geq 0\right) =x\left( \frac{8}{\pi \left\langle \Delta \tau \right\rangle ^{2}} \right) ^{2}\int_{0}^{\infty }d\zeta \;J_{0}\left( \frac{8\zeta x}{\pi \left\langle \Delta \tau \right\rangle ^{2}}\right) \left( \textrm{sech}\zeta \right) \sqrt{\left( \zeta \textrm{tanh}\zeta \right) }. \end{aligned}$$
with $\mathcal {P}_{\Delta \tau }\left ( x<0\right ) =0$ and $\mathcal {P} _{\left \vert \vec {\tau }_{\omega \bot }\right \vert }\left ( x<0\right ) =0$.Next, we turn to the $4$-dimensional modal space, conducting simulations with the following parameters: vacuum wavelength $\lambda =1.55$ $\mu m$, ideal-mode ERIs $n_{1}^{\left ( 0\right ) }=$ $n_{2}^{\left ( 0\right ) }=1.47778$, $n_{3}^{\left ( 0\right ) }=n_{4}^{\left ( 0\right ) }=1.47777$ and EGRIs $n_{\mathrm {G},j}^{\left ( 0\right ) }=n_{j}^{\left ( 0\right ) }$, correlation length $L_{C}=50$ $m$, total fiber length $L=100L_{C}$, and “random kick” factors
$$\begin{aligned} \gamma _{11} &=\gamma _{22}=\gamma _{33}=\gamma _{44}=\gamma _{12\left( \textrm{Re}\right) }=\gamma _{34\left( \textrm{Re}\right) }=\gamma _{0}=10^{-5}, \\ \gamma _{13\left( \textrm{Re}\right) } &=\gamma _{14\left( \textrm{Re}\right) }=\gamma _{23\left( \textrm{Re}\right) }=\gamma _{24\left( \textrm{Re}\right) }=0.5\gamma _{0},\;\gamma _{j<j^{\prime }\left( \textrm{Im}\right) }=0. \end{aligned}$$
In Fig. 4, the energy distribution $P_{j}=\bar {\rho }_{jj}$ and the coherence $\left \vert C_{j\neq j^{\prime }}\right \vert =\left \vert \bar {\rho }_{{j}\neq j^{\prime }}\right \vert$ are plotted for the above two initial states. Like in the $2$D case (Fig. 1), the energy turns equally distributed among the free eigenmodes and the coherence (if existing initially) gets destroyed during propagation under the random perturbation. The mean value ($\left \langle \Delta \tau _{\max }\right \rangle$) and rms ($\sqrt { \left \langle \Delta \tau _{\max }^{2}\right \rangle }$) of the group-delay difference ($\Delta \tau _{\max }\equiv \tau _{\mathrm {p}_{1}}-\tau _{ \mathrm {p}_{4}}$) vs $z-z_{\mathrm {I}}$ are shown in Fig. 5, and we numerically find that, similar to $\sqrt {\left \langle \Delta \tau ^{2}\right \rangle }$ and $\left \langle \Delta \tau \right \rangle$ in Eq. (47), $\sqrt {\left \langle \Delta \tau _{\max }^{2}\right \rangle }$ and $\left \langle \Delta \tau _{\max }\right \rangle$ also vary essentially as $\propto \left ( e^{-\frac {z-z_{\mathrm {I}}}{L_{C}}}+\frac {z-z_{\mathrm {I}} }{L_{C}}-1\right ) ^{\frac {1}{2}}$, only with the ratio $\frac {\sqrt { \left \langle \Delta \tau _{\max }^{2}\right \rangle }}{\left \langle \Delta \tau _{\max }\right \rangle }$ appearing to be significantly less than $\sqrt { \frac {3\pi }{8}}$.
The simulated PDFs for $\Delta \tau _{\max }$, $\Delta \tau _{\max ,\omega }\equiv \frac {\partial }{\partial \omega }\Delta \tau _{\max }$, and $T_{ \mathrm {p}}$ [Eq. (29)] of the whole fiber are depicted in Fig. 6(a)-(c). Although currently no analytical formula has been rigorously achieved, via numerical fitting we show that the PDFs can be described by some modified version of Eqs. (48) – (49 ), i.e.,
(50)$$\begin{aligned} &\mathcal{P}_{\Delta \tau _{\max }}\left( x\geq 0\right) =\frac{32\alpha _{1}x^{2}}{\pi ^{2}\left\langle \Delta \tau _{\max }\right\rangle ^{3}}\exp \left( -\frac{4\alpha _{2}}{\pi }\left( \frac{x-\tau _{\mathrm{0}}}{ \left\langle \Delta \tau _{\max }\right\rangle }\right) ^{2}\right) , \end{aligned}\\$$
$$\begin{aligned}\mathcal{P}_{\Delta \tau _{\max ,\omega }}\left( x\right) =\frac{2}{\alpha _{3}^{2}\left\langle \Delta \tau _{\max }\right\rangle ^{2}}\textrm{sech} ^{2}\left( \frac{4x}{\alpha _{3}^{2}\left\langle \Delta \tau _{\max }\right\rangle ^{2}}\right) ,\\ &\mathcal{P}_{T_{\mathrm{p}}}\left( x\geq T_{\mathrm{1}}^{\left( 0\right) }\right)\\,\end{aligned}$$
(51)$$\begin{aligned} &=\alpha _{4}\left( x-T_{\mathrm{1}}^{\left( 0\right) }\right) \left( \frac{ 8}{\pi \left\langle \Delta \tau _{\max }\right\rangle ^{2}}\right) ^{2}\int_{0}^{\infty }d\zeta \;J_{0}\left( \frac{8\zeta \left( x-T_{\mathrm{2 }}^{\left( 0\right) }\right) }{\pi \alpha _{5}^{2}\left\langle \Delta \tau _{\max }\right\rangle ^{2}}\right) \left( \textrm{sech}\zeta \right) \sqrt{ \left( \zeta \textrm{tanh}\zeta \right) }, \end{aligned}$$
with $\mathcal {P}_{\Delta \tau _{\max }}\left ( x<0\right ) =0$, $\mathcal {P} _{T_{\mathrm {p}}}\left ( x<T_{\mathrm {1}}^{\left ( 0\right ) }\right ) =0$, and $$\begin{aligned} \alpha _{1} &=0.72,\;\alpha _{2}=10,\;\alpha _{3}=0.62,\;\alpha _{4}=0.40,\;\alpha _{5}=0.89, \\ \tau _{\mathrm{0}} &=72\textrm{ ps},\;T_{\mathrm{1}}^{\left( 0\right) }=1.8\times 10^{3}\textrm{ ps}^{2},T_{\mathrm{2}}^{\left( 0\right) }=4.6\times 10^{3}\textrm{ ps}^{2}. \end{aligned}$$
Finally, to get a deeper insight into the influence of the perturbation strength, we calculate $\left \langle \Delta \tau _{\max }\right \rangle$ of the whole fiber for various $\gamma _{0}$’s and show the result in Fig. 7 as the red-solid curve. One can see that, as $\gamma _{0}$ goes up from $0$, $\left \langle \Delta \tau _{\max }\right \rangle$ first decreases from the ideal-fiber value, $\left ( n_{1}^{\left ( 0\right ) }-n_{4}^{\left ( 0\right ) }\right ) \frac {L}{c}=167$ ps, to a minimum of $62$ ps at $\gamma _{0}=4.5\times 10^{-6}$, and then increases monotonically for larger $\gamma _{0}$. This can be understood as follows: Without perturbation ($\gamma _{0}=0$), light propagating in the ideal fiber accumulates an intrinsic (determinstic) $\Delta \tau _{\max ,0}$ of $167$ ps; for relatively small $\gamma _{0}$, random fluctuations of the principle states $\left \vert \mathrm {p}_{j}\right \rangle$ and group delays $\tau _{\mathrm {p}_{j}}$ “average out” some determinstic accumulation of $\Delta \tau _{\max }$ and reduce the overall $\left \langle \Delta \tau _{\max }\right \rangle$, and hence $\left \langle \Delta \tau _{\max }\right \rangle$ decreases as $\gamma _{0}$ goes up; for sufficiently large $\gamma _{0}$, the perturbation-induced random-walk accumulation of $\Delta \tau _{\max }$ overwhelms the deterministic part, such that the overall $\left \langle \Delta \tau _{\max }\right \rangle$ rises linearly with $\gamma _{0}$; for intermediate $\gamma _{0}$, the “average-out” and the random-walk “accumulation” compete to yield a minimum $\left \langle \Delta \tau _{\max }\right \rangle =62$ ps at $\gamma _{0}=4.5\times 10^{-6}$. The scenario can be compared to the situation of zero intrinsic $\Delta \tau _{\max ,0}$. For example, setting $n_{j}^{\left ( 0\right ) }=n_{\mathrm {G },j}^{\left ( 0\right ) }=1.47777$ for all $j=1,2,3,4$, and other parameters identical to those adopted above, we simulate $\left \langle \Delta \tau _{\max }\right \rangle$ vs $\gamma _{0}$ and plot the curve in Fig. 7 as the blue-dashed one, which exhibits an essentially linear increase of $\left \langle \Delta \tau _{\max }\right \rangle$ with $\gamma _{0}$. The random “average-out vs accumulation” competition mechanism opens an interesting possibility to optimize the modal dispersion of a long fiber ($L\gg L_{C}$) with a non-zero intrinsic $\Delta \tau _{\max }$ by manipulating the random perturbation strength.
Based on the random “average-out vs accumulation” understanding, together with numerical-fitting and physical intuition, we again find an analytical formula for $\left \langle \Delta \tau _{\max }\right \rangle$ as function of $\gamma _{0}$
(52)$$\left\langle \Delta \tau _{\max }\right\rangle =\sqrt{\tau _{\max ,0}^{2}\left( 1+\delta \right) ^{-1}\left( e^{-\left( \frac{\gamma _{0}}{ \sigma _{1}}\right) ^{2}}+\delta \right) +\left( \frac{\gamma _{0}}{\sigma _{2}}\right) ^{2}},$$
with $\delta =0.09,$ $\sigma _{1}=2.2\times 10^{-6},$ $\sigma _{2}=1.55\times 10^{-7}$ ps$^{-1}$. This analytical formula is also plotted in Fig. 7 as the black-dash-dotted curve, which well fits the simulated one (red-solid). We also numerically simulate $\left \langle \Delta \tau _{\max }\right \rangle$ vs $\gamma _{0}$ for various deterministic $\Delta \tau _{\max ,0}$, and the results (not presented here) verify the validity of Eq. (52), with $\delta ,$ $\sigma _{2}$ independent of $\Delta \tau _{\max ,0}$ and $\sigma _{1}$ proportional to $\Delta \tau _{\max ,0}$. It is easy to understand that $\left ( \frac {\gamma _{0}}{\sigma _{2}}\right ) ^{2}$ represents the pure random-walk accumulation and $\sigma _{2}$ should not depend on $\Delta \tau _{\max ,0}$. Indeed, we plot $\left \langle \Delta \tau _{\max }\right \rangle =\frac {\gamma _{0}}{\sigma _{2}}$ (with $\sigma _{2}=1.55\times 10^{-7}$ ps$^{-1}$) in Fig. 7 as the red asterisks, and it excellently agrees with the simulated curve (blue-dashed) for $\Delta \tau _{\max ,0}=0$. The exponential function (with $\sigma _{1}$ proportional to $\Delta \tau _{\max ,0}$) in Eq. (52) is qualitatvely consistent with the walking distance of an object randomly side-kicked to deviate from the given direction, and thus understandably represents the “average-out” part of $\left \langle \Delta \tau _{\max }\right \rangle .$ The existence of $\delta$ seems to indicate some residue of the average-out process, although the interpretation is not as clear as the previous two terms. From Eq. (52), we can show that the $\gamma _{0}$ which gives the minimum $\left \langle \Delta \tau _{\max }\right \rangle$ is determined by $e^{-\left ( \sigma _{1}^{-1}\gamma _{0,\min }\right ) ^{2}}=\left ( 1+\delta \right ) \sigma _{2}^{-2}\left ( \sigma _{1}\tau _{\max ,0}^{-1}\right ) ^{2}$. Further noting that $\sigma _{1}\tau _{\max ,0}^{-1}$ is independent of $\tau _{\max ,0}$, one has $\gamma _{0,\min }=2.05\sigma _{1}.$Since $\sigma _{1}\propto$ $\Delta \tau _{\max ,0}$, $\gamma _{0,\min }$ is also proportional to $\Delta \tau _{\max ,0}$. Therefore, once the $\gamma _{0,\min }$ for a certain $\Delta \tau _{\max ,0}$ is known (e.g., for $\Delta \tau _{\max ,0}=167$ ps, $\gamma _{0,\min }=4.5\times 10^{-6}$ in Fig. 7), the $\gamma _{0,\min }$ for other $\Delta \tau _{\max ,0}$ can be easily determined. Finally, we remark that a major distinction between the above $\left \langle \Delta \tau _{\max }\right \rangle$-optimization mechanism and the group-delay reduction via strong mode coupling in Refs. [24,31,32] is that the later mainly addresses the random average-out effect of the perturbation while our treatment above simultaneously takes into account the competition from the random-walk accumulation.7. Summary
We have proposed and developed a density-matrix formalism for modal coupling and dispersion in mode-division multiplexing communications systems. The density-matrix formalism is simple and straightforward in modal spaces of arbitrary dimension $J$. We established the differential evolution equations and/or concatenation rules for the density matrix operator $\hat {\rho }$ and other elements (e.g., $\hat {\Omega }$, $\hat {\Omega }_{\omega }$) that characterize essential modal properties. Moreover, we constructed a statistical model ready for numerical analysis on the stochastic propagation of light in randomly perturbed fibers, and numerically simulated the modal properties of a $4$-mode ($J=4$) fiber under random perturbation. We found via numerical treatment that the PDFs of some key modal-dispersion quantities can be described by analytical formulae generalized from the PMD theory, and it should be desirable in future study to rigorously derive these formulae from the density-matrix formalism. We further found that for a long fiber ($L\gg L_{C}$) with an intrinsic non-zero group-delay difference $\Delta \tau _{\max ,0}$, the “average-out” and random-walk “accumulation” effects of the perturbation would compete to determine the overall $\left \langle \Delta \tau _{\max }\right \rangle$ and yield a minimum $\left \langle \Delta \tau _{\max }\right \rangle$ at certain perturbation strength, which raises the interesting possibility to optimize the modal dispersion by manipulation of the random perturbation.
Appendix
A. Correspondance of density-matrix and Stokes-vector formalisms in 2D modal Space
Here, we formulate the connection between the density-matrix formalism and the Stokes-vector formalism in the $2$-dimensional ($2$D) modal space. To this end, we first introduce the Pauli operators in the $2$D Jones space
$$ \hat{\sigma}_{1}=\left\vert 1\right\rangle \left\langle 1\right\vert -\left\vert 2\right\rangle \left\langle 2\right\vert ,\;\hat{\sigma} _{2}=\left\vert 1\right\rangle \left\langle 2\right\vert +\left\vert 2\right\rangle \left\langle 1\right\vert ,\;\hat{\sigma}_{3}=-i\left\vert 1\right\rangle \left\langle 2\right\vert +i\left\vert 2\right\rangle \left\langle 1\right\vert , $$
with
$\left \{ \left \vert 1\right \rangle ,\left \vert 2\right \rangle \right \}$ being a given basis. The Pauli operators satisfy
(53)$$\sigma _{\mu }^{\dagger }=\hat{\sigma}_{\mu },\;\hat{\sigma}_{\mu }\hat{ \sigma}_{\nu }=\delta _{\mu \nu }\hat{I}+i\sum_{\xi =1}^{3}\varepsilon _{\mu \nu \xi }\hat{\sigma}_{\xi },\;\textrm{Tr}\left( \hat{\sigma}_{\mu }\right) =0,$$
where
$\delta _{\mu =\nu }=1$,
$\delta _{\mu \neq \nu }=0,$ and
$\varepsilon _{\mu \nu \xi }=1$ for
$\left ( \mu \nu \xi \right ) =$ even permutation of
$\left ( 123\right )$,
$\varepsilon _{\mu \nu \xi }=-1$ for
$\left ( \mu \nu \xi \right ) =$ odd permutation of
$\left ( 123\right )$, and
$\varepsilon _{\mu \nu \xi }=0$ for others. Any operator
$\hat {O}$ in the
$2$D Jones space can be expanded as the superposition of
$\hat {\sigma }_{\mu }$ and unit operator
$\hat {I},$ (54)$$\hat{O}=\frac{1}{2}\hat{I}\textrm{Tr}\left( \hat{O}\right) +\frac{1}{2} \sum_{\mu =1}^{3}\hat{\sigma}_{\mu }\textrm{Tr}\left( \hat{\sigma}_{\mu }\hat{O }\right) .$$
From Eqs. (2), (3) and (14), we can verify that the Stokes vector and the density matrix are connected via the Pauli operators as
(55)$$S_{\mu =1,2,3}=\textrm{Tr}\left( \hat{\rho}\hat{\sigma}_{\mu }\right) ,\;\hat{ \rho}=\frac{1}{2}\hat{I}+\frac{1}{2}\sum_{\mu =1}^{3}S_{\mu }\hat{\sigma} _{\mu },$$
and similarly, it follows from Eqs. (
6) and (
13) that the birefringence vector
$\vec {\beta }$ and the PMD vector can be written in terms of
$\hat {Z}$ and
$\hat {\Omega }$ as
(56)$$\beta _{\mu }=\textrm{Tr}\left( \hat{Z}\hat{\sigma}_{\mu }\right) ,\;\tau _{\mu }=\textrm{Tr}\left( \hat{\Omega}\hat{\sigma}_{\mu }\right) .$$
Furthermore, using Eqs. (
53), (
55) and (
56), we can cast Eqs. (
17) and (
21) into Eq. (
5) and
$$ \frac{\partial \vec{\tau}}{\partial z}=\frac{\partial \vec{\beta}}{\partial \omega }+\vec{\beta}\times \vec{\tau}, $$
respectively.
The depolarization in the PMD theory is characterized by the “perpendicular component” of $\vec {\tau } _{\omega }\equiv \frac {\partial }{\partial \omega }\vec {\tau }$ [8],
$$ \vec{\tau}_{\omega \perp }=\Delta \tau \frac{\partial }{\partial \omega } \vec{S}_{\mathrm{P}}^{\left( +\right) }, $$
where
$\vec {S}_{\mathrm {P}}^{\left ( +\right ) }\equiv \frac {\vec {\tau }}{ \left \vert \vec {\tau }\right \vert }$ is the Stokes vector corresponding to the principal state
$\left \vert \mathrm {p}_{+}\right \rangle$ of larger group delay
$\tau _{\mathrm {p}_{+}}$. According to Eq. (
55), one has
$\tau _{\omega \perp ,\mu =1,2,3}=\Delta \tau$Tr
$\left ( \hat {\rho }_{ \mathrm {p}_{+},\omega }\hat {\sigma }_{\mu }\right )$ [
$\hat {\rho }_{\mathrm {p} _{\pm },\omega }\equiv \frac {\partial }{\partial \omega }\hat { \rho }_{\mathrm {p}_{\pm },\omega }$,
$\hat {\rho }_{\mathrm {p}_{\pm }}\equiv \left \vert \mathrm {p}_{\pm }\right \rangle \left \langle \mathrm {p}_{\pm }\right \vert$], which further yields the modulus of
$\vec {\tau }_{\omega \perp }$ as
$$ \left\vert \vec{\tau}_{\omega \perp }\right\vert =\Delta \tau \sqrt{ \sum_{\mu =1}^{3}\tau _{\omega \perp ,\mu }\tau _{\omega \perp ,\mu }} =\Delta \tau \sqrt{\textrm{Tr}\left( \hat{\rho}_{\mathrm{p}_{+},\omega }\sum_{\mu =1}^{3}\hat{\sigma}_{\mu }\textrm{Tr}\left( \hat{\rho}_{\mathrm{p} _{+},\omega }\hat{\sigma}_{\mu }\right) \right) }=\Delta \tau \sqrt{2\textrm{Tr }\hat{\rho}_{\mathrm{p}_{+},\omega }^{2}}, $$
where, in the last step, we have applied Eq. (
54) and Tr
$\hat { \rho }_{\mathrm {p}_{+},\omega }=\frac {\partial }{\partial \omega }$Tr
$\hat { \rho }_{\mathrm {p}_{+}}=0$. Defining
$\left \vert \hat {\rho }_{\mathrm {p}_{\pm },\omega }\right \vert =\sqrt {2\textrm {Tr}\hat {\rho }_{\mathrm {p}_{\pm },\omega }^{2}}$, and noting
$\left \vert \hat {\rho }_{\mathrm {p}_{-},\omega }\right \vert =\left \vert \frac {\partial }{\partial \omega }\left ( \hat {I}- \hat {\rho }_{\mathrm {p}_{+}}\right ) \right \vert =\left \vert \hat {\rho }_{ \mathrm {p}_{+},\omega }\right \vert$, we can re-write
$\left \vert \vec {\tau } _{\omega \perp }\right \vert$ into an alternative form that is more suitable for generalization to higher-dimensional modal spaces,
$$ \left\vert \vec{\tau}_{\omega \perp }\right\vert =\Delta \tau \frac{ \left\vert \hat{\rho}_{\mathrm{p}_{+},\omega }\right\vert +\left\vert \hat{ \rho}_{\mathrm{p}_{-},\omega }\right\vert }{2}. $$
B. Derivation of Eq. (38)
To derive Eq. (38), we divide the fiber $\left [ z_{\mathrm {I}},z \right ]$ into small segments $\left [ z_{m},z_{m+1}\right ]$ ($m=0,1,2,\ldots ,M-1$), with $z_{0}=z_{\mathrm {I}}$, $z_{M}=z$, and $\delta z=z_{m+1}-z_{m}$. From Eqs. (32) and (9), it follows that
$$\begin{aligned} &\hat{U}\left( z_{m+1},z_{m},\omega \right) =\hat{U}\left( z_{m+1},z_{0},\omega \right) \hat{U}^{\dagger }\left( z_{m},z_{0},\omega \right) \\ &=\hat{U}\left( z_{m+1},z_{0},\omega \right) \left( \hat{U}^{\dagger }\left( z_{m+1},z_{0},\omega \right) -U_{z}^{\dagger }\left( z_{m+1},z_{0},\omega \right) \delta z\right) =\hat{I}-i\hat{Z}\left( z_{m+1},\omega \right) \delta z, \\ &\hat{U}_{\omega }\left( z_{m+1},z_{m},\omega \right) =\frac{\partial }{ \partial \omega }\left( \hat{I}-i\hat{Z}\left( z_{m+1},\omega \right) \delta z\right) =-i\hat{Z}_{\omega }\left( z_{m+1},\omega \right) \delta z, \end{aligned}$$
where we have used Eqs. (
19) and (
22), and neglected terms of second or higher orders in
$\delta z$. Insertion of the above equations into Eq. (
11) leads to
$$ \hat{\Omega}\left( z_{m+1},z_{m},\omega \right) =i\left( -i\hat{Z}_{\omega }\left( z_{m+1},\omega \right) \delta z\right) \left( \hat{I}+i\hat{Z}\left( z_{m+1},\omega \right) \delta z\right) =\hat{Z}_{\omega }\left( z_{m+1},\omega \right) \delta z, $$
which we substitute into Eq. (
37) to get (with
$m^{\prime }=M-m+1$, and noting
$z=z_{M}$,
$z_{\mathrm {I}}=z_{0}$)
$$\begin{aligned} \hat{\Omega}\left( z,z_{\mathrm{I}},\omega \right) &=\sum_{m^{\prime }=1}^{M}\delta z\hat{U}\left( z,z_{m^{\prime }},\omega \right) \hat{Z} _{\omega }\left( z_{m^{\prime }},\omega \right) \hat{U}^{\dagger }\left( z,z_{m^{\prime }},\omega \right) \\ &\quad\rightarrow \int_{z_{\mathrm{I}}}^{z}dz\;\hat{U}\left( z,z^{\prime },\omega \right) \hat{Z}_{\omega }\left( z^{\prime },\omega \right) \hat{U} ^{\dagger }\left( z,z^{\prime },\omega \right) \textrm{ \ }\left( \textrm{as } \delta z\rightarrow 0\right) . \end{aligned}$$
This completes the derivation of Eq. (
38).
Funding
Natural Science Foundation of Shandong Province (ZR2018MA044, ZR2013AL007); National Natural Science Foundation of China (61501214, 11404156).
Disclosures
The authors declare no conflicts of interest.
References
1. S. Berdague and P. Facq, “Mode Division Multiplexing in Optical Fibers,” Appl. Opt. 21(11), 1950–1955 (1982). [CrossRef]
2. R. Ryf, S. Randel, A. H. Gnauck, C. Bolle, A. Sierra, S. Mumtaz, M. Esmaeelpour, E. C. Burrows, R. J. Essiambre, P. J. Winzer, and D. W. Peckham, “Mode-Division Multiplexing Over 96 km of Few-Mode Fiber Using Coherent 6×6 MIMO Processing,” J. Lightwave Technol. 30(4), 521–531 (2012). [CrossRef]
3. L. W. Luo, N. Ophir, C. P. Chen, L. H. Gabrielli, C. B. Poitras, K. Bergmen, and M. Lipson, “WDM-Compatible Mode-Division Multiplexing on a Silicon Chip,” Nat. Commun. 5(1), 3069 (2014). [CrossRef]
4. D. Zhou, C. Sun, Y. Lai, Y. Yu, and X. Zhang, “Integrated Silicon Multifunctional Mode-Division Multiplexing System,” Opt. Express 27(8), 10798–10805 (2019). [CrossRef]
5. G. J. Foschini and C. D. Poole, “Statistical Theory of Polarization Dispersion in Single Mode Fibers,” J. Lightwave Technol. 9(11), 1439–1456 (1991). [CrossRef]
6. J. P. Gordon and H. Kogelnik, “PMD Fundamentals: Polarization Mode Dispersion in Optical Fibers,” Proc. Natl. Acad. Sci. U. S. A. 97(9), 4541–4550 (2000). [CrossRef]
7. H. Kogelnik, R. M. Jopson, and L. E. Nelson, “Polarization-Mode Dispersion,” in Optical Fiber Telecommunications IV-B (Academic, 2002), pp. 725–861
8. J. N. Damask, Polarization Optics in Telecommunications Vol. 101 (Springer Science & Business Media, 2004).
9. K. P. Ho and J. M. Kahn, “Statistics of Group Delays in Multimode Fiber with Strong Mode Coupling,” J. Lightwave Technol. 29(21), 3119–3128 (2011). [CrossRef]
10. K. P. Ho and J. M. Kahn, “Mode Coupling and its Impact on Spatially Multiplexed Systems,” in Optical Fiber Telecommunications VI, I. P. Kaminow, T. Li, and A. E. Willner, eds. (Elsevier, 2013).
11. J. Wang, J. Y. Yang, I. M. Fazal, N. Ahmed, Y. Yan, H. Huang, Y. X. Ren, Y. Yue, S. Dolinar, M. Tur, and A. E. Willner, “Terabit Free-Space Data Transmission Employing Orbital Angular Momentum Multiplexing,” Nat. Photonics 6(7), 488–496 (2012). [CrossRef]
12. N. Bozinovic, Y. Yue, Y. Ren, M. Tur, P. Kristensen, H. Huang, A. E. Willner, and S. Ramachandran, “Terabit-scale Orbital Angular Momentum Mode Division Multiplexing in Fibers,” Science 340(6140), 1545–1548 (2013). [CrossRef]
13. R. M. Nejad, K. Allahverdyan, P. Vaity, S. Amiralizadeh, C. Brunet, Y. Messaddeq, S. LaRochelle, and L. A. Rusch, “Mode Division Multiplexing Using Orbital Angular Momentum Modes over 1.4-km Ring Core Fiber,” J. Lightwave Technol. 34(18), 4252–4258 (2016). [CrossRef]
14. J. Liu, G. X. Zhu, J. W. Zhang, Y. H. Wen, X. Wu, Y. F. Zhang, Y. J. Chen, X. L. Cai, Z. H. Li, Z. Y. Hu, J. B. Zhu, and S. Y. Yu, “Mode Division Multiplexing Based on Ring Core Optical Fibres,” IEEE J. Quantum Electron. 54(5), 0700118 (2018). [CrossRef]
15. J. Yang, H. Liu, F. Pang, J. Wen, H. Zheng, L. Chen, X. He, Y. Shang, N. Chen, Y. Li, and T. Wang, “All-Fiber Multiplexing and Transmission of High-Order Circularly Polarized Orbital Angular Momentum Modes with Mode Selective Couplers,” IEEE Photonics J. 11(3), 1–9 (2019). [CrossRef]
16. L. Allen, M. W. Beijersbergen, R. J. Spreeuw, and J. P. Woerdman, “Orbital Angular Momentum of Light and the Transformation of Laguerre-Gaussian Laser Modes,” Phys. Rev. A 45(11), 8185–8189 (1992). [CrossRef]
17. A. M. Yao and M. J. Padgett, “Orbital Angular Momentum: Origins, Behavior and Applications,” Adv. Opt. Photonics 3(2), 161–204 (2011). [CrossRef]
18. D. L. Andrews and M. Babiker, eds. The Angular Momentum of Light (Cambridge University, 2012).
19. C. Antonelli, A. Mecozzi, M. Shtaif, and P. J. Winzer, “Stokes-Space Analysis of Modal Dispersion in Fibers with Multiple Mode Transmission,” Opt. Express 20(11), 11718–11733 (2012). [CrossRef]
20. I. Roudas and J. Kwapisz, “Stokes Space Representation of Modal Dispersion,” IEEE Photonics J. 9(5), 7203715 (2017). [CrossRef]
21. G. M. Fernandes, N. J. Muga, and A. N. Pinto, “Digital Monitoring and Compensation of MDL Based on Higher-Order Poincaré Spheres,” Opt. Express 27(14), 19996–20011 (2019). [CrossRef]
22. C. Antonelli, A. Mecozzi, M. Shtaif, N. Fontaine, H. S. Chen, and R. Ryf, “Stokes-Space Analysis of Modal Dispersion of SDM Fibers with Mode-Dependent Loss: Theory and Experiments,” J. Lightwave Technol. 38(7), 1668–1677 (2020). [CrossRef]
23. K. P. Ho and J. M. Kahn, “Linear Propagation Effects in Mode-Division Multiplexing Systems,” J. Lightwave Technol. 32(4), 614–628 (2014). [CrossRef]
24. S. Ö. Arik, K. Ho, and J. M. Kahn, “Delay Spread Reduction in Mode-Division Multiplexing: Mode Coupling Versus Delay Compensation,” J. Lightwave Technol. 33(21), 4504–4512 (2015). [CrossRef]
25. P. A. M. Dirac, The Principles of Quantum Mechanics (Oxford University, 1981).
26. G. Milione, D. A. Nolan, and R. R. Alfano, “Determining Principal Modes in a Multimode Optical Fiber Using the Mode Dependent Signal Delay Method,” J. Opt. Soc. Am. B 32(1), 143–149 (2015). [CrossRef]
27. G. M. Fernandes, N. J. Muga, and A. N. Pinto, “Space-Demultiplexing Based on Higher-Order Poincaré Spheres,” Opt. Express 25(4), 3899 (2017). [CrossRef]
28. I. Roudas, J. Kwapisz, and D. A. Nolan, “Optimal Launch States for the Measurement of Principal Modes in Optical Fibers,” J. Lightwave Technol. 36(20), 4915–4931 (2018). [CrossRef]
29. A. Andrusier, M. Shtaif, C. Antonelli, and A. Mecozzi, “Assessing the Effects of Mode-Dependent Loss in Space-Division Multiplexed Systems,” J. Lightwave Technol. 32(7), 1317–1322 (2014). [CrossRef]
30. P. K. A. Wai and C. R. Menyuk, “Polarization Mode Dispersion, Decorrelation, and Diffusion in Optical Fibers with Randomly Varying Birefringence,” J. Lightwave Technol. 14(2), 148–157 (1996). [CrossRef]
31. S. Ö. Arik, K. P. Ho, and J. M. Kahn, “Group Delay Management and Multi-Input Multi-Output Signal Processing in Mode-Division Multiplexing Systems,” J. Lightwave Technol. 34(11), 2867–2880 (2016). [CrossRef]
32. F. Ferreira, D. Fonseca, A. Lobato, B. Inan, and H. Silva, “Reach Improvement of Mode Division Multiplexed Systems Using Fiber Splices,” IEEE Photonics Technol. Lett. 25(12), 1091–1094 (2013). [CrossRef]
Open Access
Add to BibSonomy
