---
abstract: 'We consider a -based transmission scheme, where data is embedded into the imaginary part of the nonlinear discrete spectrum. Inspired by probabilistic amplitude shaping, we propose a scheme as a means to increase the data rate of the system. We exploit the fact that for an -based transmission scheme, the pulses in the time domain are of unequal duration by transmitting them with a dynamic symbol interval and find a capacity-achieving distribution. The scheme shapes the information symbols according to the capacity-achieving distribution and transmits them together with the parity symbols at the output of a encoder, suitably modulated, via time-sharing. We furthermore derive an achievable rate for the proposed scheme. We verify our results with simulations of the discrete-time model as well as with simulations.'
author:
- |
    Andreas Buchberger, Alexandre Graell i Amat, \
    Vahid Aref,  and Laurent Schmalen,  [^1][^2][^3][^4]
title: |
    Probabilistic Eigenvalue Shaping for\
    Nonlinear Fourier Transform Transmission
---

Discrete spectrum, nonlinear Fourier transform (NFT), probabilistic shaping, soliton communication.

Introduction {#sec:introduction}
============

propagation in optical fibers is severely impaired by nonlinear effects that should be either compensated or utilized for the design of the communication system. The  [@yousefi2014nftI-III] provides a method to transform a signal from the time domain into a nonlinear frequency domain (spectrum), where the channel acts as a multiplicative filter on the signal. The nonlinear spectrum consists of a continuous and a discrete part. Both parts can be used to transmit information, either separately or jointly, and several schemes have been presented in theory and practice [@yousefi2014nftI-III; @dong2015; @aref2015ecoc; @aref2016ecoc; @geisler2016ecoc; @hari2016]. However, very little is known so far about the of the received signal in the nonlinear spectral domain when it is contaminated by channel noise.

In [@shevchenko2018tcom], a simplified communication system modulating only the imaginary part of the eigenvalues in the discrete nonlinear spectrum was presented. For this scheme, an approximation for the conditional of the channel can be obtained in closed form. In general, for a given channel, the capacity-achieving distribution is not known and is often different from the conventional distribution with equispaced signal points and uniform signaling. Hence, some form of shaping is required [@forney1984]. Two popular methods of shaping are probabilistic shaping and geometric shaping. In geometric shaping, the capacity-achieving distribution is mimicked by optimizing the position of the constellation points for equiprobable signaling [@sun1993_geometric_shaping] whereas probabilistic shaping uses uniformly spaced constellation points and approximates the capacity-achieving distribution by assigning different probabilities to different constellation points [@forney1984].

The main drawback of probabilistic shaping is its practical implementation. An abundance of probabilistic shaping schemes have been presented, most suffering from high decoding complexity, low flexibility in adapting the spectral efficiency, or error propagation. For a literature review on probabilistic shaping, we refer the reader to [@bocherer2015_bw_efficient_rate_matched_ldpc Section II].

Recently, a new scheme called has been proposed in [@bocherer2015_bw_efficient_rate_matched_ldpc]. Compared to other shaping schemes, yields high flexibility and close-to-capacity performance over a wide range of spectral efficiencies for the channel while still allowing bit-metric decoding. Although originally introduced for the channel, can be applied to other channels with a symmetric capacity-achieving input distribution assuming a sufficiently high spectral efficiency.

In this paper, we consider a similar -based transmission scheme to the one presented in [@shevchenko2018tcom], where data is embedded into the imaginary part of the nonlinear discrete spectrum. As a means to increase the data rate, we demonstrate that the concept of can be adapted to this -based transmission system. In particular, we propose a scheme, enabling similar low complexity and bit-metric decoding as . We take advantage of the dependence of the pulse length on the data for the -based transmission system and transmit each pulse as soon as the previous one has been transmitted rather than with a fixed interval as in [@shevchenko2018tcom], yielding increased data rate. Accordingly, we find the capacity-achieving input distribution, maximizing the time-scaled . For ease of notation, we refer to the maximized as capacity noting that it is in fact the constrained capacity of a system transmitting first-order solitons. The scheme then shapes the information symbols according to the capacity-achieving distribution by a . The information symbols are also encoded by a encoder and the parity symbols at the output of the encoder are suitably modulated. The resulting sequence of modulated symbols and the sequence at the output of the are transmitted via time-sharing. We further derive an achievable rate for such a scheme. We demonstrate via discrete-time Monte-Carlo and simulations, that performs at around $\SI{2}{\decibel}$ from capacity using off-the-shelf codes. The proposed scheme yields a significant improvement of up to twice the data rate compared to an unshaped system as in [@shevchenko2018tcom].

It is important to note that although first-order solitons do not outperform conventional coherent systems due to their spectrally inefficient pulse shape compared to a Nyquist pulse shape, they have some other advantages. For instance, the first-order soliton transmission does not require compensation or as dispersion and nonlinearity are balanced and hence compensated. This work attempts to approach the limits of current NFT-based systems. To improve the spectral efficiency further, one should use higher-order solitons as well as the continuous part of the nonlinear spectrum together [@aref2016ecoc]. However, the channel equalization will not be as easy as the one for the first-order solitons and the channel model is not yet fully known.

The remainder of the paper is organized as follows. In [Section ]{}, we describe pulse propagation in an optical fiber and the -based transmission scheme. In [Section ]{}, we optimize the input distribution and in [Section ]{}, we introduce and describe the proposed scheme and derive an achievable rate. In [Section ]{}, we present numerical results for , both from Monte-Carlo simulation and simulation, and in [Section ]{} we draw some conclusions.

Notation: The following notation is used throughout the paper. ${\mathfrak{R}\{\cdot\}}$ and ${\mathfrak{I}\{\cdot\}}$ denote the real and the imaginary part of a complex number, respectively, and $\jmath=\sqrt{-1}$ denotes the imaginary unit. Vectors are typeset in bold, e.g., $\bm{x}$, are capitalized, e.g., $X$, and hence vectors of are capitalized bold, e.g., $\bm{X}$. The of an $X$ is written as $p_X(x)$ and its expectation as ${\mathbb{E}_{X}\!{\mathopen{}\mathclose\bgroup\originalleft}\{{x}{\aftergroup\egroup\originalright}\}
}$. The conditional of $Y$ given $X$ is denoted as $p_{Y|X}(y|x)$. The of an $X$ is denoted by $P_X(x)$. The transpose of a vector or matrix is given as $(\cdot )^\mathsf{T}$. A set is denoted by a capitalized Greek letter, e.g., $\Lambda$, and its cardinality by $|\Lambda|$. We write $\log_a(\cdot)$ for the logarithm of base $a$ and $\ln(\cdot)$ for the natural logarithm.

Nonlinear Fourier Transform-based Transmission System {#sec:system}
=====================================================

Pulse Propagation and the Nonlinear Fourier Transform {#subsec:NFT}
-----------------------------------------------------

Pulse propagation in optical fibers is governed by a partial differential equation, the stochastic , $$\begin{gathered}
        \jmath\frac{\partial u(\tau,\ell)}{\partial \ell} + \jmath\frac{\alpha}{2}u(\tau,\ell) - \frac{\beta_2}{2}\frac{\partial^2 u(\tau,\ell)}{\partial \tau^2} + \gamma u(\tau,\ell) |u(\tau,\ell)|^2  \\= n(\tau,\ell ),\nonumber$$ where $u(\tau,\ell)$ denotes the envelope of the electrical field as a function of the position $\ell$ along the fiber and time $\tau$, $\alpha$ the attenuation, $\beta_2$ the second order dispersion, $\gamma$ the nonlinearity parameter, and $n(\tau, \ell)$ is a white Gaussian process in time and in space with spectral density $\sigma_0^2$. The spectral density depends on the system and for distributed Raman amplification is given as $\sigma_0^2 = \alpha K_T h \nu_0$, where $K_T$ is the temperature-dependent phonon occupancy factor, and $h\nu_0$ is the average photon energy [@shevchenko2018tcom]. A general closed-form solution of the stochastic does not exist. In some special cases, e.g., for noisefree and lossless fibers, special solutions like, e.g., solitons, exist. Furthermore, we consider the in normalized form in the focusing regime, i.e., $\beta_2 < 0$, under the assumption of ideal distributed Raman amplification, i.e., $\alpha = 0$, $$\begin{aligned}
    \jmath\frac{\partial q(t,z)}{\partial z}  + \frac{\partial^2 q(t,z)}{\partial t^2} + 2 q(t,z) |q(t,z)|^2  =  0,
    \label{eqn:nlse_losless_noisefree}\end{aligned}$$ where $t=\tau / \sqrt{|\beta_2| L/2}$, $z=\ell / L$, $q = u\sqrt{\gamma L}/\sqrt{2}$, and $L$ is the length of the fiber. In this case, the is an integrable partial differential equation for which a pair of operators, called Lax pair, can be found. The eigenvalues of such an operator remain invariant during noiseless propagation and the Lax pair can be used to solve the partial differential equation. Solutions of [()]{} can be uniquely represented in terms of its eigenvalues via the so-called . For a given position $z$, the of a signal $q(t)$ (we drop the position $z$ for simplicity of presentation) with support on the time interval $t\in [t_1,t_2]$, is calculated by solving the partial differential equation $$\begin{aligned}
\frac{\partial \bm{v}(t,\lambda)}{\partial t} = \begin{pmatrix}
-\jmath\lambda & q(t)\\ 
- q(t)^*&  \jmath\lambda
\end{pmatrix}\bm{v}(t,\lambda),\label{eqn:nft_dgl}\end{aligned}$$ where $\bm{v}(t,\lambda)=\begin{pmatrix} v_1(t,\lambda) & v_2(t,\lambda) \end{pmatrix}$ is the eigenvector of the auxiliary operator, with boundary conditions $$\begin{aligned}
\bm{v}^{(1)}(t,\lambda) &\rightarrow {\begin{pmatrix} 0 & 1 \end{pmatrix}^\mathsf{T}} {\operatorname{e}^{\jmath\lambda t}}, &\text{as } t &\rightarrow t_2\\
\bm{v}^{(2)}(t,\lambda) &\rightarrow {\begin{pmatrix} 1 & 0 \end{pmatrix}^\mathsf{T}} {\operatorname{e}^{-\jmath\lambda t}}, &\text{as } t &\rightarrow t_1,\end{aligned}$$ and $\lambda$ is the spectral component. Solving [()]{} gives rise to the continuous and discrete nonlinear spectrum $$\begin{aligned}
    \hat{q}(\lambda ) &= \frac{b(\lambda)}{a(\lambda)}, \lambda\in{\mathbb{R}}& \tilde{q}(\lambda_i) &= \frac{b(\lambda_i)}{\mathrm{d}a(\lambda) / \mathrm{d}\lambda|_{\lambda=\lambda_i}}, \lambda_i\in{\mathbb{C}}^{+},\end{aligned}$$ respectively, where $a(\lambda ) = \lim_{t\rightarrow t_2} v_1^{(2)}(t,\lambda)    {\operatorname{e}^{\jmath\lambda t}}$, $b(\lambda ) = \lim_{t\rightarrow t_2} v_2^{(2)}(t,\lambda)   {\operatorname{e}^{-\jmath\lambda t}}$, and $\lambda_i$ are the zeros of $a(\lambda)$, $\lambda_i\in{\mathbb{C}}^{+}$, a finite set of isolated complex zeros, referred to as eigenvalues. Hence, the represents the signal in the nonlinear spectral domain, where the influence of the channel on the signal is a multiplicative filter.

As a counterpart to the that transforms a signal from the time domain to the nonlinear spectral domain, the transforms a signal from the nonlinear spectral domain to the time domain. For an in-depth mathematical description of the , we refer the interested reader to [@yousefi2014nftI-III].

Soliton Transmission
--------------------

![Block diagram of the -based system.[]{data-label="fig:system:block_diagram_nft"}](block_diagram_nft)

As in [@shevchenko2018tcom], we embed information in the imaginary part of the discrete spectrum, also referred to as eigenvalues. Hence, the input of the channel is an $X\in\Lambda=\{\lambda_1,\ldots,\lambda_M\}$, where $\Lambda$ is the set of eigenvalues, $\lambda_i$ is the $i$th eigenvalue, and $M$ is the order of the modulation. The eigenvalues $\{\lambda_i\}$ are assumed to be ordered in ascending order by their imaginary parts. Furthermore, the output of the channel is an $Y\in\Psi$, where $\Psi = \{y\in\mathbb{C}: {\mathfrak{R}\{y\}}=0, {\mathfrak{I}\{y\}}\geq 0\}$. A block diagram is depicted in [Fig.]{}. The information embedded in a single eigenvalue $\lambda\in\Lambda$ is transformed to a time-domain signal $q(t,0)$ via the where the transmitter is located at position $z=0$ along the fiber. At position $z=1$, the receiver calculates the discrete spectrum $\psi\in\Psi$ from the received signal $q(t,1)$ via the . The time-domain signal corresponds to first order solitons, i.e., $$\begin{aligned}
    q(t,0) &= 2{\mathfrak{I}\{\lambda\}} {\operatorname{sech}{\mathopen{}\mathclose\bgroup\originalleft}(2{\mathfrak{I}\{\lambda\}} t{\aftergroup\egroup\originalright})}.\end{aligned}$$ For the to be valid, the signal must have finite support, i.e., before transmitting the next pulse, the previous one must have returned to zero. As the pulses in general have infinite tails, we truncate them when they fall below a threshold close to zero. We define the pulse over the smallest support containing a fraction $(1-\delta ) $ of the energy of the pulse and hence, we can formally define the pulse width as follows.

The pulse width of $\lambda$ is defined as the smallest support containing a fraction $(1-\delta ) $ of the energy of the pulse, $$\begin{aligned}
    T(\lambda, \delta ) &\triangleq \frac{1}{2{\mathfrak{I}\{\lambda\}}} \ln{\mathopen{}\mathclose\bgroup\originalleft}(\frac{2}{\delta} -1{\aftergroup\egroup\originalright}),\end{aligned}$$ where $0 < \delta < 1$.

The value of the cutoff parameter $\delta$ must be chosen in a way such that soliton-soliton interactions are negligible. For longer transmission distances, $\delta$ decreases, i.e., the pulses must be spaced further apart. Furthermore, the condition $$\begin{aligned}
    {\operatorname{e}^{-2{\mathfrak{I}\{\lambda\}} \Delta(\lambda, \delta )}} 
    &= {\operatorname{e}^{-\ln{\mathopen{}\mathclose\bgroup\originalleft}(\frac{2}{\delta} -1{\aftergroup\egroup\originalright})}} \ll 1
    \label{eqn:guard_time_condition}\end{aligned}$$ must be fulfilled [@shevchenko2018tcom].

At this point, it is important to comment on the memorylessness of the system emanating from the absence of soliton-soliton interactions. A pulse train of well-separated first order solitons was investigated in [@shevchenko2018tcom] for launch powers of $\SI{-1.5}{\dBm}$ and $\SI{1.45}{\dBm}$ and transmission over $\SI{500}{\kilo\meter}$ and $\SI{2000}{\kilo\meter}$. It was shown via simulations that the correlation between the symbols at the receiver is essentially zero, concluding that the channel is indeed memoryless in the transmission range of $\SI{500}{\kilo\meter}$ to $\SI{2000}{\kilo\meter}$ and transmit power range of $\SI{-1.5}{\dBm}$ to $\SI{1.45}{\dBm}$ for which the model [()]{} is applicable. While this approach is not a rigorous proof, the results indicate that memorylessness is a valid assumption. Although the transmission scheme is different in [@shevchenko2018tcom], the underlying condition that any two pulses need to be sufficiently separated is the same. Hence, we can treat the -based transmission system in this work as a memoryless channel.

In a practical system, we assume distributed Raman amplification and noise with received power spectral density $\sigma^2$ to compensate for the lossy fiber and be able to use the to relate the input and the output. The conditional of such a system has been derived via a perturbative approach and the Fokker-Planck equation method [@derevyanko2005] and is used to design a communication system in [@shevchenko2018tcom]. It is given by

$$\begin{aligned}
    p_{Y|X}(\psi|\lambda) &= \frac{2}{\sigma^2}\sqrt{\frac{{\mathfrak{I}\{\psi\}}}{{\mathfrak{I}\{\lambda\}}}}{\operatorname{e}^{-2\frac{{\mathfrak{I}\{\lambda\}}+{\mathfrak{I}\{\psi\}}}{\sigma^2}}}{\operatorname{I}_{1}{\mathopen{}\mathclose\bgroup\originalleft}(\frac{4\sqrt{{\mathfrak{I}\{\lambda\}}{\mathfrak{I}\{\psi\}}}}{\sigma^2}{\aftergroup\egroup\originalright})},
\label{eqn:pdf_yx_sac}\end{aligned}$$

where $\psi$ is the received symbol as in [Fig.]{}, and ${\operatorname{I}_{1}{\mathopen{}\mathclose\bgroup\originalleft}(\cdot{\aftergroup\egroup\originalright})}$ is the modified Bessel function of the first kind of order one. The power spectral density of the received noise $\sigma^2$ is normalized and relates to real world units as $\sigma^2 = \gamma\sqrt{L^3}\sigma_0^2 / {\mathopen{}\mathclose\bgroup\originalleft}(\sqrt{2|\beta_2|}{\aftergroup\egroup\originalright})$. The is defined as $\mathrm{SNR}\triangleq 4{\mathbb{E}_{X}\!{\mathopen{}\mathclose\bgroup\originalleft}\{{{\mathfrak{I}\{\lambda\}}}{\aftergroup\egroup\originalright}\}
}/\sigma^2$. It is important to note that the model [()]{} assumes the noise intensity to be small such that it can be treated as a perturbation to the soliton. Hence, the model is only applicable if the signal energy is not the same order as that of the noise. Furthermore, for very high signal powers, [()]{} is no longer valid either since the impact of the inelastic scattering effects (i.e., stimulated Raman or Brillouin scattering) is not considered within the 1st-order perturbation approach. For a detailed derivation of the model, we refer the reader to [@derevyanko2005].

In [@shevchenko2018tcom], the shortest possible symbol interval is defined by the pulse duration of $\lambda_1$, i.e., the longest pulse. However, this tends to be inefficient since especially for short pulses, the guard interval between two consecutive pulses is longer than necessary and thereby limits the data rate. Here, we exploit the effect of varying pulse lengths and transmit each pulse as soon as the previous one has returned to zero. This concept is depicted in [Fig.]{}, where pulse sequences with fixed and varying symbol interval are compared. The figure clearly shows the advantage of a varying pulse interval and also demonstrates the aforementioned inefficiencies. The data rate of a system with varying symbol intervals depends on the distribution of the data. Thus, we define the average symbol interval as follows.

![Comparison of a pulse sequence with static symbol intervals and dynamic symbol intervals.[]{data-label="fig:system:plot_comparison_spacing"}](plot_comparison_spacing)

\[def:avg\_duration\] The average symbol interval is $$\begin{aligned}
\bar{T}(X) &\triangleq \sum_{k=1}^{M} p_X (\lambda_k) {T{\mathopen{}\mathclose\bgroup\originalleft}(\lambda_k{\aftergroup\egroup\originalright})} = {\mathbb{E}_{X}\!{\mathopen{}\mathclose\bgroup\originalleft}\{{{T{\mathopen{}\mathclose\bgroup\originalleft}(\lambda{\aftergroup\egroup\originalright})}}{\aftergroup\egroup\originalright}\}
}.\end{aligned}$$

In [@shevchenko2018tcom], only eigenvalues with an imaginary part larger than zero are used. We extend this by allowing ${\mathfrak{I}\{\lambda\}}=0$. In the time domain, this results in a pulse with amplitude zero, i.e., we do not transmit anything. We define its corresponding duration as the same as the duration of the shortest pulse, ${T{\mathopen{}\mathclose\bgroup\originalleft}(\lambda=0{\aftergroup\egroup\originalright})} \triangleq {T{\mathopen{}\mathclose\bgroup\originalleft}(\lambda_M{\aftergroup\egroup\originalright})}$.

As any practical system can handle only a maximum peak power and a maximum bandwidth, we enforce a peak power constraint which relates to a maximum eigenvalue constraint. Especially in systems with lumped amplification and , such a constraint is required as eigenvalues fluctuate depending on their amplitude, which decreases the performance [@zafrulla2003].

We note that the varying symbol interval introduces additional challenges on detection. In particular, an erroneously detected symbol may lead to error propagation, insertion errors (detection of symbols when none was transmitted), deletion errors (not detecting a transmitted symbol), or the loss of synchronization. To calculate the capacity, however, we neglect these effects. Hence, the results can be seen as an upper bound on the performance.

Capacity Achieving Distribution {#sec:optimum_input_distribution}
===============================

From [Fig.]{}, it is intuitive that pulses with short duration should be transmitted more frequently than pulses with long duration. However, shorter pulses are more perturbed by noise than longer pulses. Hence, the optimal input distribution to the channel as described by the conditional [()]{} is not the conventional uniform distribution. The channel capacity is obtained by maximizing the , $$\begin{aligned}
    {\mathbb{I}(X;Y)} &\triangleq {\mathbb{E}_{X,Y}\!{\mathopen{}\mathclose\bgroup\originalleft}\{{\log_2{\mathopen{}\mathclose\bgroup\originalleft}({\frac{p_{Y|X}(Y|X)}{\sum_{\tilde{\lambda}\in\Lambda}p_{Y|X}(Y|\tilde{\lambda}) p_X(\tilde{\lambda})}}{\aftergroup\egroup\originalright})}{\aftergroup\egroup\originalright}\}
}\end{aligned}$$ over all possible input distributions $p_X(\lambda )$. Here, due to the variable transmission duration, we need to consider the under a variable cost constraint $\bar{T}(\cdot)$ [@verdu1990], $$\begin{aligned}
    {\mathsf{I}(X;Y)} &\triangleq \frac{{\mathbb{I}(X;Y)}}{\bar{T}(X)}.
    \label{eqn:mutual_information_time_scaled}\end{aligned}$$ To emphasize that the cost of a symbol is its corresponding pulse duration, we refer to the in the form of [()]{} as time-scaled . We can therefore define the capacity as $$\begin{aligned}
{\mathsf{C}}&\triangleq \sup_{p_X(\lambda)} {\mathsf{I}(X;Y)}
\label{eqn:def_cstar}\end{aligned}$$ where we set the supremum to zero if the set of distributions therein is empty. The capacity-achieving distribution, denoted by $p_X^*(\lambda )$, is in the set for which the supremum is non-zero.

As the ${\mathbb{I}(X;Y)}$ is concave in $p_X(\lambda )$ and $\bar{T}(X)$ is linear in $p_X(\lambda )$ and positive, the time-scaled ${\mathsf{I}(X;Y)}$ is quasiconcave [@stancu-minasian_frac_programming Table 2.5.2]. We can solve [()]{} and obtain the corresponding capacity-achieving distribution numerically.

![Optimal distribution for different .[]{data-label="fig:optDist:example_opt_dist"}](optimal_distribution_combined_sufigures "fig:") -0.5cm

Exemplary results of the capacity-achieving distribution are shown in [Fig.]{}. We note that the lowest and highest amplitudes are always used with equal and high probability. For low , only these are used, i.e., is optimal. Furthermore, the capacity-achieving distribution is discrete and is of exponential-like shape with the exception of a point mass at zero as it can be seen in [Fig.]{}.

Note that ${\mathsf{C}}$ assumes memorylessness, which does not necessarily hold due to the variable symbol interval. Hence, ${\mathsf{C}}$ is, in fact, the constraint capacity under the assumption of a memoryless channel and the constraint of transmitting only first-order solitons. However, for notational simplicity, we refer to it simply as capacity with its corresponding capacity-achieving distribution.

In the case of a noiseless channel, it is possible to derive a closed form solution to [()]{} under the assumption of a finite discretization.

\[lemma:noiseless\] Let $\lambda_1, \lambda_2, \ldots, \lambda_M$ be $M \geq 2$ eigenvalues with $0 \leq {\mathfrak{I}\{\lambda_1\}} < {\mathfrak{I}\{\lambda_2\}} < \ldots < {\mathfrak{I}\{\lambda_M\}}$ and let ${T{\mathopen{}\mathclose\bgroup\originalleft}(\lambda_k{\aftergroup\egroup\originalright})}$ be the time of transmitting a pulse with eigenvalue $\lambda_k$. Let $r$ be the unique real positive root of the polynomial $\sum_{k=1}^Mx^{-T(\lambda_k)} - 1$. Then, in the noiseless case, the capacity is obtained as $${\mathsf{C}}= \log_2(r)$$ and the capacity-achieving distribution is given by $$\begin{aligned}
 P_X^\diamond(\lambda_k) &= {\operatorname{e}^{-\ln(r){T{\mathopen{}\mathclose\bgroup\originalleft}(\lambda_k{\aftergroup\egroup\originalright})}}}, & k&=1,\ldots,M.
\label{eqn:opt_dist_noiseless}\end{aligned}$$

Suppose that the $k$-th eigenvalue is transmitted with probability $P_k$. For any fixed average symbol interval $\bar{T}(X) = \sum_kP_kT(\lambda_k)$, where ${T{\mathopen{}\mathclose\bgroup\originalleft}(\lambda_M{\aftergroup\egroup\originalright})} \leq \bar{T}(X) \leq {T{\mathopen{}\mathclose\bgroup\originalleft}(\lambda_1{\aftergroup\egroup\originalright})}$, we are interested in the distribution that maximizes the entropy while leading to the average symbol duration $\bar{T}(X)$. It is known that this distribution takes the form [@cover_inf_theory Ch. 12] $$\begin{aligned}
    P_k &= \frac{{\operatorname{e}^{-\theta {T{\mathopen{}\mathclose\bgroup\originalleft}(\lambda_k{\aftergroup\egroup\originalright})}}}}{\xi(\theta)}
    \label{eqn:maximizer}\end{aligned}$$ where $\xi(\theta) = \sum_i{\operatorname{e}^{-\theta T(\lambda_i)}}$ ensures that $\sum_k P_k = 1$ and $\theta$ has to be selected such that $\sum_kP_kT(\lambda_k) = \bar{T}$. In the noiseless case, the MI is given by ${\mathbb{I}(X;Y)} = {\mathbb{H}(X)}$. The entropy ${\mathbb{H}(X)}$ then is $$\begin{aligned}
    {\mathbb{H}(X)} =: {\mathbb{H}(\theta)} &= -\sum_{k=1}^M P_k \log_2{\mathopen{}\mathclose\bgroup\originalleft}(\frac{{\operatorname{e}^{-\theta {T{\mathopen{}\mathclose\bgroup\originalleft}(\lambda_k{\aftergroup\egroup\originalright})}}}}{\xi(\theta)}{\aftergroup\egroup\originalright})\\
                                     &= \frac{1}{{\ln(2)}}\sum_{k=1}^M P_k (\theta {T{\mathopen{}\mathclose\bgroup\originalleft}(\lambda_k{\aftergroup\egroup\originalright})} + \ln(\xi(\theta)))\\
                                     &= \frac{\theta \bar{T}(X)}{{\ln(2)}} + \log_2(\xi(\theta)).\end{aligned}$$ The time-scaled MI hence takes the form $$\begin{aligned}
 {\mathsf{I}(X;Y)} &=  \frac{\theta}{\ln(2)} + \frac{\log_2(\xi(\theta))}{\bar{T}(X)} \\
 &= \frac{\theta}{\ln(2)} + \frac{\log_2(\xi(\theta))}{\sum_kP_kT(\lambda_k)} \\
 &= \frac{\theta}{\ln(2)} + \frac{\log_2(\sum_k{\operatorname{e}^{-\theta T(\lambda_k)}})\sum_k{\operatorname{e}^{-\theta T(\lambda_k)}}}{\sum_k{\operatorname{e}^{-\theta T(\lambda_k)}}T(\lambda_k)}\,. \end{aligned}$$ In order to maximize ${\mathsf{I}(X;Y)}$, we find the optimal parameter $\theta$ by setting $\xi(\theta) = 1$. This can be seen by setting the derivative of ${\mathsf{I}(X;Y)}$ to zero, with $$\begin{gathered}
\frac{\partial}{\partial\theta} {\mathsf{I}(X;Y)} = \\
 \log_2{\mathopen{}\mathclose\bgroup\originalleft}(\sum_k{\operatorname{e}^{-\theta T(\lambda_k)}}{\aftergroup\egroup\originalright}){\mathopen{}\mathclose\bgroup\originalleft}(\frac{\sum_k{\operatorname{e}^{-\theta T(\lambda_k)}}}{\sum_kT(\lambda_k){\operatorname{e}^{-\theta T(\lambda_k)}}}{\aftergroup\egroup\originalright})^2\mathop{\mathrm{var}}(T(\lambda))\, ,\end{gathered}$$ where $\mathop{\mathrm{var}}(T(\lambda))$ denotes the variance of $T(\lambda_k)$ for the given $\theta$. By assumption, as all $T(\lambda_k)$ are different, the middle part of this expression is strictly positive and $\mathop{\mathrm{var}}(T(\lambda)) > 0$. Hence, it is easy to see that this derivative can only be zero if $\sum_k{\operatorname{e}^{-\theta T(\lambda_k)}} = 1$. The optimal $\theta$ is hence found by setting $\xi(\hat{\theta}) = 1$. Consider the polynomial $$f(x) = \sum_{k=1}^Mx^{-T(\lambda_k)} - 1.$$ As this polynomial is monotonically decreasing for positive $x$, with $\lim_{x\to 0^+}f(x) = +\infty$ and $\lim_{x\to+\infty}f(x)=-1$, $f(x)$ has exactly one positive real root. Let $r$ be the unique positive real root of $f(x)$. Then $\hat{\theta} = \ln(r)$. Inserting $\hat{\theta}$ into ${\mathsf{I}(X;Y)}$ and  proves the lemma.

![Time-scaled of the optimal distribution for linearly spaced constellations with $M$ points (colored with markers solid), and of a system as in [@shevchenko2018tcom] (with markers dotted). As a reference, the capacity ${\mathsf{C}}$ is plotted as well (black solid without markers). For the cutoff parameter, $\delta=0.005$ was used.[]{data-label="fig:shaping:mi_comparison"}](mi_comparison)

We clearly see that [()]{} is of exponential shape with an additional point mass at zero. Furthermore, we note that the shape of the distribution is mostly caused by the variable pulse duration. The noise then determines the optimal location and optimal number of constellation points.

For a transmission system, the is an upper bound on the achievable rate. In [Fig.]{} we evaluate the time-scaled for various input distributions for a cutoff parameter $\delta=0.005$. The capacity is depicted with a black solid line. To reduce the complexity of implementation, we constrain the constellation $\Lambda$ to $M$ linearly spaced points from $\lambda_1=0$ to $\lambda_M$, i.e., $$\begin{aligned}
    \lambda_i &= (i-1)\frac{\lambda_M}{M-1} \text{ for } i=1,\ldots,M,\end{aligned}$$ and plot the corresponding time-scaled in colored solid lines with markers. We note that the time-scaled is very close to the capacity curve until it saturates. Increasing the modulation order $M$ shows significant increase in the time-scaled . For comparison purposes, we also plot the time-scaled for a system with fixed symbol duration and conventional uniform distribution on a linearly spaced constellation as in [@shevchenko2018tcom]. We observe that the rate saturates at very low values and that increasing the modulation order $M$ shows only slight improvement.

![Block diagram of the scheme.[]{data-label="fig:shaping:system_block_diagram_pes"}](block_diagram_pes)

![image](parity_constellation_combined_subfigure)

Probabilistic Eigenvalue Shaping {#sec:shaping_time_share}
================================

In the previous section, we observed a significant gap between the time-scaled of the system in [@shevchenko2018tcom] and the capacity. This gap is referred to as shaping gap. In order to close it, we propose a system as shown in [Fig.]{}, inspired by  [@bocherer2015_bw_efficient_rate_matched_ldpc].

In the scheme, the sequence of uniformly distributed data bits is mapped to a sequence of positive amplitudes distributed half Gaussian by a . The binary image of this sequence is encoded by a systematic code, resulting in uniformly distributed parity bits, which are then used to map the sequence of half Gaussian distributed symbols to a stream of Gaussian distributed symbols.

As the capacity-achieving distribution $p_X^*(\lambda)$ is not symmetric, cannot be directly applied here. However, in order to keep the benefits of , we wish to apply the before the . We describe in the following with reference to [Fig.]{}. The binary data sequence $\bm{u}$ of length $k_\mathsf{s}$ bits is mapped by the to a sequence of eigenvalues $\bm{\lambda}\in\Lambda^{n_\mathsf{s}}$ of length $n_\mathsf{s}$ distributed according to $p_X^*(\lambda)$. The can be used for that purpose [@schulte2016_ccdm]. It is asymptotically optimal as its rate $R_\mathsf{s}$ approaches the entropy of the desired channel input $X$, $$\begin{aligned}
    R_\mathsf{s} &=\frac{k_\mathsf{s}}{n_\mathsf{s}} \rightarrow {\mathbb{H}(X)} \text{  as  } n_\mathsf{s} \rightarrow\infty.\end{aligned}$$ For large block sizes, the gap between $R_\mathsf{s}$ and ${\mathbb{H}(X)}$ is sufficiently small and can be neglected. Note that some of the possible eigenvalues may occur with probability zero.

We consider the modulation order $M$ to be a power of two such that we can define its binary image. The binary image of $\bm{\lambda}$, $\bm{bi}(\bm{\lambda} )$, is then encoded by a systematic encoder with information block length $k_\mathsf{c}$, code length $n_\mathsf{c}$, and rate $R_\mathsf{c}=\frac{k_\mathsf{c}}{n_\mathsf{c}}$. The code is denoted by $\mathcal{C}$, with $|\mathcal{C}|=2^{k_\mathsf{c}}$. The parity bits at the output of the encoder are mapped to a sequence of eigenvalues $\bm{\lambda}_\mathsf{par}\in \Lambda_\mathsf{par}$ with modulation order $M_\mathsf{par}=|\Lambda_\mathsf{par}|$ and $\Lambda_\mathsf{par}\subseteq\Lambda$ by the block $\bm{s}(\cdot )$ in [Fig.]{} such that they are uniformly distributed.

Assuming that a high code rate $R_\mathsf{c}$ is used, we accept a small penalty with respect to the optimal channel input distribution and transmit $\bm{\lambda}$ and $\bm{\lambda}_\mathsf{par}$ via time-sharing. The major difference of compared to is the fact that the channel input distribution is not the optimal distribution due to the time-sharing with the sequence $\bm{\lambda}_\mathsf{par}$. Consequently, this causes a performance degradation. However, is highly flexible as the spectral efficiency can be adapted by the and the code rate $R_\mathsf{c}$, and a single code can be used. Note that every eigenvalue is protected by the code as is performed after the and decoding and demapping can be performed independently. Thus, shares these advantages with .

We wish for a high code rate $R_\mathsf{c}$ to keep the performance degradation due to the time-sharing low. More precisely, we wish to maximize the number of symbols distributed according to $p_X^*(\lambda)$. The ratio between information symbols and coded symbols, denoted by $R_\mathsf{ts}$, is an indication for the expected performance degradation, $$\begin{aligned}
R_\mathsf{ts} &= \frac{\frac{n_\mathsf{c} R_\mathsf{c}}{\log_2(M)}}{\frac{n_\mathsf{c}R_\mathsf{c}}{\log_2(M)} + \frac{n_\mathsf{c}(1-R_\mathsf{c})}{\log_2(M_\mathsf{par})}}\nonumber\\
&= \frac{R_\mathsf{c}\log_2(M_\mathsf{par})}{\log_2(M)(1-R_\mathsf{c})+R_\mathsf{c}\log_2(M_\mathsf{par})}.
\label{eqn:shaping:R_ts}\end{aligned}$$

Parity symbols
--------------

The parity symbols at the output of the code encoder are uniformly distributed. In [Fig.]{}, we observed that with uniform signaling, i.e., $\Lambda_\mathsf{par}=\{\lambda_1, \lambda_M\}$ and $M_\mathsf{par}=2$, is optimal for low as it achieves capacity and performs reasonably well for high . However, we note from [Fig.]{} that for a higher order modulation, even with uniform signaling, higher rates are possible. Hence, here we consider a scenario where $M_\mathsf{par}>2$. We further increase the rate by only using a subset of $\Lambda$ and by picking the eigenvalues such that they are not uniformly spaced.

Consider the information symbol alphabet $\Lambda=\{\lambda_1,\ldots,\lambda_8\}$ with $M=8$. For the $\Lambda_\mathsf{par}$, we could pick $\Lambda_\mathsf{par}=\{\lambda_1, \lambda_6, \lambda_7,\lambda_8\}$ with $p_X(\lambda)=\{0.25, 0.25, 0.25, 0.25\}$ and $M_\mathsf{par}=4$.

To find the function $\bm{s}(\cdot )$ that maps the parity symbols onto $\lambda\in\Lambda_\mathsf{par}$, we use a greedy algorithm as described in [Algorithm]{}. It starts with , i.e., $\Lambda_\mathsf{par}=\{\lambda_1, \lambda_M\}$. For each of the remaining symbols $\lambda\in\Lambda\setminus\Lambda_\mathsf{par}$, it calculates the time-scaled of $\lambda\cup\Lambda_\mathsf{par}$, finds the symbol $\lambda$ for which the time-scaled of $\lambda\in\Lambda\setminus\Lambda_\mathsf{par}$ is maximized, and adds it to $\Lambda_\mathsf{par}$. All symbols with a greater or equal imaginary part than $\lambda$ are removed, i.e., the eigenvalues $\{\lambda'\in\Lambda : {\mathfrak{I}\{\lambda'\}}\geq{\mathfrak{I}\{\lambda\}}\}$ are removed. This process is repeated until there are no symbols left. We then choose the set of symbols that gives the highest time-scaled as $\Lambda_\mathsf{par}$. We note that this procedure does not guarantee an optimal solution. However, for $M=\{4,8\}$ an exhaustive search gives the same result as that of [Algorithm]{}.

Constellation $\Lambda$ Constellation $\Lambda_\mathsf{par}$ $\Lambda_\mathrm{placed} = \{\lambda_1, \lambda_M\}$ $\Lambda_\mathsf{par} = \{\lambda_1, \lambda_M\}$

$\Lambda_\mathrm{not\,placed} = \Lambda \setminus \Lambda_\mathsf{par}$ Calculate ${\mathsf{I}(\Lambda_\mathrm{placed} \cup \lambda_i)}$ $\lambda_\mathrm{max} := \arg\max {\mathsf{I}(\cdot)}$ $\Lambda_\mathrm{placed} = \Lambda_\mathrm{placed} \cup \lambda_\mathrm{max}$ $\Lambda_\mathsf{par} = \Lambda_\mathrm{placed}$ $\Lambda_\mathrm{not\,placed} = \Lambda_\mathrm{not\,placed}\setminus\{\lambda:\lambda\in\Lambda_\mathrm{not\,placed}, {\mathfrak{I}\{\lambda\}}\geq{\mathfrak{I}\{\lambda_\mathrm{max}\}}\}$

$\Lambda_\mathsf{par}$

In [Fig.]{}, we show $\Lambda_\mathsf{par}$ for different modulation orders and . For $M=4$, we note that for low gives the best result. Increasing the results in a third level being added. The same behavior is observed for $M=8$. Compared to $M=4$, the third level is introduced at a slightly lower . This results from the fact that for $M=8$, different constellation points are available. For $M=16$, we note that again a third level appears when increasing the . When further increasing it, this third level moves to an eigenvalue with larger imaginary part and consequently a fourth level at an eigenvalue with lower imaginary part appears. This behavior can be observed repeatedly. To map the binary parity bits to the constellation points, we require $M_\mathsf{par}$ to be a power of two. As this is not always the case (see [Fig.]{}), we pick the largest power of two that is smaller or equal than the number of constellation points given by [Algorithm]{}.

Achievable Rate of Probabilistic Eigenvalue Shaping {#sec:achievable_rate}
---------------------------------------------------

To characterize the performance of , we derive the achievable rate of , denoted by $\mathsf{R}_\mathsf{ps}$. We assume that the channel is memoryless and that the decoder performs bit-metric decoding.

The achievable rate of is $$\begin{aligned}
\mathsf{R}_\mathsf{ps} &= R_\mathsf{ts}{\mathopen{}\mathclose\bgroup\originalleft}({\mathbb{H}(X)}  - \sum\limits_{i=1}^{m}{\mathbb{H}(X^{\mathsf{B}}_i|Y^{\mathsf{B}}_i)} {\aftergroup\egroup\originalright}) \nonumber\\ 
&\quad + {\mathopen{}\mathclose\bgroup\originalleft}(1-R_\mathsf{ts}{\aftergroup\egroup\originalright}){\mathopen{}\mathclose\bgroup\originalleft}(m_\mathsf{par} - \sum\limits_{i=1}^{m_\mathsf{par}}{\mathbb{H}(X^{\mathsf{B}}_{\mathsf{par},i}|Y^{\mathsf{B}}_{\mathsf{par},i})}{\aftergroup\egroup\originalright}).\label{eqn:air}\end{aligned}$$

The achievable rate for has been derived in [@bocherer2017_achievable_rate_ps]. For a system employing time-sharing, the resulting achievable rate is the average of the achievable rate of the two transmission schemes.

![Achievable rates for different code rates with $\Lambda_\mathsf{par}$ according to [Algorithm]{} for a cutoff parameter $\delta=0.005$.[]{data-label="fig:shaping:achievable_rate_air"}](air_time_share)

In [Fig.]{}, we plot the capacity and the achievable rate [()]{} for different code rates $R_\mathsf{c}=\{$$1/4$, $1/3$, $2/5$, $1/2$, $3/5$, $2/3$, $3/4$, $4/5$, $5/6$, $8/9$, $9/10\}$ and modulation orders for a cutoff parameter $\delta=0.005$. $\Lambda_\mathsf{par}$ and hence $M_\mathsf{par}$ are chosen according to the results of [Algorithm]{}. For each modulation order, we notice that the curves cross at a certain . For below this point, the lowest code rate (corresponding to the highest curve) gives the best performance whereas for above this point, the highest code rate (corresponding to the highest curve) gives the best performance. We note the influence of time-sharing, which results in a gap between the achievable rate and capacity. The gap increases for lower code rates $R_\mathsf{c}$ as the channel input distribution deviates more from the optimal one.

Numerical Evaluation {#sec:performance}
====================

In this section, we evaluate the performance of the scheme via discrete-time Monte-Carlo and simulations. For the mapping $\bm{bi}(\cdot)$ (see [Fig.]{}), we use Gray labeling. Also, for the , we use the binary codes of the standard with code length $n_\mathsf{c}=64800$ and code rates $R_\mathsf{c}=\{$$1/4$, $1/3$, $2/5$, $1/2$, $3/5$, $2/3$, $3/4$, $4/5$, $5/6$, $8/9$, $9/10\}$. For the parity symbols, we use the constellation arising from [Algorithm]{}, depicted in [Fig.]{}.

Detection {#sec:detection}
---------

For the simulation, we simulate a continuous signal and hence, we require a detector. We use the following method to deal with the variable pulse durations: We set a threshold $\theta$ sufficiently higher than the noise. Once the magnitude of the signal rises above $\theta$, we save the time as $t_\mathsf{start}$ and when the magnitude of the signal falls below $\theta$, we save the time as $t_\mathsf{end}$. We then extend the interval bounded by $t_\mathsf{start}$ and $t_\mathsf{end}$, i.e., $\tilde{t}_\mathsf{start} = t_\mathsf{start} - \delta_t$ and $\tilde{t}_\mathsf{end} = t_\mathsf{end} + \delta_t$. Calculating the over the interval $[\tilde{t}_\mathsf{start}, \tilde{t}_\mathsf{end}]$ using the spectral method [@yousefi2014nftI-III Part II, Section IV] and only considering the imaginary part of the discrete eigenvalue gives the received symbol $y$. This approach requires that the is sufficiently high. As the model has the same requirement due to the perturbation approach, this requirement is fulfilled.

It may happen that due to noise, a received pulse never rises above the threshold $\theta$. In this case, the shortest duration is assumed (i.e., the duration of the pulse with amplitude zero). This scenario can be avoided by choosing the threshold sufficiently lower than the lowest amplitude. Furthermore, due to the shape of the capacity-achieving distribution, lower amplitudes are less likely, hence preventing this scenario.

To find the best threshold, we tested the performance for different values of $\theta$ and found that the performance of a threshold at $75 \%$ of the lowest non-zero amplitude of the constellation works best. We observed that small deviations of the threshold do not affect the performance significantly whereas setting the threshold too high (missing symbols with low amplitude) or to low (detecting a symbol where there is none) leads to performance degradation. Furthermore, we assume synchronization sequences spread sufficiently far apart in order not to impact the rate. We assume synchronization to be ideal such that it is guaranteed that error propagation is limited.

Numerical Results
-----------------

We perform Monte-Carlo simulations of the discrete-time model [()]{} and show the results in [Fig.]{}, where we plot the transmission rate at a of $10^{-5}$ for $M=4,8$ and $16$. The highest transmission rate for each modulation order corresponds to the highest code rate $R_\mathsf{c}$. We notice that the gap to capacity for $M=4$ is smaller than for $M=8$ and $M=16$. If we consider $\Delta M = M - M_\mathsf{par}$, i.e., the difference of the modulation order of $\Lambda$ and $\Lambda_\mathsf{par}$, we note that for a low $M$, $\Delta M$ is low was well. For example, for $M=4$, $\Delta M \leq 2$. Hence, the rate loss due to time-sharing is small. For $M=16$, the gap to capacity is smaller than for $M=8$. Considering the relevant range, we note that $\Delta M$ is smaller for $M=16$ than for $M=8$ and thus explaining the smaller rate loss.

We also simulated the transmission over a fiber using simulations transmitting a train of solitons. We consider a with parameters as in [Table]{} and two different amplification schemes, distributed Raman amplification and lumped amplification using . For both schemes, the peak power constraint is chosen such that the effect of the can be neglected, i.e., $\lambda_\mathsf{max}=2\jmath$. We employ the detection schemes as described in and choose the cutoff-parameter $\delta=0.005$, i.e., $99.5\,\%$ of the energy is contained in the pulse, for which the condition [()]{} is fulfilled. This then leads to a similar cutoff parameter as in [@shevchenko2018tcom]. For each modulation order $M=\{4,8,16\}$, we determine the furthest distance over which we achieve a of less than $10^{-5}$ and consider the rate gain compared to an unshaped system as in [@shevchenko2018tcom]. This results for $M=\{4,8,16\}$ in transmission over $\SI{3200}{\kilo\meter}$, $\SI{3040}{\kilo\meter}$, and $\SI{2960}{\kilo\meter}$ at a rate gain of $\SI{20}{\percent}$, $\SI{26}{\percent}$, and $\SI{95}{\percent}$, respectively. The results do not differ for distributed and lumped amplification as this is ensured by the peak power constraint.

Conclusion {#sec:conclusion}
==========

In this paper, we presented a probabilistic shaping scheme for an -based transmission system embedding information in the imaginary part of the discrete spectrum. It shapes the information symbols according to the capacity-achieving distribution and transmits them via time-sharing together with the uniformly distributed, suitably modulated parity symbols. We exploited the fact that the pulses of the signal in the time domain are of unequal length to improve the data rate compared to [@shevchenko2018tcom]. We used the time-scaled and derived the capacity-achieving distribution in closed form for the noiseless case and numerically in the general case. We showed that significantly improves the performance of an -based transmission scheme, and can almost double the data rate. As a possible extension of our work, the continuous spectrum can be used to increase the spectral efficiency [@aref2016ecoc].

![Performance of time sharing with parity symbols according to [Algorithm]{}. The rate points correspond to a performance at $\text{BER}=10^{-5}$. The highest transmission rate for each modulation order corresponds to the highest code rate $R_\mathsf{c}$.[]{data-label="fig:performance:mi_monte_carlo_perf_geometric"}](mi_monte_carlo_perf)

Acknowledgments {#acknowledgments .unnumbered}
===============

The authors would like to thank the anonymous reviewers for their feedback and comments which helped to improve this paper significantly. Especially, we would like to acknowledge one of the reviewers for proposing an elegant way to prove [Lemma ]{}, which is included in this paper.

[10]{} \[1\][\#1]{} url@samestyle \[2\][\#2]{} \[2\][[l@\#1=l@\#1\#2]{}]{}

M. I. Yousefi and F. R. Kschischang, “Information transmission using the nonlinear fourier transform, part [I]{}-[III]{},” *[IEEE]{} Trans. Inf. Theory*, vol. 60, no. 7, pp. 4312–4369, Jul. 2014.

Z. Dong, S. Hari, T. Gui, K. Zhong, M. I. Yousefi, C. Lu, P. K. A. Wai, F. R. Kschischang, and A. P. T. Lau, “Nonlinear frequency division multiplexed transmissions based on [NFT]{},” *[IEEE]{} Photon. Technol. Lett.*, vol. 27, no. 15, pp. 1621–1623, Aug. 2015.

V. Aref, H. B[ü]{}low, K. Schuh, and W. Idler, “Experimental demonstration of nonlinear frequency division multiplexed transmission,” in *Proc. 41st Eur. Conf. Opt. Commun. ([ECOC]{})*, Valencia, Spain, Sep. 2015, pp. 1–3.

V. Aref, S. T. Le, and H. B[ü]{}low, “Demonstration of fully nonlinear spectrum modulated system in the highly nonlinear optical transmission regime,” in *Proc. 42nd Eur. Conf. Opt. Commun. ([ECOC]{})*, D[ü]{}sseldorf, Germany, Sep. 2016, pp. 1–3.

A. Geisler and C. Schaeffer, “Experimental nonlinear frequency division multiplexed transmission using eigenvalues with symmetric real part,” in *Proc. 42nd Eur. Conf. Opt. Commun. ([ECOC]{})*, D[ü]{}sseldorf, Germany, Sep. 2016, pp. 1–3.

S. Hari, M. I. Yousefi, and F. R. Kschischang, “Multieigenvalue communication,” *J. Lightw. Technol.*, vol. 34, no. 13, pp. 3110–3117, Jul. 2016.

N. A. Shevchenko, S. A. Derevyanko, J. E. Prilepsky, A. Alvarado, P. Bayvel, and S. K. Turitsyn, “Capacity lower bounds of the noncentral chi-channel with applications to soliton amplitude modulation,” *[IEEE]{} Trans. Commun.*, to appear.

G. D. Forney, R. Gallager, G. Lang, F. Longstaff, and S. Qureshi, “Efficient modulation for band-limited channels,” *[IEEE]{} J. Sel. Areas Commun.*, vol. 2, no. 5, pp. 632–647, Sep. 1984.

F.-W. Sun and H. C. A. van Tilborg, “Approaching capacity by equiprobable signaling on the [G]{}aussian channel,” *[IEEE]{} Trans. Inf. Theory*, vol. 39, no. 5, pp. 1714–1716, Sep. 1993.

G. B[ö]{}cherer, F. Steiner, and P. Schulte, “Bandwidth efficient and rate-matched low-density parity-check coded modulation,” *[IEEE]{} Trans. Commun.*, vol. 63, no. 12, pp. 4651–4665, Dec. 2015.

S. A. Derevyanko, S. K. Turitsyn, and D. A. Yakushev, “Fokker-planck equation approach to the description of soliton statistics in optical fiber transmission systems,” *J. Opt. Soc. Am. B*, vol. 22, no. 4, pp. 743–752, Apr. 2005.

M. Zafruullah, M. Waris, and M. K. Islam, “Simulation and design of [EDFA]{}s for long-haul soliton based communication systems,” in *Proc. Asia-Pacific Conf. Commun. ([APCC]{})*, Penang, Malaysia, Sep. 2003.

S. Verd[ú]{}, “On channel capacity per unit cost,” *[IEEE]{} Trans. Inf. Theory*, vol. 36, no. 5, pp. 1019–1030, Sep. 1990.

I. M. Stancu-Minasian, *Fractional Programming*, 1st ed.1em plus 0.5em minus 0.4emDordrecht, The Netherlands: Kluwer Academic Publishers, 1997.

T. M. Cover and J. A. Thomas, *Elements of Information Theory*, 2nd ed.1em plus 0.5em minus 0.4emHoboken, NJ, USA: Wiley, 2006.

P. Schulte and G. B[ö]{}cherer, “Constant composition distribution matching,” *[IEEE]{} Trans. Inf. Theory*, vol. 62, no. 1, pp. 430–434, Jan. 2016.

G. B[ö]{}cherer, “Achievable rates for probabilistic shaping,” *ArXiv e-prints*, Jul. 2017. \[Online\]. Available: <https://arxiv.org/abs/1707.01134>

[^1]: This work was funded by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 676448.

[^2]: A. Buchberger is with the Department of Electrical Engineering, Chalmers University of Technology, Gothenburg, SE-412 96, Sweden and Nokia Bell Labs, Lorenzstr. 10, 70435 Stuttgart, Germany, e-mail: andreas.buchberger@chalmers.se.

[^3]: A. Graell i Amat is with the Department of Electrical Engineering, Chalmers University of Technology, Gothenburg, SE-412 96, Sweden, e-mail: alexandre.graell@chalmers.se.

[^4]: V. Aref and L. Schmalen are with Nokia Bell Labs, Lorenzstr. 10, 70435 Stuttgart, Germany, e-mail: {firstname.lastname}@nokia-bell-labs.com.
