# Joint Power and Gain Allocation in MDM-WDM Optical Communication Networks
Based on Enhanced Gaussian Noise Model

M. A. Amirabadi  School of Electrical Engineering, Iran University of Science
and Technology (IUST), Tehran, 1684613114 Iran  M. H. Kahaei  School of
Electrical Engineering, Iran University of Science and Technology (IUST),
Tehran, 1684613114 Iran  Corresponding author<EMAIL_ADDRESS>S. A.
Nezamalhosseini  School of Electrical Engineering, Iran University of Science
and Technology (IUST), Tehran, 1684613114 Iran

###### Abstract

Achieving reliable communication over different channels and modes is one of
the main goals of Mode Division Multiplexing-Wavelength Division Multiplexing
(MDM-WDM) communication networks. The reliability can be described by minimum
Signal to Noise Ratio (SNR) margin which dependents on launched power,
Multimode-Erbium Doped Fiber Amplification (MM-EDFA) gain, and MMF
nonlinearity. In this paper, an analytical model for MMF nonlinearity is
derived based on Enhanced Gaussian Noise (EGN) model formulation by
considering carrier phase estimation and the first four dispersion terms. The
proposed EGN model is verified through the split step Fourier method
simulation. Considering a multi-node linear network, the joint optimized power
and gain allocation based on minimum SNR margin maximization is formulated.
The practical constraints including MM-EDFA saturation power and maximum gain
are considered and the problem is solved by using convex optimization. Three
scenarios are assumed including best equal power, optimized power, and joint
optimized power and gain. In the first scenario, equal powers are considered
for different channels and modes with equal MM-EDFA gain in all spans. It is
worth mentioning that the MM-EDFA gain is equal to span loss. In the second
scenario, different powers are allocated to different channels and modes with
equal MM-EDFA gain in all spans. In the third case, allocated powers to each
channel and mode are optimized. Moreover, the MM-EDFA gain for each span is
optimized separately. In MDM-single channel systems, simulation results
demonstrate that joint power and gain optimization leads to
$0.61\leavevmode\nobreak\ dB$ and $0.62\leavevmode\nobreak\ dB$ minimum SNR
margin improvement compared to optimal power and best equal power allocation,
respectively. In SMF-WDM systems considering the joint optimal power and gain
allocation, the minimum SNR margin improvements is about
$1.72\leavevmode\nobreak\ dB$ and $1.71\leavevmode\nobreak\ dB$ in comparison
with optimal power and best equal power allocation, respectively.

††journal: ao

##  1 Introduction

Mode Division Multiplexing (MDM) over Multimode Fibers (MMF) or multicore
fibers has emerged as a possible solution for overcoming the data-rate crunch
in optical communication networks. Theoretically, deploying an MMF propagating
$D$ spatial modes would increase capacity $D$ times [1]. The combination of
MDM with Wavelength Division Multiplexing (WDM) and polarization division
multiplexing schemes further increases the data rate.

The MDM-WDM systems suffer from both MMF linear and nonlinear effects. The MMF
linear effects include attenuation, chromatic and modal dispersion, and linear
coupling. The MMF nonlinear impairments include Kerr-effect-based nonlinearity
and nonlinear coupling [1]. The chromatic and modal dispersion, as well as the
linear coupling, can be mitigated by Multi-Input Multi-Output (MIMO) Digital
Signal Processing (DSP) techniques [2]. However, the MMF nonlinear effects are
the main limitations towards the practical implementation of MDM-WDM systems.

During the last decade, several theoretical works have considered MMF
nonlinearity analysis which focused mainly on numerical simulations [3] and
analytical predictions [4] combined with experimental verifications [5]. The
first step towards analyzing the MMF nonlinearity is solving the Manakov
equation [3]. The Split Step Fourier Method (SSFM) solves the Manakov equation
through many successive numerical simulation steps. The SSFM has high
computational complexity. The perturbation-based methods are the mostly used
analytical models that approximately solve the Manakov equation [6],[7], and
result in analytical formulations for predicting MMF nonlinearity. The
Gaussian Noise (GN) model is the most practical perturbation-based model which
describes the nonlinear effects by an additive Gaussian noise source [8], [9].
The GN model is accurate only for Gaussian shaped constellations in long-range
multi-span links.

In the first part of the paper, the Enhanced GN (EGN) model is derived by
considering Carrier Phase Estimation (CPE) and the first four dispersion terms
in order to improve the model accuracy. Moreover, the accuracy of the proposed
model is validated by SSFM. The proposed EGN formulation provides a clear
relationship between MMF nonlinearity with launched power of different
channels and modes as well as the MM-Erbium Doped Fiber Amplifier (MM-EDFA)
gain. Furthermore, the MM-EDFA produces the ASE noise which can be modeled by
an additive white Gaussian noise. Therefore, the joint power and gain
allocation problem can be composed by considering the ASE noise and nonlinear
noise formulations.

Achieving reliable communication is one of the main goals of MDM-WDM systems
which can be described by the minimum SNR margin. Therefore, the joint
optimized power and gain allocation considering minimum SNR margin
maximization is presented in the second part of the paper. Some practical
constraints including MM-EDFA saturation power and MM-EDFA maximum gain are
considered. The problem is solved considering a multi-node linear network by
using convex optimization approaches. The rest of this paper is organized as
follows. Section II presents the system and signal model, Section III provides
the EGN model formulation, Section IV presents the problem statement, Section
V brings the simulation results, and Section VI concludes the paper.

Figure 1: Schematic diagram of the MDM-WDM system.

##  2 System and signal model

###  2.1 System model

Fig. 1 shows the considered MDM-WDM system in which the input data is a
multiplexing of $N_{ch}$ channels, $D$ spatial modes, and $2$ polarization
modes. This link has $N_{s}$ spans with length $L_{s}$, combined by an MM-EDFA
at the end of each span for compensating the optical fiber loss. In order to
minimize the MMF nonlinear effect, both modal and chromatic dispersion are not
compensated [8]. The Kerr nonlinearity produced by inter/intra channel and
mode interactions is considered. It is assumed that the MMF imperfections
causing linear coupling among spatial modes be so weak, therefore, the linear
coupling effect is negligible [3]. Moreover, a MIMO DSP at the receiver is
considered for compensating MMF linear effects. Following, a CPE is used to
recover the average carrier phase rotation [7] due to the MMF nonlinear
effect.

###  2.2 Signal model

The following Ket notation represents the time domain of the optical launched
signal into the MMF link as

$\displaystyle\ket{A(0,t)}=\sum_{i_{1}=-\infty}^{\infty}\sum_{i_{2}=1}^{N_{ch}}\sum_{i_{3}=1}^{2D}\zeta_{\mathbf{i}}W_{Tx}^{i_{2},i_{3}}(t-i_{1}T_{i_{2}})e^{j2\pi
f_{i_{2}}t}\ket{i_{3}},$ (1)

where $\zeta_{\mathbf{i}}$ is the digital symbol (for example, QPSK) at time
index $i_{1}$, WDM channel index $i_{2}$, polarization-spatial mode index
$i_{3};i_{3}=1,…,2D$, and $\mathbf{i}=[i_{1},i_{2},i_{3}]$. Moreover,
$W_{Tx}^{i_{2},i_{3}}(t-i_{1}T_{i_{2}})$ is the transmitted pulse at time
index $i_{1}$. Moreover, $i_{2}$ and $i_{3}$ represent WDM channel and
polarization-spatial mode index, respectively. $T_{i_{2}}$ shows the symbol
duration, and $f_{i_{2}}$ is the carrier frequency. $\ket{i_{3}}$ represents a
one-hot vector wherein the $i_{3}$th element is one and the other elements are
zero. Actually, $i_{3}$ is used to denote the polarization-spatial mode of the
propagated signal. In other words, $i_{3}\triangleq p^{{}^{\prime}},p$, where
$p^{{}^{\prime}};p^{{}^{\prime}}=x,y$ is the polarization mode index and
$p;p=1,…,D$ represents the spatial mode index. It is obvious that $i_{3}$
takes values between $1$ and $2D$, since each spatial mode is a multiplexing
of $2$ polarization modes. Therefore, the time domain of the optical launched
signal into the MMF can be expressed as

$\displaystyle\ket{A(0,t)}=\sum_{i_{3}=1}^{2D}A_{i_{3}}(0,t)\ket{i_{3}}\triangleq\begin{bmatrix}A_{1}(0,t)\\\
A_{2}(0,t)\\\ …\\\ A_{2D-1}(0,t)\\\
A_{2D}(0,t)\end{bmatrix}\triangleq\begin{bmatrix}A_{x,1}(0,t)\\\
A_{y,1}(0,t)\\\ …\\\ A_{x,D}(0,t)\\\ A_{y,D}(0,t)\end{bmatrix},$ (2)

and

$\displaystyle\bra{A(0,t)}\triangleq\begin{bmatrix}A_{1}^{*}(0,t)&A_{2}^{*}(0,t)&…&A_{2D-1}^{*}(0,t)&A_{2D}^{*}(0,t)\end{bmatrix}$
(3)
$\displaystyle\triangleq\begin{bmatrix}A_{x,1}^{*}(0,t)&A_{y,1}^{*}(0,t)&…&A_{x,D}^{*}(0,t)&A_{y,D}^{*}(0,t)\end{bmatrix},$

where $A_{i_{3}}(0,t)\triangleq A_{p^{{}^{\prime}},p}(0,t)$ represents the
time domain of the propagated signal in $i_{3}$th polarization-spatial mode
(i.e., $p^{{}^{\prime}}$th polarization mode and $p$th spatial mode) and can
be expressed as

$\displaystyle\ket{A_{i_{3}}(0,t)}=\sum_{i_{1}=-\infty}^{\infty}\sum_{i_{2}=1}^{N_{ch}}\zeta_{\mathbf{i}}W_{Tx}^{i_{2},i_{3}}(t-i_{1}T_{i_{2}})e^{j2\pi
f_{i_{2}}t}\ket{i_{3}},$ (4)

Considering (1), the propagated signal in the frequency domain can be
expressed by

$\displaystyle\ket{\tilde{A}(0,f)}=\sum_{\mathbf{i}}\zeta_{\mathbf{i}}\ket{\tilde{W}_{Tx}^{\mathbf{i}}(f)},$
(5)

where

$\displaystyle\ket{\tilde{W}_{Tx}^{\mathbf{i}}(f)}\triangleq\tilde{W}_{Tx}^{\mathbf{i}}(f)\ket{i_{3}}=\tilde{W}_{Tx}^{i_{2},i_{3}}(f-f_{i_{2}})e^{-j2\pi(f-f_{i_{2}})i_{1}T_{i_{2}}}\ket{i_{3}}.$
(6)

$\displaystyle
n_{\mathbf{i}}(f)=-j\sum_{\mathbf{{\mathbf{k,m,n}}}}\zeta_{\mathbf{k}}^{*}\zeta_{\mathbf{m}}\zeta_{\mathbf{n}}\iiint_{-\infty}^{\infty}\ket{\eta(f,f_{1},f_{2})}\braket{\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{1}+f_{2})}{\tilde{W}_{Tx}^{\mathbf{m}}(f+f_{2})}\braket{g^{\mathbf{i}}(f)}{\tilde{W}_{Tx}^{\mathbf{n}}(f+f_{1})}R_{\mathbf{i}}df_{1}df_{2}df,$
(17) $\displaystyle
n_{\mathbf{i}}(f)=-j\sum_{\mathbf{k,m,n}}\zeta_{\mathbf{k}}^{*}\zeta_{\mathbf{m}}\zeta_{\mathbf{n}}\iiint_{-\infty}^{\infty}\ket{\eta(f,f_{1},f_{2})}\tilde{W}_{Tx}^{{\mathbf{k}}^{*}}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{m}}(f+f_{2})g^{{\mathbf{i}}^{*}}(f)\tilde{W}_{Tx}^{\mathbf{n}}(f+f_{1})\braket{k_{3}}{m_{3}}\braket{i_{3}}{n_{3}}R_{\mathbf{i}}df_{1}df_{2}df,$
(18) $\displaystyle
E[n_{\mathbf{i}}(f)n_{\mathbf{i}}^{*}(f)]=\sum_{\mathbf{kmnljo}}E[\zeta_{\mathbf{k}}^{*}\zeta_{\mathbf{m}}\zeta_{\mathbf{n}}\zeta_{\mathbf{l}}\zeta_{\mathbf{j}}^{*}\zeta_{\mathbf{o}}^{*}]\idotsint_{\infty}^{\infty}\braket{\eta(f,f_{1},f_{2})}{\eta(v,v_{1},v_{2})}\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{m}}(f+f_{2})$
(19) $\displaystyle
g^{\mathbf{i}*}(f)\tilde{W}_{Tx}^{\mathbf{n}}(f+f_{1})\tilde{W}_{Tx}^{\mathbf{j}*}(v+v_{2})\tilde{W}_{Tx}^{\mathbf{l}}(v+v_{1}+v_{2})\tilde{W}_{Tx}^{\mathbf{o}*}(v+v_{1})g^{\mathbf{i}}(v)\braket{k_{3}}{m_{3}}\braket{i_{3}}{n_{3}}\braket{j_{3}}{l_{3}}\braket{o_{3}}{i_{3}}R_{\mathbf{i}}df_{1}df_{2}dfdv_{1}dv_{2}dv.$
$\displaystyle\sigma_{EGN,i_{2},p}^{2}=\sum_{s=1}^{N_{s}}\prod_{n=1}^{s-1}(G_{n}\leavevmode\nobreak\
L_{n})^{3}\prod_{n=s}^{N_{s}}(G_{s}\leavevmode\nobreak\
L_{s})\sum_{q=1}^{D}\bigg{[}3/4\sum_{k_{2},m_{2},n_{2}}\kappa_{1}^{(k_{2})}\kappa_{1}^{(m_{2})}\kappa_{1}^{(n_{2})}P_{k_{2},q}P_{m_{2},q}P_{n_{2},p}X_{l_{2},p}^{a}(k_{2},m_{2},n_{2},q)+$
(25) $\displaystyle
1/4\sum_{k_{2},n_{2}}\kappa_{2}^{(k_{2})}\kappa_{1}^{(n_{2})}(P_{k_{2},q}^{2}P_{n_{2},p}5X_{l_{2},p}^{b}(k_{2},k_{2},n_{2},q)+P_{k_{2},p}P_{k_{2},q}P_{n_{2},q}X_{l_{2},p}^{c}(k_{2},n_{2},k_{2},q))+$
$\displaystyle
1/4\sum_{n_{2}}\kappa_{3}^{(n_{2})}P_{n_{2},q}^{2}P_{n_{2},p}X_{l_{2},p}^{d}(n_{2},n_{2},n_{2},q)\bigg{]}$

##  3 EGN model formulation

The Manakov equation for the considered MDM-WDM link can be expressed as [7]

$\displaystyle\frac{\partial\ket{A_{i_{3}}(z,t)}}{\partial
z}=\mathcal{L}+\mathcal{N}$ (7)

where the linear and nonlinear effects of MMF can be shown as [3]

$\displaystyle\mathcal{L}=-\frac{\alpha_{i_{3}}}{2}\ket{A_{i_{3}}(z,t)}+j{\beta_{0}}_{i_{3}}\ket{A_{i_{3}}(z,t)}-{\beta_{1}}_{i_{3}}\frac{\partial\ket{A_{i_{3}}(z,t)}}{\partial
t}$ (8)
$\displaystyle-j\frac{{\beta_{2}}_{i_{3}}}{2}\frac{\partial^{2}\ket{A_{i_{3}}(z,t)}}{\partial
t^{2}}-\frac{{\beta_{3}}_{i_{3}}}{6}\frac{\partial^{3}\ket{A_{i_{3}}(z,t)}}{\partial
t^{3}}$
$\displaystyle\mathcal{N}=j\sum_{k_{3},m_{3},n_{3}=1}^{2D}\tilde{\gamma}_{i_{3}k_{3}m_{3}n_{3}}\bigg{(}\braket{A_{n_{3}}(z,t)}{A_{m_{3}}(z,t)}\ket{A_{k_{3}}^{*}(z,t)}$
(9)
$\displaystyle+\braket{A_{k_{3}}^{*}(z,t)}{A_{m_{3}}(z,t)}\ket{A_{n_{3}}(z,t)}\bigg{)},$

$\alpha_{i_{3}}$ is the attenuation of $i_{3}$th polarization-spatial mode.
$\tilde{\gamma}_{i_{3}i_{3}i_{3}i_{3}}=\frac{8}{9}\gamma
f_{i_{3}i_{3}i_{3}i_{3}}$, and
$\tilde{\gamma}_{i_{3}k_{3}m_{3}n_{3}}=\frac{4}{3}\gamma
f_{i_{3}k_{3}m_{3}n_{3}}$ where $\gamma$ is the Kerr nonlinearity coefficient,
with
$f_{i_{3}k_{3}m_{3}n_{3}}=\frac{A_{eff}}{\sqrt{I_{i_{3}}I_{k_{3}}I_{m_{3}}I_{n_{3}}}}\iint{F_{i_{3}}(x,y)F_{k_{3}}(x,y)F_{m_{3}}(x,y)F_{n_{3}}(x,y)dxdy}$
is the nonlinear coupling coefficient between modes $i_{3}$, $k_{3}$, $m_{3}$
and $n_{3}$ where $F_{i_{3}}(x,y)$ is the spatial profile of $i_{3}$th mode,
$I_{i_{3}}=\iint{F_{i_{3}}^{2}(x,y)dxdy}$, $A_{eff}$ is the effective area of
the fundamental mode [8].

The first-order perturbation approximation of Manakov equation solution
expresses the received signal as [7]

$\displaystyle\ket{A(z,t)}\simeq
e^{\mathcal{L}z}\ket{A(0,t)}+\int_{0}^{z}e^{\mathcal{L}(z-\xi)}\mathcal{N}\bigg{(}e^{\mathcal{L}\xi}\ket{A(0,t)}\bigg{)}d\xi,$
(10)

After compensating the linear effects of MMF using the MIMO DSP, the received
signal can be expressed as

$\displaystyle\ket{A_{R}}\simeq\ket{A(0,t)}+\int_{0}^{z}e^{-\mathcal{L}\xi}\mathcal{N}(e^{\mathcal{L}\xi}\ket{A(0,t)})d\xi,$
(11)

The $\mathcal{L}$ can be described in frequency domain by its Fourier
transform $\mathcal{F}(e^{\mathcal{L}z})=\ket{e^{\nu(z,f)}}$ where

$\displaystyle\nu_{i_{3}}(z,f)=-\int_{0}^{z}(\alpha_{i_{3}}(\xi)+j\beta_{i_{3}}(\xi,f))d\xi,$
(12)

Moreover, $\beta_{i_{3}}(z,f)$ can be calculated as

$\displaystyle\beta_{i_{3}}(z,f)=\beta_{0_{i_{3}}}+\beta_{1_{i_{3}}}(2\pi
f)+\frac{\beta_{2_{i_{3}}}}{2}(2\pi f)^{2}+\frac{\beta_{3_{i_{3}}}}{6}(2\pi
f)^{3}.$ (13)

The second term in (11) represents the nonlinear noise which by considering
(12) can be simplified in the frequency domain as

$\displaystyle\ket{\tilde{n}(f)}=-j\iint_{-\infty}^{-\infty}\ket{\eta(f,f_{1},f_{2})}\braket{\tilde{A}(f+f_{1}+f_{2})}{\tilde{A}(f+f_{2})}$
(14) $\displaystyle\ket{\tilde{A}(f+f_{1})}df_{1}df_{2},$

where

$\displaystyle\eta_{i_{3}}(f,f_{1},f_{2})\triangleq\sum_{k_{3},m_{3},n_{3}}\tilde{\gamma}_{i_{3},k_{3},m_{3},n_{3}}$
(15)
$\displaystyle\int_{0}^{z}e^{\nu_{n_{3}}(\xi,f+f_{1})+\nu_{m_{3}}(\xi,f+f_{2})+\nu_{k_{3}}^{*}(\xi,f+f_{1}+f_{2})-\nu_{i_{3}}(\xi,f)}d\xi.$

The matched filtering deployment on the received signal can be described by

$\displaystyle\int_{-\infty}^{\infty}\braket{g^{\mathbf{i}}(f)}{\tilde{A}(f)}R_{\mathbf{i}}df,$
(16)

where $\ket{g^{\mathbf{i}}(f)}=g^{\mathbf{i}}(f)\ket{i_{3}}$ is the spectral
shape of transmitted pulse on $i_{2}$th channel and $i_{3}$th polarization-
spatial mode which has been normalized such that
$\int_{-\infty}^{+\infty}g^{\mathbf{i}}(f)df=1$, $R_{\mathbf{i}}$ is the
symbol rate of transmitted pulse on $i_{2}$th channel and $i_{3}$th
polarization-spatial mode. Accordingly, the nonlinear noise takes the form of
(17) at the receiver which can be written as (18).

The variance of the nonlinear noise is derived in (19), which depends on six
infinite integration/summations where the following Poisson summation helps
dropping some interactions/summations:

$\displaystyle\sum_{k=-\infty}^{\infty}e^{jkf2\pi
T}=T\sum_{k=-\infty}^{\infty}\delta(f-k/T).$ (20)

By considering the $sinc$ pulses with finite bandwidth interacting with only
one Dirac’s delta, summation over time index can be in (19). Furthermore,
equating the arguments with equal atoms results in

$\displaystyle f+f_{1}+f_{2}=v+v_{1}+v_{2}$ (21) $\displaystyle
f+f_{2}=v+v_{2}$ $\displaystyle f+f_{1}=v+v_{1},$

and accordingly

$\displaystyle f=v$ (22) $\displaystyle f_{1}=v_{1}$ $\displaystyle
f_{2}=v_{2},$

which yields to simplify three integrals with $v,v_{1}$, and $v_{2}$.

$\displaystyle
M_{i_{2},p}=\Bigg{(}P_{i_{2},p}\prod_{n=1}^{N_{s}}(G_{n}\leavevmode\nobreak\
L_{n})/\bigg{(}\prod_{n=1}^{N_{s}}(G_{n}\leavevmode\nobreak\
L_{n})(F(G_{BA}-1)h\nu B_{i_{2}})+\sum_{s=1}^{N_{s}}[(F(G_{s}-1)h\nu
B_{i_{2}})\prod_{n=s+1}^{N_{s}}(G_{n}\leavevmode\nobreak\
L_{n})]+\sum_{s=1}^{N_{s}}\prod_{n=1}^{s-1}(G_{n}\leavevmode\nobreak\
L_{n})^{3}$ (28) $\displaystyle\prod_{n=s}^{N_{s}}(G_{s}\leavevmode\nobreak\
L_{s})\sum_{q=1}^{D}\bigg{[}3/4\sum_{k_{2},m_{2},n_{2}}\kappa_{1}^{(k_{2})}\kappa_{1}^{(m_{2})}\kappa_{1}^{(n_{2})}P_{k_{2},q}P_{m_{2},q}P_{n_{2},p}X_{l_{2},p}^{a}(k_{2},m_{2},n_{2},q)+1/4\sum_{k_{2},n_{2}}\kappa_{2}^{(k_{2})}\kappa_{1}^{(n_{2})}(P_{k_{2},q}^{2}P_{n_{2},p}$
$\displaystyle
5X_{l_{2},p}^{b}(k_{2},k_{2},n_{2},q)+P_{k_{2},p}P_{k_{2},q}P_{n_{2},q}X_{l_{2},p}^{c}(k_{2},n_{2},k_{2},q))+1/4\sum_{n_{2}}\kappa_{3}^{(n_{2})}P_{n_{2},q}^{2}P_{n_{2},p}X_{l_{2},p}^{d}(n_{2},n_{2},n_{2},q)\bigg{]}+\sigma_{RxN}^{2}\bigg{)}\Bigg{)}/SNR^{req}_{i_{2},p},$

The Kerr nonlinearity is cubic, therefore, the product
$E[n_{\mathbf{i}}(f)n_{\mathbf{i}}^{*}(f)]$ depends on the product of six
atoms. Note that only combinations with an equal number of conjugate/non-
conjugate pairs are non-zero. In addition, it should be noted that
$\zeta_{\mathbf{i}}$ are independent and identically distributed random
variables with zero mean and unit variance. The considered CPE removes the
average phase at the receiver ($\phi$) [7]. In a perturbative frame, this
corresponds to work with the following nonlinear interference noise

$\displaystyle
n_{\mathbf{i}}^{{}^{\prime}}=n_{\mathbf{i}}+j\phi\zeta_{\mathbf{i}}.$ (23)

The nonlinear noise variance in EGN model can be interpreted as a summation
over the Second-Order Noise (SON) which is usually called the GN contribution,
Fourth-Order Noise (FON), and Higher-Order Noise (HON) variances, $i.e.$

$\displaystyle\sigma_{EGN}^{2}=\sigma_{GN}^{2}+\sigma_{FON}^{2}+\sigma_{HON}^{2}.$
(24)

$\displaystyle\max\limits_{G_{i},P_{i_{2},p}}\leavevmode\nobreak\
\min\limits_{i_{2},p,i}\leavevmode\nobreak\
P_{i_{2},p}\prod_{n=1}^{N_{s}}(G_{n}\leavevmode\nobreak\
L_{n})/\bigg{(}\prod_{n=1}^{N_{s}}(G_{n}\leavevmode\nobreak\
L_{n})(F(G_{BA}-1)h\nu B_{i_{2}})+\sum_{s=1}^{N_{s}}[(F(G_{s}-1)h\nu
B_{i_{2}})\prod_{n=s+1}^{N_{s}}(G_{n}\leavevmode\nobreak\
L_{n})]+\sum_{s=1}^{N_{s}}\prod_{n=1}^{s-1}(G_{n}\leavevmode\nobreak\
L_{n})^{3}$ (29) $\displaystyle\prod_{n=s}^{N_{s}}(G_{s}\leavevmode\nobreak\
L_{s})\sum_{q=1}^{D}\bigg{[}3/4\sum_{k_{2},m_{2},n_{2}}\kappa_{1}^{(k_{2})}\kappa_{1}^{(m_{2})}\kappa_{1}^{(n_{2})}P_{k_{2},q}P_{m_{2},q}P_{n_{2},p}X_{l_{2},p}^{a}(k_{2},m_{2},n_{2},q)+1/4\sum_{k_{2},n_{2}}\kappa_{2}^{(k_{2})}\kappa_{1}^{(n_{2})}(P_{k_{2},q}^{2}P_{n_{2},p}$
$\displaystyle
5X_{l_{2},p}^{b}(k_{2},k_{2},n_{2},q)+P_{k_{2},p}P_{k_{2},q}P_{n_{2},q}X_{l_{2},p}^{c}(k_{2},n_{2},k_{2},q))+1/4\sum_{n_{2}}\kappa_{3}^{(n_{2})}P_{n_{2},q}^{2}P_{n_{2},p}X_{l_{2},p}^{d}(n_{2},n_{2},n_{2},q)\bigg{]}+\sigma_{RxN}^{2}\bigg{)}\frac{1}{SNR^{req}_{i_{2},p}}$
$\displaystyle
s.t.\left\\{\begin{array}[]{ll}P_{i_{2},p,s}=P_{i_{2},p,s-1}G_{s}\leavevmode\nobreak\
L_{s}\\\ \sum_{i_{2},p}P_{i_{2},p}\prod_{n=1}^{s-1}(G_{n}\leavevmode\nobreak\
L_{n})\leq P_{sat}^{MM-EDFA_{s}}\\\ G_{s}\leq G_{max}^{MM-
EDFA_{s}}\end{array}\right.$
$\displaystyle\min\limits_{P_{i_{2},p},G_{i}}\leavevmode\nobreak\
\max\limits_{i_{2},p,i}\leavevmode\nobreak\
SNR^{req}_{i_{2},p}\bigg{(}\prod_{n=1}^{N_{s}}(G_{n}\leavevmode\nobreak\
L_{n})(F(G_{BA}-1)h\nu B_{i_{2}})+\sum_{s=1}^{N_{s}}[(F(G_{s}-1)h\nu
B_{i_{2}})\prod_{n=s+1}^{N_{s}}(G_{n}\leavevmode\nobreak\
L_{n})]+\sum_{s=1}^{N_{s}}\prod_{n=1}^{s-1}(G_{n}\leavevmode\nobreak\
L_{n})^{3}\prod_{n=s}^{N_{s}}(G_{s}\leavevmode\nobreak\ L_{s})$ (30)
$\displaystyle\sum_{q=1}^{D}\bigg{[}3/4\sum_{k_{2},m_{2},n_{2}}\kappa_{1}^{(k_{2})}\kappa_{1}^{(m_{2})}\kappa_{1}^{(n_{2})}P_{k_{2},q}P_{m_{2},q}P_{n_{2},p}X_{l_{2},p}^{a}(k_{2},m_{2},n_{2},q)+1/4\sum_{k_{2},n_{2}}\kappa_{2}^{(k_{2})}\kappa_{1}^{(n_{2})}(P_{k_{2},q}^{2}P_{n_{2},p}5X_{l_{2},p}^{b}(k_{2},k_{2},n_{2},q)$
$\displaystyle+P_{k_{2},p}P_{k_{2},q}P_{n_{2},q}X_{l_{2},p}^{c}(k_{2},n_{2},k_{2},q))+1/4\sum_{n_{2}}\kappa_{3}^{(n_{2})}P_{n_{2},p}^{2}P_{n_{2},q}X_{l_{2},p}^{d}(n_{2},n_{2},n_{2},q)\bigg{]}+\sigma_{RxN}^{2}\bigg{)}/\bigg{(}P_{i_{2},p}\prod_{n=1}^{N_{s}}(G_{n}\leavevmode\nobreak\
L_{n})\bigg{)}$ $\displaystyle
s.t.\left\\{\begin{array}[]{ll}P_{i_{2},p,s}=P_{i_{2},p,s-1}G_{s}\leavevmode\nobreak\
L_{s}\\\ \sum_{i_{2},p}P_{i_{2},p}\prod_{n=1}^{s-1}(G_{n}\leavevmode\nobreak\
L_{n})\leq P_{sat}^{MM-EDFA_{s}}\\\ G_{s}\leq G_{max}^{MM-
EDFA_{s}}\end{array}\right.$

It is shown in Appendix A that the nonlinear noise variance of the whole link
becomes equal to (25). Therefore, the nonlinear noise variance of the whole
link becomes as (25) where $P_{i_{2},p}\triangleq P_{i_{2},p,1}$ is the
launched power at the first span, and
$P_{i_{2},p,s}=P_{i_{2},p,s-1}(G_{s}\leavevmode\nobreak\ L_{s})$ is the
launched power at the $s$th span. Note that the power at each span input is
the multiplication of the power of the previous span by loss and MM-EDFA gain
of that span.

##  4 Problem statement

The SNR margin available between the operating conditions and the error
correction threshold allows for system aging, fiber repairs, and transient
events. Extra margin ensures that service level agreements have been met for
the lifetime of the system. The channel/mode with the lowest SNR margin is the
most failure probability, thus minimum SNR margin maximization would minimize
this probability. Based on the proposed notations in the previous section, the
SNR of $i_{2}$th channel and $p$th mode can be expressed as [10]

$SNR_{i_{2},p}=\frac{P_{i_{2},p}\prod_{n=1}^{N_{s}}(G_{n}\leavevmode\nobreak\
L_{n})}{\sigma_{ASE}^{2}+\sigma_{EGN,i_{2},p}^{2}+\sigma_{RxN}^{2}},$ (26)

where $\sigma_{RxN}^{2}$ is the receiver noise power. Moreover, the variance
of ASE noise in the receiver can be expressed as

$\displaystyle\sigma_{ASE}^{2}=\prod_{n=1}^{N_{s}}(G_{n}\leavevmode\nobreak\
L_{n})(F(G_{BA}-1)h\nu B_{i_{2}})+$ (27)
$\displaystyle\sum_{s=1}^{N_{s}}[(F(G_{s}-1)h\nu
B_{i_{2}})\prod_{n=s+1}^{N_{s}}(G_{n}\leavevmode\nobreak\ L_{n})],$

$\displaystyle\min\limits_{\hat{P}_{i_{2},p},g_{i}}\leavevmode\nobreak\
\max\limits_{i_{2},p,i}\leavevmode\nobreak\
log(SNR^{req}_{i_{2},p})+log\bigg{(}\prod_{n=1}^{N_{s}}(e^{g_{n}}L_{n})(F(G_{BA}-1)h\nu
B_{i_{2}})+\sum_{s=1}^{N_{s}}[(F(e^{g_{s}}-1)h\nu
B_{i_{2}})\prod_{n=s+1}^{N_{s}}(e^{g_{n}}L_{n})]+\sum_{q=1}^{D}\sum_{k_{2},m_{2},n_{2}=1}^{N_{ch}}$
(31)
$\displaystyle\kappa_{1}^{(k_{2})}\kappa_{1}^{(m_{2})}\kappa_{1}^{(n_{2})}e^{\hat{P}_{k_{2},p}+\hat{P}_{m_{2},q}+\hat{P}_{n_{2},q}}3X_{l_{2},p}^{a}(k_{2},m_{2},n_{2},q)+\sum_{q=1}^{D}\sum_{k_{2},n_{2}=1}^{N_{ch}}\kappa_{2}^{(k_{2})}\kappa_{1}^{(n_{2})}(e^{2\hat{P}_{k_{2},q}+\hat{P}_{n_{2},p}}5X_{l_{2},p}^{b}(k_{2},k_{2},n_{2},q)+e^{\hat{P}_{k_{2},p}+\hat{P}_{k_{2},q}+\hat{P}_{n_{2},q}}$
$\displaystyle
X_{l_{2},p}^{c}(k_{2},n_{2},k_{2},q))+\sum_{q=1}^{D}\sum_{n_{2}=1}^{N_{ch}}\kappa_{3}^{(n_{2})}e^{2\hat{P}_{n_{2},q}+\hat{P}_{n_{2},p}}X_{l_{2},p}^{d}(n_{2},n_{2},n_{2},q)+\sigma_{RxN}^{2}\bigg{)}-(\hat{P}_{i_{2},p}+\sum_{n=1}^{N_{s}}g_{n}log(L_{n}))$
$\displaystyle
s.t.\left\\{\begin{array}[]{ll}\hat{P}_{i_{2},p}^{s}=\hat{P}_{i_{2},p}^{s-1}+g_{s}+log(L_{s})\\\
\sum_{i_{2},p}e^{\hat{P}_{i_{2},p}}\prod_{n=1}^{s-1}(e^{g_{n}}L_{n})\leq
P_{sat}^{MM-EDFA_{s}}\\\ g_{s}\leq log(G_{max}^{MM-
EDFA_{s}})\end{array}\right.$
$\displaystyle\min\limits_{\beta,\hat{P}_{i_{2},p},g_{i}}\leavevmode\nobreak\
\beta$ (32) $\displaystyle
s.t.\left\\{\begin{array}[]{ll}log(SNR^{req}_{i_{2},p})+log\bigg{(}\prod_{n=1}^{N_{s}}(e^{g_{n}}L_{n})(F(G_{BA}-1)h\nu
B_{i_{2}})+\sum_{s=1}^{N_{s}}[(F(e^{g_{s}}-1)h\nu
B_{i_{2}})\prod_{n=s+1}^{N_{s}}(e^{g_{n}}L_{n})]+\sum_{q=1}^{D}\sum_{k_{2},m_{2},n_{2}=1}^{N_{ch}}\\\
\kappa_{1}^{(k_{2})}\kappa_{1}^{(m_{2})}\kappa_{1}^{(n_{2})}e^{\hat{P}_{k_{2},p}+\hat{P}_{m_{2},q}+\hat{P}_{n_{2},q}}3X_{l_{2},p}^{a}(k_{2},m_{2},n_{2},q)+\sum_{q=1}^{D}\sum_{k_{2},n_{2}=1}^{N_{ch}}\kappa_{2}^{(k_{2})}\kappa_{1}^{(n_{2})}(e^{2\hat{P}_{k_{2},q}+\hat{P}_{n_{2},p}}5X_{l_{2},p}^{b}(k_{2},k_{2},n_{2},q)\\\
+e^{\hat{P}_{k_{2},p}+\hat{P}_{k_{2},q}+\hat{P}_{n_{2},q}}X_{l_{2},p}^{c}(k_{2},n_{2},k_{2},q))+\sum_{q=1}^{D}\sum_{n_{2}=1}^{N_{ch}}\kappa_{3}^{(n_{2})}e^{2\hat{P}_{n_{2},q}+\hat{P}_{n_{2},p}}X_{l_{2},p}^{d}(n_{2},n_{2},n_{2},q)+\sigma_{RxN}^{2}\bigg{)}-(\hat{P}_{i_{2},p}+\\\
\sum_{n=1}^{N_{s}}g_{n}log(L_{n}))\leq\beta\\\
\hat{P}_{i_{2},p}^{s}=\hat{P}_{i_{2},p}^{s-1}+g_{s}+log(L_{s})\\\
\sum_{i_{2},p}e^{\hat{P}_{i_{2},p}}\prod_{n=1}^{s-1}(e^{g_{n}}L_{n})\leq
P_{sat}^{MM-EDFA_{s}}\\\ g_{s}\leq log(G_{max}^{MM-
EDFA_{s}})\end{array}\right.$
$\displaystyle\min\limits_{\beta,\hat{P}_{l},g_{s}}\leavevmode\nobreak\ \beta$
(33) $\displaystyle
s.t.\left\\{\begin{array}[]{ll}\bigg{[}log(SNR^{req}_{l})+log\bigg{(}\prod_{n=1}^{N_{s}}(e^{g_{n}}L_{n})(F(G_{BA}-1)h\nu
B_{i_{2}})+\sum_{s=1}^{N_{s}}[(F(e^{g_{s}}-1)h\nu
B_{i_{2}})\prod_{n=s+1}^{N_{s}}(e^{g_{n}}L_{n})]+\sum_{l_{1},l_{2},l_{3}=1}^{DN_{ch}}\kappa_{1}^{{{}^{\prime}}^{(l_{1})}}\\\
\kappa_{1}^{{{}^{\prime}}^{(l_{2})}}\kappa_{1}^{{{}^{\prime}}^{(l_{3})}}e^{\hat{P}_{l_{1}}+\hat{P}_{l_{2}}+\hat{P}_{l_{3}}}3H_{l}^{a}(l_{1},l_{2},l_{3})+\sum_{l_{1},l_{2},l_{3}=1}^{DN_{ch}}\kappa_{2}^{{{}^{\prime}}^{(l_{1})}}\kappa_{1}^{{{}^{\prime}}^{(l_{2})}}(e^{2\hat{P}_{l_{1}}+\hat{P}_{l_{2}}}5H_{l}^{b}(l_{1},l_{1},l_{2})+e^{\hat{P}_{l_{1}}+\hat{P}_{l_{2}}+\hat{P}_{l_{3}}}H_{l}^{c}(l_{1},l_{2},l_{3}))\\\
+\sum_{l_{1},l_{2}=1}^{DN_{ch}}\kappa_{3}^{{{}^{\prime}}^{(l_{1})}}e^{2\hat{P}_{l_{1}}+\hat{P}_{l_{2}}}H_{l}^{d}(l_{1},l_{1},l_{2})+\sigma_{RxN}^{2}\bigg{)}-(\hat{P}_{l}+\sum_{n=1}^{N_{s}}g_{n}log(L_{n}))\bigg{]}\leq\beta\\\
\hat{P}_{l}^{s}=\hat{P}_{l}^{s-1}+g_{s}+log(L_{s})\\\
\sum_{l}e^{\hat{P}_{l}}\prod_{n=1}^{s-1}(e^{g_{n}}L_{n})\leq P_{sat}^{MM-
EDFA_{s}}\\\ g_{s}\leq log(G_{max}^{MM-EDFA_{s}})\end{array}\right.$

where $F$ is the amplifier noise figure, $G_{BA}$ is the booster amplifier
gain, $G_{s}$ is the MM-EDFA amplifier gain compensating the loss in each
fiber span, $h$ is Plank constant, and $\nu$ is the central frequency. The SNR
margin of $i_{2}$th channel and $p$th mode can be defined as (28) where
$M_{i_{2},p}$ denotes SNR margin of $i_{2}$th channel and $p$th mode, and
$SNR^{req}_{i_{2},p}$ is the required SNR of $i_{2}$th channel and $p$th mode.
Therefore, the minimum SNR margin maximization problem can be expressed as
(29) where the second constraint means that the total power at the $s$-th MM-
EDFA should be less than saturation power of the $s$-th MM-EDFA, and the third
constraint means that $s$-th MM-EDFA gain should be less than the maximum
possible gain. The optimization problem (29) is equivalent to the optimization
problem (30), as the min-max of a function is equivalent to the max-min of its
inverse. (30) is a non-convex optimization problem. To solve this issue, by
replacing $P_{i_{2},p},G_{i}$ with $e^{\hat{P}_{i_{2},p}},e^{g_{i}}$ in (30)
noting that $log(x)$ is a monotonic function in $x$, we get (31) with the same
minimum as (30). By defining the slack variable $\beta$ (31) can be rewritten
as (32).

We use the gradient descent algorithm in vector form to solve (32). This is
performed by introducing a vector $\mathbf{p}$ of dimension $DN_{ch}$ whose
elements $P_{l};l=1,2,...,DN_{ch}$ are given by $P_{n,m}$,
$n=1,2,...,N_{ch},\leavevmode\nobreak\ m=1,2,...D$. In order to incorporate
the values of $B_{n}$, we use a vector with the same dimension as $\mathbf{p}$
defined as
$\mathbf{B}=[B_{1},B_{1},\ldots,B_{1},B_{2},B_{2},\ldots,B_{2},\ldots,B_{N_{ch}},B_{N_{ch}},\ldots,B_{N_{ch}}]$
in which each $B_{n}$ has been repeated $D$ times. Also let $X$ be a
$N_{ch}\times D\times N_{ch}\times N_{ch}\times N_{ch}\times D$ dimensional
tensor with elements $X_{l_{2},p}^{(.)}(k_{2},m_{2},n_{2},q)$. To match the
latter dimensions with $\mathbf{p}$, we define a $N_{ch}D\times N_{ch}D\times
N_{ch}D\times N_{ch}D$ tensor, $H$, whose elements,
$H_{l}^{(.)}(l_{1},l_{2},l_{3})$, are equal to
$X_{l_{2},p}^{(.)}(k_{2},m_{2},n_{2},q)$ in different subscripts. Thus, (32)
can be expressed as (33). (33) is a convex optimization problem (see Appendix
B) and can be solved by many methods, e.g., Bisection method [11]. The
Bisection method converts the main problem into a feasibility problem by
selecting a region and choosing a candidate for the objective function. The
feasibility problem can be solved using the Lagrangian method [12]. In each
step, the region boundaries are updated based on the obtained solution for the
feasibility problem from the previous step. In this manner, the Bisection
method is converged to the optimum objective.

$\displaystyle\beta+\sum_{l=1}^{DN_{ch}}\lambda_{l}\Bigg{(}\bigg{[}log(SNR^{req}_{l})+log\bigg{(}\prod_{n=1}^{N_{s}}(e^{g_{n}}L_{n})(F(G_{BA}-1)h\nu
B_{i_{2}})+\sum_{s=1}^{N_{s}}[(F(e^{g_{s}}-1)h\nu
B_{i_{2}})\prod_{n=s+1}^{N_{s}}(e^{g_{n}}L_{n})]+\sum_{l_{1},l_{2},l_{3}=1}^{DN_{ch}}\kappa_{1}^{{{}^{\prime}}^{(l_{1})}}\kappa_{1}^{{{}^{\prime}}^{(l_{2})}}$
(34)
$\displaystyle\kappa_{1}^{{{}^{\prime}}^{(l_{3})}}e^{\hat{P}_{l_{1}}+\hat{P}_{l_{2}}+\hat{P}_{l_{3}}}3H_{l}^{a}(l_{1},l_{2},l_{3})+\sum_{l_{1},l_{2},l_{3}=1}^{DN_{ch}}\kappa_{2}^{{{}^{\prime}}^{(l_{1})}}\kappa_{1}^{{{}^{\prime}}^{(l_{2})}}(e^{2\hat{P}_{l_{1}}+\hat{P}_{l_{2}}}5H_{l}^{b}(l_{1},l_{1},l_{2})+e^{\hat{P}_{l_{1}}+\hat{P}_{l_{2}}+\hat{P}_{l_{3}}}H_{l}^{c}(l_{1},l_{2},l_{3}))+\sum_{l_{1},l_{2}=1}^{DN_{ch}}\kappa_{3}^{{{}^{\prime}}^{(l_{1})}}e^{2\hat{P}_{l_{1}}+\hat{P}_{l_{2}}}$
$\displaystyle
H_{l}^{d}(l_{1},l_{1},l_{2})+\sigma_{RxN}^{2}\bigg{)}-(\hat{P}_{l}+\sum_{n=1}^{N_{s}}g_{n}log(L_{n}))\bigg{]}-\beta\Bigg{)}+\sum_{s=1}^{N_{s}}\mu_{s}\bigg{(}\sum_{l}e^{\hat{P}_{l}}\prod_{n=1}^{s-1}(e^{g_{n}}L_{n})-P_{sat}^{MM-
EDFA_{s}}\bigg{)}+\sum_{s=1}^{N_{s}}\nu_{s}(g_{s}-$ $\displaystyle
log(G_{max}^{MM-EDFA_{s}}))$ $\displaystyle\inf\limits_{\hat{P}_{l},g_{i}}{}$
$\displaystyle\beta+\sum_{l=1}^{DN_{ch}}\lambda_{l}\Bigg{(}\bigg{[}log(SNR^{req}_{l})+log\bigg{(}\prod_{n=1}^{N_{s}}(e^{g_{n}}L_{n})(F(G_{BA}-1)h\nu
B_{i_{2}})+\sum_{s=1}^{N_{s}}[(F(e^{g_{s}}-1)h\nu
B_{i_{2}})\prod_{n=s+1}^{N_{s}}(e^{g_{n}}L_{n})]+\sum_{l_{1},l_{2},l_{3}=1}^{DN_{ch}}\kappa_{1}^{{{}^{\prime}}^{(l_{1})}}$
(35)
$\displaystyle\kappa_{1}^{{{}^{\prime}}^{(l_{2})}}\kappa_{1}^{{{}^{\prime}}^{(l_{3})}}e^{\hat{P}_{l_{1}}+\hat{P}_{l_{2}}+\hat{P}_{l_{3}}}3H_{l}^{a}(l_{1},l_{2},l_{3})+\sum_{l_{1},l_{2},l_{3}=1}^{DN_{ch}}\kappa_{2}^{{{}^{\prime}}^{(l_{1})}}\kappa_{1}^{{{}^{\prime}}^{(l_{2})}}(e^{2\hat{P}_{l_{1}}+\hat{P}_{l_{2}}}5H_{l}^{b}(l_{1},l_{1},l_{2})+e^{\hat{P}_{l_{1}}+\hat{P}_{l_{2}}+\hat{P}_{l_{3}}}H_{l}^{c}(l_{1},l_{2},l_{3}))$
$\displaystyle+\sum_{l_{1},l_{2}=1}^{DN_{ch}}\kappa_{3}^{{{}^{\prime}}^{(l_{1})}}e^{2\hat{P}_{l_{1}}+\hat{P}_{l_{2}}}H_{l}^{d}(l_{1},l_{1},l_{2})+\sigma_{RxN}^{2}\bigg{)}-(\hat{P}_{l}+\sum_{n=1}^{N_{s}}g_{n}log(L_{n}))\bigg{]}-\beta\Bigg{)}+\sum_{s=1}^{N_{s}}\mu_{s}\bigg{(}\sum_{l}e^{\hat{P}_{l}}\prod_{n=1}^{s-1}(e^{g_{n}}L_{n})-P_{sat}^{MM-
EDFA_{s}}\bigg{)}$
$\displaystyle+\sum_{s=1}^{N_{s}}\nu_{s}(g_{s}-log(G_{max}^{MM-EDFA_{s}}))$
$\displaystyle\Delta_{\lambda_{l}}=\bigg{[}log(SNR^{req}_{l})+log\bigg{(}\prod_{n=1}^{N_{s}}(e^{g_{n}}L_{n})(F(G_{BA}-1)h\nu
B_{i_{2}})+\sum_{s=1}^{N_{s}}[(F(e^{g_{s}}-1)h\nu
B_{i_{2}})\prod_{n=s+1}^{N_{s}}(e^{g_{n}}L_{n})]+\sum_{l_{1},l_{2},l_{3}=1}^{DN_{ch}}\kappa_{1}^{{{}^{\prime}}^{(l_{1})}}\kappa_{1}^{{{}^{\prime}}^{(l_{2})}}\kappa_{1}^{{{}^{\prime}}^{(l_{3})}}$
(36) $\displaystyle
e^{\hat{P}_{l_{1}}+\hat{P}_{l_{2}}+\hat{P}_{l_{3}}}3H_{l}^{a}(l_{1},l_{2},l_{3})+\sum_{l_{1},l_{2},l_{3}=1}^{DN_{ch}}\kappa_{2}^{{{}^{\prime}}^{(l_{1})}}\kappa_{1}^{{{}^{\prime}}^{(l_{2})}}(e^{2\hat{P}_{l_{1}}+\hat{P}_{l_{2}}}5H_{l}^{b}(l_{1},l_{1},l_{2})+e^{\hat{P}_{l_{1}}+\hat{P}_{l_{2}}+\hat{P}_{l_{3}}}H_{l}^{c}(l_{1},l_{2},l_{3}))+\sum_{l_{1},l_{2}=1}^{DN_{ch}}\kappa_{3}^{{{}^{\prime}}^{(l_{1})}}e^{2\hat{P}_{l_{1}}+\hat{P}_{l_{2}}}$
$\displaystyle
H_{l}^{d}(l_{1},l_{1},l_{2})+\sigma_{RxN}^{2}\bigg{)}-(\hat{P}_{l}+\sum_{n=1}^{N_{s}}g_{n}log(L_{n}))\bigg{]}-\beta$
$\displaystyle\Delta_{\mu_{s}}=\sum_{l}e^{\hat{P}_{l}}\prod_{n=1}^{s-1}(e^{g_{n}}L_{n})-P_{sat}^{MM-
EDFA_{s}}$ $\displaystyle\Delta_{\nu_{s}}=g_{s}-log(G_{max}^{MM-EDFA_{s}})$
Algorithm 1 Bisection Method to Solve Convex Optimization Problem (33).

1:Initialization: upper bound $u=100$, and lower bound $l=-10$

2:$\beta\leftarrow u$

3:Solve convex problem (33) by Lagrangian method

4:if $\hat{\emph{P}}^{*(t)}==NAN$ then break

5:$\beta\leftarrow l$

6:Solve convex problem (33) by Lagrangian method

7:if$\hat{\emph{P}}^{*(t)}==NAN$ then break

8:while $u-l\leq\epsilon$ do

9: $\beta\leftarrow(u+l)/2$

10: Solve convex problem (33) by Lagrangian method

11: if $\hat{\emph{P}}^{*(t)}==NAN$ then $l\leftarrow\beta$ else
$u\leftarrow\beta$

Algorithm 1 summarizes the Bisection method for solving (33). The first step
in this algorithm is defining appropriate upper ( $u$) and lower ( $l$) bounds
for the search region for $\beta$. The defined upper bound is assigned to
$\beta$, and the problem (30) is converted to a feasibility problem. The
feasibility problem is solved using the Lagrangian method. If the defined
upper bound be lower than the optimal solution for $\beta$, the feasibility
problem would not have a solution. In other words, the feasible set is empty
and a higher upper bound should be used. The defined lower bound can be tested
and adjusted in the same manner.

In the second step, the upper bound $(u+l)/2$ is assigned to $\beta$. The same
as the first step, it is tested whether $(u+l)/2$ is the upper or lower bound
of the feasible set. Thereby the upper and lower bounds are updated. The
second step is repeated until convergence.

At each iteration of Algorithm 1, the feasibility problem is solved by the
Lagrange duality method as summarized in Algorithm 2. Furthermore, the second
constraint can be relaxed since it is satisfied by the objective function. The
Lagrangian function of (33) is given by (34) where $\lambda_{l},\mu_{s}$, and
$\nu_{s}\in R^{+}$ are the Lagrangian multipliers. Accordingly, the Lagrangian
dual function of (33) can be expressed as (35).

(35) is a convex problem with respect to $\hat{P}_{l},g_{i}$, since the dual
problem is a convex optimization problem [12]. Note that at each iteration of
Algorithm 2, $\lambda_{l},\mu_{s}$, and $\nu_{s}$ are updated based on the
derivative of (35) with respect to $\lambda_{l},\mu_{s}$, and $\nu_{s}$ which
are shown in (36).

Algorithm 2 Lagrangian Duality Method to Solve the Convex Problem (35).

1:Initialization: iteration counter $t=0$, step size parameter $a>0,b>0,c>0$,
and $\lambda^{(0)}\succeq 0,\mu^{(0)}\succeq 0,\nu^{(0)}\succeq 0$

2:while achieving convergence do

3: Solve convex problem (35) with fixed $\lambda,\mu,and\nu$, and obtain
optimal power $\hat{\emph{P}}^{*(t)}$, and optimal gain $\emph{g}^{*(t)}$

4: $\lambda^{(t+1)}=\bigg{[}\lambda^{(t)}\textendash
a\bigg{(}\Delta_{\lambda}\bigg{)}\bigg{]}^{+}$

5: $\mu^{(t+1)}=\bigg{[}\mu^{(t)}\textendash
b\bigg{(}\Delta_{\mu}\bigg{)}\bigg{]}^{+}$

6: $\nu^{(t+1)}=\bigg{[}\nu^{(t)}\textendash
c\bigg{(}\Delta_{\nu}\bigg{)}\bigg{]}^{+}$

7: update $t=t+1$

Figure 2: a) Nonlinear noise power and b) SNR based on the proposed EGN model
and the SSFM simulation versus launched power of each channel and mode for
MDM-WDM ($D=3,N_{ch}=3$) system.

Remark 1: In Algorithm 2, (35) is solved at each iteration as a function of
$\hat{P}_{l}$ using the gradient descent algorithm which will converge to its
optimum solution due to the convexity of the problem. This procedure is
repeated by Algorithm 1 in the "While loop", by which the minimum SNR margin
is improved successively until convergence to the maximum value. Note that
Algorithm 1 will stop searching while the difference between upper and lower
bound becomes less than $\epsilon$.

##  5 Simulation results

###  5.1 Accuracy of the proposed EGN model formulation

Table 1: Simulation parameters. Coefficient | Symbol | Value
---|---|---
Symbol rate | $R_{n}$ | $32\leavevmode\nobreak\ GBaud$
Channel Spacing | $B_{w}$ | $50\leavevmode\nobreak\ GHz$
Channel bandwidth | $B_{n}$ | $32\leavevmode\nobreak\ GHz$
Span length | $L_{s}$ | $80\leavevmode\nobreak\ km$
Nonlinearity coefficient | $\gamma$ | $1.3\leavevmode\nobreak\ 1/(W.km)$
Center frequency | $\nu$ | $1550\leavevmode\nobreak\ nm$
Noise figure | $F$ | $6\leavevmode\nobreak\ dB$
MM-EDFA maximum gain | $G_{max}^{EDFA}$ | $30\leavevmode\nobreak\ dB$
MM-EDFA saturation power | $P_{sat}^{EDFA}$ | $25\leavevmode\nobreak\ dBm$
Receiver noise power | $\sigma_{Rxn}^{2}$ | $-28\leavevmode\nobreak\ dBm$
Booster amplifier gain | $G_{BA}$ | $20\leavevmode\nobreak\ dB$
Table 2: Nonlinear coupling coefficient $(f_{pq})$ [4]. mq | LP01 | LP11a | LP11b | LP02 | LP21a | LP21b
---|---|---|---|---|---|---
LP01 | 1 | 0.661 | 0.661 | 0.734 | 0.455 | 0.455
LP11a | 0.660 | 1.053 | 1.053 | 0.369 | 0.608 | 0.608
LP11b | 0.660 | 1.053 | 1.053 | 0.369 | 0.608 | 0.608
LP02 | 0.731 | 0.369 | 0.369 | 0.964 | 0.335 | 0.335
LP21a | 0.455 | 0.608 | 0.608 | 0.335 | 0.930 | 0.930
LP21b | 0.455 | 0.608 | 0.608 | 0.335 | 0.930 | 0.930
Table 3: Attenuation ($\alpha_{p}$ ($dB/km$)), and dispersion terms ($\beta_{1_{p}}$($ns/km$), $\beta_{2_{p}}$ ($ps^{2}/km$), and $\beta_{3_{p}}$ ($ps^{3}/km$)) [4]. | LP01 | LP11a | LP11b | LP02 | LP21a | LP21b
---|---|---|---|---|---|---
$\alpha_{p}$ | $0.226$ | $0.226$ | $0.226$ | $0.226$ | $0.226$ | $0.226$
$\beta_{1_{p}}$ | 0 | 6.5 | 6.5 | 9.9 | 12 | 12
$\beta_{2_{p}}$ | 31.86 | 34.8 | 34.8 | 2.93 | 26.51 | 26.51
$\beta_{3_{p}}$ | 0.1452 | 0.1452 | 0.1452 | 0.1452 | 0.1452 | 0.1452
Table 4: Lightpath number propagated by each channel and mode. SMF-WDM |  MDM-single channel
---|---
Channel | lightpath # | Mode | lightpath #
$1,2$ | $L1$ | LP01 | $L1$
$3,4$ | $L2$ | Lp11a | $L2$
$5,6$ | $L3$ | Lp11b | $L3$
$7,8$ | $L4$ | Lp02 | $L4$
$9,10$ | $L5$ | Lp21a | $L5$
$11$ | $L6$ | Lp21b | $L6$
Figure 3: A $4$-node linear network with $6$ lightpaths.

In this section, the signal propagation is simulated by approximating the
output of the Manakov equation using the well-known SSFM method with
logarithmic step-size [13] in the Python/Tensorflow environment. The
simulation parameters and their values are presented in Tables 1, 2, and 3.

Figs 2a and 2b respectively show nonlinear noise power and SNR based on the
proposed EGN model and the SSFM simulation versus launched power of each
channel and mode for MDM-WDM ( $D=3,N_{ch}=3$) system. Note that only the SNR
of the central channel is plotted, and QPSK modulation is considered in SSFM.
As seen in Fig. 2a, the proposed EGN model and SSFM simulation are in close
agreement in all power ranges. However, the GN model overestimates the
nonlinear noise power of SSFM simulation with QPSK modulation, since the
nonlinear noise power of the Gaussian constellation is higher than QPSK. As
seen in Fig. 2b, a close agreement can be seen between the proposed EGN model,
GN model, and SSFM simulation in the linear region, since the linear effects
are dominant in this region. However, a fixed gap appears between the GN model
with the proposed EGN mode and SSFM simulation in the nonlinear region.

###  5.2 Minimum SNR margin maximization

In this section the joint optimized power and gain allocation is performed
based on minimum SNR margin maximization. Three scenarios are considered
including a) best equal power, b) optimized power, and c) joint optimized
power and gain. In the first scenario, equal powers are considered for
different channels and modes with equal MM-EDFA gain in all spans. It is worth
mentioning that the MM-EDFA gain is equal to span loss. In the second
scenario, different powers are allocated to different channels and modes with
equal MM-EDFA gain in all spans. In the third case, allocated powers to each
channel and mode are different. Moreover, the MM-EDFA gain for each span is
optimized separately. The SMF-WDM ( $D=1,N_{ch}=11$), and MDM-single channel (
$D=6,N_{ch}=1$) systems are considered [3], [8].

The point-to-point links often have a homogeneous set where different
channels/modes experience the same interacting channels/modes. However, in
multi-node linear networks, channel/modes may propagate different distances
thus accumulate different nonlinear noise, experience fragmentation/partial
utilization thereby see different interacting channels/modes and observe
different MM-EDFA gains. The $4$-node linear network [15] with $6$ lightpaths
shown in Fig. 3 is considered in this paper for joint power and gain
allocation. The lightpath number propagated by each channel and mode is
presented in Table 4. Note that the BPSK modulation with
$5.5\leavevmode\nobreak\ dB$ required SNR is considered [14].

Figure 4: a) Channel power and b) SNR margin versus channel number, and c) MM-
EDFA gain versus span number, for joint optimized power and gain allocation,
optimized power allocation and best equal power allocation, considering SMF-
WDM system.

Figs 4a and 4b respectively depict channel power and SNR margin versus channel
number, and Fig 4c shows MM-EDFA gain versus span number, for joint optimized
power and gain allocation, optimized power allocation, and best equal power
allocation, considering SMF-WDM system. The obtained minimum SNR margins are
$14.89\leavevmode\nobreak\ dB$ for best equal power allocation,
$15.89\leavevmode\nobreak\ dB$ for optimized power allocation, and
$16.51\leavevmode\nobreak\ dB$ for joint optimized power and gain allocation.
Therefore, the obtained improvements of joint optimized power and gain
allocation over optimized power allocation and best equal power allocation are
$0.62$ and $1.62\leavevmode\nobreak\ dB$, respectively. The joint power and
amplifier gain allocation observes higher degrees of freedom and thus obtains
more improvements over optimized power allocation and best equal power
allocation. Central channel indices observe higher nonlinear noise (lower SNR
margins). Therefore, they should be allocated higher power to have a reliable
link. Moreover, propagating channels in longer lightpath should be allocated
higher power, since the longer the lightpath propagated by a channel, the more
ASE/nonlinear noise power is added to that channel. In joint optimized power
and gain, the last MM-EDFA gain is set to its maximum possible value, this
result can be deduced from SNR formulation where all terms except the receiver
noise are scaled with $G_{N_{s}}$, and to minimize the contribution of the
receiver noise term, the maximum possible value should be chosen for
$G_{N_{s}}$. In joint optimized power and gain, the last MM-EDFA gain is set
to its maximum possible value, this result can be deduced from SNR formulation
where all terms except the receiver noise are scaled with $G_{N_{s}}$, and to
minimize the contribution of the receiver noise term, the maximum possible
value should be chosen for $G_{N_{s}}$.

Figure 5: a) Mode power and b) SNR margin versus mode number, and c) MM-EDFA
gain versus span number, for joint optimized power and gain allocation,
optimized power allocation and best equal power allocation considering MDM-
single channel system.

Figs 5a and 5b respectively depict mode power and SNR margin versus mode
number, and Fig 5c shows MM-EDFA gain versus span number, for joint optimized
power and gain allocation, optimized power allocation, and best equal power
allocation considering MDM-single channel system. The obtained minimum SNR
margins are $15.63\leavevmode\nobreak\ dB$ for best equal power allocation,
$16.74\leavevmode\nobreak\ dB$ for optimized power allocation, and
$17.35\leavevmode\nobreak\ dB$ for joint optimized power and gain allocation.
Therefore, the obtained improvements of joint optimized power and gain
allocation over optimized power allocation and best equal power allocation are
$0.61$ and $1.72\leavevmode\nobreak\ dB$, respectively. The main difference
between different modes of the same channel is their spatial profiles. The
LP11a/b mode has a larger spatial profile and therefore, a higher overlap with
the other modes. Therefore, this mode has more nonlinear noise power and lower
SNR margin. Accordingly, it should be allocated higher power than the other
modes. Neither the allocated powers nor the MM-EDFA gains are not symmetric,
since different modes have different nonlinear noise power which is not
symmetric due to the nonlinear coupling.

##  6 Conclusion

Achieving reliable communication over different channels and modes is one of
the main goals in MDM-WDM networks, and is generally quantified through
minimum SNR margin. In this paper, the EGN model formulation is derived for
MDM-WDM systems for the first time and verified through SSFM simulations.
Based on the proposed EGN model, the joint optimized power and MM-EDMA gain is
proposed considering the minimum SNR margin maximization over different
channels and modes in a multi-node linear network. It is shown that that the
joint optimized power and MM-EDMA gain allocation results in minimum SNR
margin improvement compared to optimized power allocation and the best equal
power allocation. For instance, considering MDM-single channel and single mode
fiber-WDM systems, joint power and gain optimization improves the minimum SNR
margin $0.61\leavevmode\nobreak\ dB$ and $1.72\leavevmode\nobreak\ dB$
compared to optimized power allocation, respectively. Moreover, in these
equations we have

## Appendix A: The nonlinear noise variance of $i_{2}$th channel and $p$th
mode

The GN, FON, and HON terms in EGN model can be expressed as (37), (38), and
(39), respectively.

$\displaystyle\sigma_{GN,i_{2},i_{3}}^{2}=\sum_{{\mathbf{k,m,n}}}\kappa_{1}^{({\mathbf{k}})}\kappa_{1}^{({\mathbf{m}})}\kappa_{1}^{({\mathbf{n}})}\Bigg{(}\iiint_{-\infty}^{\infty}\braket{\eta(f,f_{1},f_{2})}{\eta(f,f_{1},f_{2})}\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{m}}(f+f_{2})g^{\mathbf{i}*}(f)\tilde{W}_{Tx}^{\mathbf{n}}(f+f_{1})\braket{k_{3}}{m_{3}}\braket{i_{3}}{n_{3}}$
(37)
$\displaystyle\tilde{W}_{Tx}^{\mathbf{m}*}(f+f_{2})\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{n}*}(f+f_{1})g^{\mathbf{i}}(f)\braket{m_{3}}{k_{3}}\braket{n_{3}}{i_{3}}df_{1}df_{2}df+\iiint_{-\infty}^{\infty}\braket{\eta(f,f_{1},f_{2})}{\eta(f,f_{1},f_{2})}\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{m}}(f+f_{2})$
$\displaystyle
g^{\mathbf{i}*}(f)\tilde{W}_{Tx}^{\mathbf{n}}(f+f_{1})\braket{k_{3}}{m_{3}}\braket{i_{3}}{n_{3}}\tilde{W}_{Tx}^{\mathbf{n}*}(f+f_{2})\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{m}*}(f+f_{1})g^{\mathbf{i}}(f)\braket{n_{3}}{k_{3}}\braket{m_{3}}{i_{3}}df_{1}df_{2}df\Bigg{)}+\Bigg{(}\iiint_{-\infty}^{\infty}$
$\displaystyle\braket{\eta(f,f_{1},f_{2})}{\eta(f,f_{1},f_{2})}\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{n}}(f+f_{2})g^{\mathbf{i}*}(f)\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{1})\braket{k_{3}}{n_{3}}\braket{i_{3}}{k_{3}}\tilde{W}_{Tx}^{\mathbf{m}*}(f+f_{2})\tilde{W}_{Tx}^{\mathbf{m}}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{n}*}(f+f_{1})g^{\mathbf{i}}(f)$
$\displaystyle\braket{m_{3}}{m_{3}}\braket{n_{3}}{i_{3}}df_{1}df_{2}df+\iiint_{-\infty}^{\infty}\braket{\eta(f,f_{1},f_{2})}{\eta(f,f_{1},f_{2})}\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{n}}(f+f_{2})g^{\mathbf{i}*}(f)\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{1})\braket{k_{3}}{n_{3}}\braket{i_{3}}{k_{3}}\tilde{W}_{Tx}^{\mathbf{n}*}(f+f_{2})$
$\displaystyle\tilde{W}_{Tx}^{\mathbf{m}}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{m}*}(f+f_{1})g^{\mathbf{i}}(f)\braket{n_{3}}{m_{3}}\braket{m_{3}}{i_{3}}df_{1}df_{2}df\Bigg{)}+\Bigg{(}\iiint_{-\infty}^{\infty}\braket{\eta(f,f_{1},f_{2})}{\eta(f,f_{1},f_{2})}\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{2})g^{\mathbf{i}*}(f)$
$\displaystyle\tilde{W}_{Tx}^{\mathbf{m}}(f+f_{1})\braket{k_{3}}{k_{3}}\braket{i_{3}}{m_{3}}\tilde{W}_{Tx}^{\mathbf{m}*}(f+f_{2})\tilde{W}_{Tx}^{\mathbf{m}}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{n}*}(f+f_{1})g^{\mathbf{i}}(f)\braket{m_{3}}{m_{3}}\braket{n_{3}}{i_{3}}df_{1}df_{2}df+\iiint_{-\infty}^{\infty}\braket{\eta(f,f_{1},f_{2})}{\eta(f,f_{1},f_{2})}$
$\displaystyle\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{2})g^{\mathbf{i}*}(f)\tilde{W}_{Tx}^{\mathbf{m}}(f+f_{1})\braket{k_{3}}{k_{3}}\braket{i_{3}}{m_{3}}\tilde{W}_{Tx}^{\mathbf{n}*}(f+f_{2})\tilde{W}_{Tx}^{\mathbf{m}}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{m}*}(f+f_{1})g^{\mathbf{i}}(f)\braket{n_{3}}{m_{3}}\braket{m_{3}}{i_{3}}df_{1}df_{2}df\Bigg{)}$
$\displaystyle\sigma_{FON,i_{2},i_{3}}^{2}=\sum_{{\mathbf{k,m,n}}}\kappa_{2}^{({\mathbf{k}})}\kappa_{1}^{({\mathbf{n}})}\Bigg{(}\iiint_{-\infty}^{\infty}\braket{\eta(f,f_{1},f_{2})}{\eta(f,f_{1},f_{2})}\bigg{|}\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{2})g^{\mathbf{i}*}(f)\tilde{W}_{Tx}^{\mathbf{n}}(f+f_{1})\braket{k_{3}}{k_{3}}\braket{i_{3}}{n_{3}}$
(38)
$\displaystyle+\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{n}}(f+f_{2})g^{\mathbf{i}*}(f)\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{1})\braket{k_{3}}{n_{3}}\braket{i_{3}}{k_{3}}\bigg{|}^{2}df_{1}df_{2}df+\iiint_{-\infty}^{\infty}\braket{\eta(f,f_{1},f_{2})}{\eta(f,f_{1},f_{2})}\bigg{|}\tilde{W}_{Tx}^{\mathbf{n}*}(f+f_{1}+f_{2})$
$\displaystyle\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{2})g^{\mathbf{i}*}(f)\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{1})\braket{n_{3}}{k_{3}}\braket{i_{3}}{k_{3}}\bigg{|}^{2}df_{1}df_{2}df\Bigg{)}+\bigg{[}\iiint_{-\infty}^{\infty}\braket{\eta(f,f_{1},f_{2})}{\eta(f,f_{1},f_{2})}\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{2})g^{\mathbf{i}*}(f)$
$\displaystyle\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{1})\braket{k_{3}}{k_{3}}\braket{i_{3}}{k_{3}}\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{2})\tilde{W}_{Tx}^{\mathbf{n}}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{n}*}(f+f_{1})g^{\mathbf{i}}(f)\braket{n_{3}}{k_{3}}\braket{n_{3}}{i_{3}}df_{1}df_{2}df+\iiint_{-\infty}^{\infty}\braket{\eta(f,f_{1},f_{2})}{\eta(f,f_{1},f_{2})}$
$\displaystyle\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{2})g^{\mathbf{i}*}(f)\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{1})\braket{k_{3}}{k_{3}}\braket{i_{3}}{k_{3}}\tilde{W}_{Tx}^{\mathbf{n}*}(f+f_{2})\tilde{W}_{Tx}^{\mathbf{n}}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{1})g^{\mathbf{i}}(f)\braket{n_{3}}{n_{3}}\braket{i_{3}}{k_{3}}df_{1}df_{2}df\bigg{]}$
$\displaystyle+\bigg{[}\iiint_{-\infty}^{\infty}\braket{\eta(f,f_{1},f_{2})}{\eta(f,f_{1},f_{2})}\tilde{W}_{Tx}^{\mathbf{n}*}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{n}}(f+f_{2})g^{\mathbf{i}*}(f)\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{1})\braket{n_{3}}{n_{3}}\braket{i_{3}}{k_{3}}\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{2})\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{1}+f_{2})$
$\displaystyle\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{1})g^{\mathbf{i}}(f)\braket{k_{3}}{k_{3}}\braket{i_{3}}{k_{3}}df_{1}df_{2}df\bigg{]}+\bigg{[}\iiint_{-\infty}^{\infty}\braket{\eta(f,f_{1},f_{2})}{\eta(f,f_{1},f_{2})}\tilde{W}_{Tx}^{\mathbf{n}*}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{2})g^{\mathbf{i}*}(f)\tilde{W}_{Tx}^{\mathbf{n}}(f+f_{1})\braket{k_{3}}{n_{3}}$
$\displaystyle\braket{i_{3}}{n_{3}}\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{2})\tilde{W}_{Tx}^{\mathbf{k}}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{k}*}(f+f_{1})g^{\mathbf{i}}(f)\braket{k_{3}}{k_{3}}\braket{i_{3}}{k_{3}}df_{1}df_{2}df\bigg{]}$
$\displaystyle\sigma_{HON,i_{2},i_{3}}^{2}=\sum_{{\mathbf{n}}}\kappa_{3}^{({\mathbf{n}})}\iiint_{-\infty}^{\infty}\braket{\eta(f,f_{1},f_{2})}{\eta(f,f_{1},f_{2})}\bigg{|}\tilde{W}_{Tx}^{\mathbf{n}*}(f+f_{1}+f_{2})\tilde{W}_{Tx}^{\mathbf{n}}(f+f_{2})g^{\mathbf{i}*}(f)\tilde{W}_{Tx}^{\mathbf{n}}(f+f_{1})\braket{n_{3}}{n_{3}}\braket{i_{3}}{n_{3}}\bigg{|}^{2}df_{1}df_{2}df$
(39)

$\displaystyle\kappa_{1}=\mu_{2}$ (40)
$\displaystyle\kappa_{2}=\mu_{4}-2\mu_{2}^{2}$
$\displaystyle\kappa_{3}=\mu_{6}-4\mu_{4}\mu_{2}+12\mu_{2}^{3},$
$\displaystyle\mu_{n}=E[|\zeta_{\mathbf{k}}|^{n}].$

Table 5: Valid combinations yielding non-zero $E[\zeta_{\mathbf{k}}^{*}\zeta_{\mathbf{m}}\zeta_{\mathbf{n}}\zeta_{\mathbf{l}}\zeta_{\mathbf{j}}^{*}\zeta_{\mathbf{o}}^{*}]$. $\mathbf{k}^{*}$ | $\mathbf{m}$ | $\mathbf{n}$ | $\mathbf{l}$ | $\mathbf{j}^{*}$ | $\mathbf{o}^{*}$ |
---|---|---|---|---|---|---
$*$ | $*$ | $*$ | $*$ | $*$ | $*$ | $HON$
$*$ | $*$ | $**$ | $*$ | $*$ | $**$ | $FONa$
$*$ | $*$ | $**$ | $*$ | $**$ | $*$ | $FONa$
$*$ | $**$ | $*$ | $*$ | $*$ | $**$ | $FONa$
$*$ | $**$ | $*$ | $*$ | $**$ | $*$ | $FONa$
$*$ | $*$ | $*$ | $**$ | $*$ | $**$ | $FONb$
$*$ | $*$ | $*$ | $**$ | $**$ | $*$ | $FONb$
$**$ | $*$ | $**$ | $*$ | $*$ | $*$ | $FONc$
$**$ | $*$ | $*$ | $**$ | $*$ | $*$ | $FONc$
$**$ | $*$ | $*$ | $**$ | $*$ | $*$ | $FONd$
$*$ | $*$ | $**$ | $+$ | $+$ | $**$ | $GNa$
$*$ | $*$ | $**$ | $+$ | $**$ | $+$ | $GNa$
$*$ | $**$ | $*$ | $+$ | $+$ | $**$ | $GNa$
$*$ | $**$ | $*$ | $+$ | $**$ | $+$ | $GNa$
$*$ | $+$ | $**$ | $*$ | $+$ | $**$ | $GNb$
$*$ | $+$ | $**$ | $*$ | $**$ | $+$ | $GNb$

Table 5 shows valid combinations yielding non-zero
$E[\zeta_{\mathbf{k}}^{*}\zeta_{\mathbf{m}}\zeta_{\mathbf{n}}\zeta_{\mathbf{l}}\zeta_{\mathbf{j}}^{*}\zeta_{\mathbf{o}}^{*}]$
where FONb, FONc, and GNa are the removed terms from the EGN model formulation
due to the CPE assumption. Moreover,
$\sigma_{EGN}^{(i_{3},odd)}(f)=\sigma_{EGN}^{(i_{3},even)}(f)$, and
$\sigma_{EGN}^{(p)}(f)=\sigma_{EGN}^{(i_{3},odd)}(f)+\sigma_{EGN}^{(i_{3},even)}(f)$.
Therefore, the GN, FON, and HON contributions of the nonlinear noise variance
of $i_{2}$th channel and $p$th mode can be written as (41), (42), and (43),
respectively.

$\displaystyle\sigma_{GN,i_{2},p}^{2}=\sum_{q=1}^{D}3\sum_{k_{2},m_{2},n_{2}}\kappa_{1}^{(k_{2})}\kappa_{1}^{(m_{2})}\kappa_{1}^{(n_{2})}\iiint_{-\infty}^{\infty}|\eta(f,f_{1},f_{2})|^{2}\bigg{(}\tilde{G}_{Tx}^{m_{2},q}(f+f_{2})\tilde{G}_{Tx}^{k_{2},q}(f+f_{1}+f_{2})\tilde{G}_{Tx}^{n_{2},p}(f+f_{1})g^{i_{2},p}(f)+$
(41)
$\displaystyle\tilde{G}_{Tx}^{n_{2},q}(f+f_{2})\tilde{G}_{Tx}^{k_{2},q}(f+f_{1}+f_{2})\tilde{G}_{Tx}^{m_{2},p}(f+f_{1})g^{i_{2},p}(f)\bigg{)}df_{1}df_{2}df$
$\displaystyle\sigma_{FON,i_{2},p}^{2}=\sum_{q=1}^{D}\sum_{k_{2},n_{2}}\kappa_{2}^{(k_{2})}\kappa_{1}^{(n_{2})}\Bigg{(}\iiint_{-\infty}^{\infty}|\eta(f,f_{1},f_{2})|^{2}5\tilde{G}_{Tx}^{k_{2},q}(f+f_{1}+f_{2})\tilde{G}_{Tx}^{k_{2},q}(f+f_{2})g^{i_{2},p}(f)\tilde{G}_{Tx}^{n_{2},p}(f+f_{1})df_{1}df_{2}df+$
(42)
$\displaystyle\iiint_{-\infty}^{\infty}|\eta(f,f_{1},f_{2})|^{2}\tilde{G}_{Tx}^{n_{2},q}(f+f_{1}+f_{2})\tilde{G}_{Tx}^{k_{2},q}(f+f_{2})g^{i_{2},p}(f)\tilde{G}_{Tx}^{k_{2},p}(f+f_{1})df_{1}df_{2}df\Bigg{)}$
$\displaystyle\sigma_{HON,i_{2},p}^{2}=\sum_{q=1}^{D}\sum_{n_{2}}\kappa_{3}^{(n_{2})}\iiint_{-\infty}^{\infty}|\eta(f,f_{1},f_{2})|^{2}\tilde{G}_{Tx}^{n_{2},q}(f+f_{1}+f_{2})\tilde{G}_{Tx}^{n_{2},q}(f+f_{2})g^{i_{2},p}(f)\tilde{G}_{Tx}^{n_{2},p}(f+f_{1})df_{1}df_{2}df$
(43)

The power spectral density of the optical launched signal can be written as
$\tilde{G}_{Tx}^{i_{2},p}(f)=P_{i_{2},p}g_{i_{2},p}(f)$. Accordingly, the GN,
FON, and HON contributions of the nonlinear noise variance of $i_{2}$th
channel and $p$th mode can be expressed as (44), (45), and (46), respectively.

$\displaystyle\sigma_{GN,i_{2},p}^{2}=\sum_{q=1}^{D}3/4\sum_{k_{2},m_{2},n_{2}}\kappa_{1}^{(k_{2})}\kappa_{1}^{(m_{2})}\kappa_{1}^{(n_{2})}\iiint_{-\infty}^{\infty}|\eta(f,f_{1},f_{2})|^{2}\bigg{(}P_{k_{2},q}P_{m_{2},q}P_{n_{2},p}g^{m_{2},q}(f+f_{2})g^{k_{2},q}(f+f_{1}+f_{2})$
(44) $\displaystyle
g^{n_{2},p}(f+f_{1})g^{i_{2},p}(f)+P_{k_{2},q}P_{m_{2},p}P_{n_{2},q}g^{n_{2},q}(f+f_{2})g^{k_{2},q}(f+f_{1}+f_{2})g^{m_{2},p}(f+f_{1})g^{i_{2},p}(f)\bigg{)}df_{1}df_{2}df$
$\displaystyle\sigma_{FON,i_{2},p}^{2}=\sum_{q=1}^{D}1/4\sum_{k_{2},n_{2}}\kappa_{2}^{(k_{2})}\kappa_{1}^{(n_{2})}\Bigg{(}P_{k_{2},q}^{2}P_{n_{2},p}5\iiint_{-\infty}^{\infty}|\eta(f,f_{1},f_{2})|^{2}g^{k_{2},q}(f+f_{1}+f_{2})g^{k_{2},q}(f+f_{2})g^{i_{2},p}(f)g^{n_{2},p}(f+f_{1})$
(45) $\displaystyle
df_{1}df_{2}df+P_{k_{2},p}P_{k_{2},q}P_{n_{2},q}\iiint_{-\infty}^{\infty}|\eta(f,f_{1},f_{2})|^{2}g^{n_{2},q}(f+f_{1}+f_{2})g^{k_{2},q}(f+f_{2})g^{i_{2},p}(f)g^{k_{2},p}(f+f_{1})df_{1}df_{2}df\Bigg{)}$
$\displaystyle\sigma_{HON,i_{2},p}^{2}=\sum_{q=1}^{D}1/4\sum_{n_{2}}\kappa_{3}^{(n_{2})}P_{n_{2},q}^{2}P_{n_{2},p}\iiint_{-\infty}^{\infty}|\eta(f,f_{1},f_{2})|^{2}g^{n_{2},q}(f+f_{1}+f_{2})g^{n_{2},q}(f+f_{2})g^{i_{2},p}(f)g^{n_{2},p}(f+f_{1})df_{1}df_{2}df$
(46)

Therefore, the nonlinear noise variance can be written as (47).

$\displaystyle\sigma_{EGN,i_{2},p}^{2}=\sum_{q=1}^{D}3/4\sum_{k_{2},m_{2},n_{2}}\kappa_{1}^{(k_{2})}\kappa_{1}^{(m_{2})}\kappa_{1}^{(n_{2})}P_{k_{2},q}P_{m_{2},q}P_{n_{2},p}X_{i_{2},p}^{a}(k_{2},m_{2},n_{2},q)+1/4\sum_{k_{2},n_{2}}\kappa_{2}^{(k_{2})}\kappa_{1}^{(n_{2})}(P_{k_{2},q}^{2}P_{n_{2},p}5X_{i_{2},p}^{b}(k_{2},k_{2},n_{2},q)$
(47)
$\displaystyle+P_{k_{2},p}P_{k_{2},q}P_{n_{2},q}X_{i_{2},p}^{c}(k_{2},n_{2},k_{2},q))+1/4\sum_{n_{2}}\kappa_{3}^{(n_{2})}P_{n_{2},q}^{2}P_{n_{2},p}X_{i_{2},p}^{d}(n_{2},n_{2},n_{2},q)$
$\displaystyle
where\left\\{\begin{array}[]{ll}X_{i_{2},p}^{a}(k_{2},m_{2},n_{2},q)=\iiint_{-\infty}^{\infty}|\eta(f,f_{1},f_{2})|^{2}g^{m_{2},q}(f+f_{2})g^{k_{2},q}(f+f_{1}+f_{2})g^{n_{2},p}(f+f_{1})g^{i_{2},p}(f)df_{1}df_{2}df\\\
X_{i_{2},p}^{b}(k_{2},k_{2},n_{2},q)=\iiint_{-\infty}^{\infty}|\eta(f,f_{1},f_{2})|^{2}g^{k_{2},q}(f+f_{1}+f_{2})g^{k_{2},q}(f+f_{2})g^{i_{2},p}(f)g^{n_{2},p}(f+f_{1})df_{1}df_{2}df\\\
X_{i_{2},p}^{c}(k_{2},n_{2},k_{2},q)=\iiint_{-\infty}^{\infty}|\eta(f,f_{1},f_{2})|^{2}g^{n_{2},q}(f+f_{1}+f_{2})g^{k_{2},q}(f+f_{2})g^{i_{2},p}(f)g^{k_{2},p}(f+f_{1})df_{1}df_{2}df\\\
X_{i_{2},p}^{d}(n_{2},n_{2},n_{2},q)=\iiint_{-\infty}^{\infty}|\eta(f,f_{1},f_{2})|^{2}g^{n_{2},q}(f+f_{1}+f_{2})g^{n_{2},q}(f+f_{2})g^{i_{2},p}(f)g^{n_{2},p}(f+f_{1})df_{1}df_{2}df\end{array}\right.$

## Appendix B: Convexity proof of (33)

The expression (48)

$\displaystyle log(\prod_{n=1}^{N_{s}}(e^{g_{n}}L_{n})(F(G_{BA}-1)h\nu
B_{i_{2}})+\sum_{s=1}^{N_{s}}[(F(e^{g_{s}}-1)h\nu
B_{i_{2}})\prod_{n=s+1}^{N_{s}}(e^{g_{n}}L_{n})]+\sum_{l_{1},l_{2},l_{3}=1}^{DN_{ch}}\kappa_{1}^{{{}^{\prime}}^{(l_{1})}}\kappa_{1}^{{{}^{\prime}}^{(l_{2})}}\kappa_{1}^{{{}^{\prime}}^{(l_{3})}}e^{\hat{P}_{l_{1}}+\hat{P}_{l_{2}}+\hat{P}_{l_{3}}}3H_{l}^{a}(l_{1},l_{2},l_{3})$
(48)
$\displaystyle+\sum_{l_{1},l_{2},l_{3}=1}^{DN_{ch}}\kappa_{2}^{{{}^{\prime}}^{(l_{1})}}\kappa_{1}^{{{}^{\prime}}^{(l_{2})}}(e^{2\hat{P}_{l_{1}}+\hat{P}_{l_{2}}}5H_{l}^{b}(l_{1},l_{1},l_{2})+e^{\hat{P}_{l_{1}}+\hat{P}_{l_{2}}+\hat{P}_{l_{3}}}H_{l}^{c}(l_{1},l_{2},l_{3}))+\sum_{l_{1},l_{2}=1}^{DN_{ch}}\kappa_{3}^{{{}^{\prime}}^{(l_{1})}}e^{2\hat{P}_{l_{1}}+\hat{P}_{l_{2}}}H_{l}^{d}(l_{1},l_{1},l_{2})\bigg{)}-(\hat{P}_{l}+\sum_{n=1}^{N_{s}}g_{n}log(L_{n}))$

is convex in $\hat{P}_{l},g_{s}$, since $log-sum-exp(x)$ is a convex function
in $x$. The constraint function of (33) the summation of some convex
functions, therefore, it is convex. The objective and constraint functions of
(33) are convex, therefore, (33) is a convex optimization problem.

Disclosures The authors declare no conflicts of interest.

## References

  * [1] F. M. Ferreira, C. S. Costa, S. Sygletos, A. D. Ellis, “Nonlinear performance of few-mode fiber links with intermediate coupling”, Journal of Lightwave Technology, Vol. 37, No. 3, pp. 989-999, 2018.
  * [2] C. Koebele, M. Salsi, L. Milord, R. Ryf, C. Bolle, P. Sillard, S. Bigo and G. Charlet , “40km transmission of five mode division multiplexed data streams at 100Gb/s with low MIMO-DSP complexity”, in Proc. 37th Eur. Conf. Opt. Commun., pp.1-3, Sept. 2011.
  * [3] S. Mumtaz, R. J. Essiambre, G. P. Agrawal, “Nonlinear Propagation in Multimode and Multicore Fibers: Generalization of the Manakov Equations”, Journal of Lightwave Technology, Vol. 31, No. 3, 398-406, 2012.
  * [4] D. Kroushkov, G. Rademacher, and K. Petermann, “Cross Mode Modulation in Multimode Fibers”, Optics Letters, Vol. 38, No. 10, pp. 1642–1644, 2013.
  * [5] A. D. Ellis, N. Mac Suibhne, F. C. G. Gunning, and S. Sygletos, “Expressions for the nonlinear transmission performance of multi-mode optical fiber”, Optics Express, Vol. 21, No. 19, pp. 22 834 – 22 846, 2013.
  * [6] M. A. Amirabadi, M. H. Kahaei, and S. A. Nezamalhosseini, “Novel suboptimal approaches for hyperparameter tuning of deep neural network [under the shelf of optical communication]”, Physical Communication, Vol. 41, pp. 101057, 2020.
  * [7] P. Serena, A. Bononi, “A time-domain extended Gaussian noise model”, Journal of Lightwave Technology, Vol. 33, No. 7, pp. 1459-1472, 2015.
  * [8] G. Rademacher, K. Petermann, “Nonlinear Gaussian noise model for multimode fibers with space-division multiplexing”, Journal of Lightwave Technology, Vol. 34, No. 9, pp. 2280-2287, 2016.
  * [9] G. Rademacher, R. S. Luís, B. J. Puttnam, H. Furukawa, R. Maruyama, K. Aikawa, …, and N. Wada, “Investigation of intermodal four-wave mixing for nonlinear signal processing in few-mode fibers”, IEEE Photonics Technology Letters, Vol. 30, No. 17, pp. 1527-1530, 2018.
  * [10] A. Carena, V. Curri, G. Bosco, P. Poggiolini, F. Forghieri, “Modeling of the impact of nonlinear propagation effects in uncompensated optical coherent transmission links”, Journal of Lightwave Technology, Vol. 30, No. 10, pp. 1524-1539, 2012.
  * [11] F. Daneshfar, E. Hosseini, “Load-frequency control in a deregulated environment based on bisection search”, Iranian Journal of Electrical and Electronic Engineering, Vol. 8, No. 4, pp. 303-310, 2012.
  * [12] S. Boyd, S. P. Boyd, L. Vandenberghe, “Convex optimization”, Cambridge university press, 2004.
  * [13] G. Bosco, A. Carena, V. Curri, R. Gaudino, P. Poggiolini, and S. Benedetto, “Suppression of spurious tones induced by the split-step method in fiber systems simulation”, IEEE Photonics Technology Letters, Vol. 12, No. 5, pp. 489-491, 2000.
  * [14] A. Raeesi, H. Rabbani, L. Beygi, S. Zokaei, “Low-complexity physical layer impairment aware spectrum assignment based on discretized Gaussian model for nonlinear noise in elastic optical networks”, Optics Communications, No. 126011, 2020.
  * [15] H. Rabbani, L. Beygi, S. Ghoshooni, H. Rabbani, E. Agrell, “Quality of Transmission Aware Optical Networking Using Enhanced Gaussian Noise Model”, Journal of Lightwave Technology, Vol. 37, No. 3, pp. 831-838, 2019.