# Frequency Logarithmic Perturbation on the Group-Velocity Dispersion
Parameter with Applications to Passive Optical Networks

Vinícius Oliari, Erik Agrell, , Gabriele Liga, , and Alex Alvarado  V. Oliari,
G. Liga, and A. Alvarado are with the Signal Processing Systems (SPS) Group,
Department of Electrical Engineering, Eindhoven University of Technology, 5600
MB Eindhoven, The Netherlands (e-mails<EMAIL_ADDRESS>a.alvarado@tue.nl). E. Agrell is with the Department of Electrical
Engineering, Chalmers University of Technology, Gothenburg SE-41296, Sweden
(e-mail: agrell@chalmers.se). This work is supported by the Netherlands
Organisation for Scientific Research (NWO) via the VIDI Grant ICONIC (project
number 15685). The work of A. Alvarado has received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (grant agreement No 757791). The work of E. Agrell has
received funding from the Swedish Research Council (VR) under Grant no.
2017-03702. The work of G. Liga is supported by the EuroTechPostdoc programme
under the European Union’s Horizon 2020 research and innovation programme
(Marie Skłodowska-Curie grant agreement No 754462).

###### Abstract

Signal propagation in an optical fiber can be described by the nonlinear
Schrödinger equation (NLSE). The NLSE has no known closed-form solution,
mostly due to the interaction of dispersion and nonlinearities. In this paper,
we present a novel closed-form approximate model for the nonlinear optical
channel, with applications to passive optical networks. The proposed model is
derived using logarithmic perturbation in the frequency domain on the group-
velocity dispersion (GVD) parameter of the NLSE. The model can be seen as an
improvement of the recently proposed regular perturbation (RP) on the GVD
parameter. RP and logarithmic perturbation (LP) on the nonlinear coefficient
have already been studied in the literature, and are hereby compared with RP
on the GVD parameter and the proposed LP model. As an application of the
model, we focus on passive optical networks. For a 20 km PON at 10 Gbaud, the
proposed model improves upon LP on the nonlinear coefficient by 1.5 dB. For
the same system, a detector based on the proposed LP model reduces the uncoded
bit-error-rate by up to 5.4 times at the same input power or reduces the input
power by 0.4 dB at the same information rate.

###### Index Terms:

Channel modeling, chromatic dispersion, Kerr nonlinearity, logarithmic
perturbation, nonlinear Schrödinger equation, optical fiber, regular
perturbation, weakly dispersive regime.

##  I Introduction

Analytical models for optical fiber transmission have been widely studied in
the literature. These models are based on the equations that govern the
optical field propagation: the nonlinear Schrödinger equation (NLSE) [1, Ch.
2] and its variants. The NLSE has no known exact solution for an arbitrary
input waveform. One of the most efficient alternatives for approximated
numerical solutions is the split step Fourier method (SSFM) [2], which
simulates the effects of fiber propagation. On the other hand, to analyse
these effects and design novel transceivers, analytical models are highly
desirable.

0LowHigh0LowHighNLPN[1,3]RP on $\beta_{2}$ [7]RP on $\gamma$ [4,5]Disp.-only
[1]LP on $\gamma$ [5,6]FLP on $\beta_{2}$ ([8] and this
paper)$|\beta_{2}|$$\gamma$ Figure 1: Regions in the $\beta_{2}$ vs $\gamma$
plane where different models are valid. This paper presents a new model for
the region to the left of and under the green curve.

Some of the NLSE’s most used analytical models are only valid under some
restricted values of two fiber parameters: the (Kerr) nonlinear coefficient
$\gamma$ and the group-velocity dispersion (GVD) parameter $\beta_{2}$. The
validity111The validity of a model is defined in this paper as the set of
parameters values in which the model waveform is sufficiently close to the
true waveform obtained from the NLSE. Details on the metric that measures the
distance between the two waveforms are given in Sec. III. of the models with
respect to these parameters is shown in Fig. 1. When $\gamma$ is equal to
zero, the NLSE admits an analytical solution, given by the dispersion-only
model [1, Ch. 3] and represented by the horizontal line at $\gamma=0$ in Fig.
1. When $\beta_{2}$ is set to zero, the NLSE also admits an analytical
solution, given by the nonlinear phase noise model (NLPN) [1, Ch. 4], [3],
represented by the vertical line at $\beta_{2}=0$ in Fig. 1. If both
parameters are nonzero, usually a perturbation approach is used [4].

A perturbation on the nonlinear coefficient $\gamma$ considers the
nonlinearities a minor effect [5, 4] and is accurate for high accumulated
dispersion. For the opposite regime, where the nonlinearities are the major
effect and the dispersion is low but nonzero, a perturbation on the GVD
parameter can be performed. The perturbation techniques usually considered for
the optical fiber propagation are regular perturbation (RP) and logarithmic
perturbation (LP) [4]. It was shown in [6] that LP converges faster to the
true NLSE’s solution than RP. RP and LP on $\gamma$ cover the area under the
yellow and blue curves in Fig. 1, respectively. Recently, we proposed RP on
$\beta_{2}$ in [7] and compared with RP on $\gamma$. RP on $\beta_{2}$
provided more accuracy in the weakly dispersive and highly nonlinear regimes,
represented by the area under the red curve in Fig. 1. A preliminary
investigation of LP on $\beta_{2}$ was reported in our recent work [8], where
LP methods on both $\gamma$ and $\beta_{2}$ were compared on the waveform
level. This paper is an extension of [8], where our first contribution is to
derive the two perturbative models on $\beta_{2}$ presented in [8]. These
models are obtained applying LP in either time or frequency. The latter
approach, which we call frequency logarithmic perturbation (FLP), is the most
accurate of the two for the $\beta_{2}$ expansion.

The FLP on $\beta_{2}$ covers the area under the green dashed curve in Fig. 1,
which improves upon RP on $\beta_{2}$. FLP should not be confused with the
frequency resolved logarithmic perturbation (FRLP) proposed in [9, 10, 11].
FRLP consists in applying LP on the frequency components of the time-domain
signal individually and using the superposition principle to combine them.
However, FRLP still applies the LP principle in the time domain, whereas FLP
we study here applies this principle in the frequency domain. In this paper,
the FLP on $\beta_{2}$ model is compared with the LP on $\gamma$, RP on
$\gamma$, and RP on $\beta_{2}$ models.

The second contribution of this paper is to demonstrate the applicability of
the proposed FLP on $\beta_{2}$ model. To this end, we consider passive
optical networks (PONs), wherein the accumulated dispersion is typically low.
We use the models to estimate optimum decision regions at the receiver. The
method used to obtain the decision regions is compared with the Parzen window
(PW) method, recently proposed in [12] for optical fiber systems. Simulations
are carried out in the C- and O-band, where the latter has a close-to-zero GVD
parameter $\beta_{2}$. At a waveform level (continuous-time), the proposed FLP
on $\beta_{2}$ is shown to outperform the other three models at powers higher
than $7$ dBm for both C- and O-band systems. At a symbol level (discrete-
time), decision regions obtained from LP on $\beta_{2}$ result in bit-error-
rates (BER) more than five times lower than the ones obtained from LP on
$\gamma$. Finally, these decision regions are analysed in a system with
forward-error-correction (FEC).

This paper is organized as follows: Section II provides the mathematical
background for the derivation of the models; Section III compares the models
in both continuous- and discrete-time; and Section IV concludes the paper.

##  II Mathematical Background

The normalized NLSE for noiseless propagation of an optical field $A$ at a
retarded time frame $t$ and distance $z$ for a single-polarization can be
approximated as [5]

$\displaystyle\dfrac{\partial A(t,z)}{\partial z}=$
$\displaystyle-\dfrac{j\beta_{2}}{2}\dfrac{\partial^{2}A(t,z)}{\partial
t^{2}}$ $\displaystyle+j\gamma\mathrm{e}^{-\alpha z}|A(t,z)|^{2}A(t,z),$ (1)

where $\alpha$ is the attenuation coefficient, $\beta_{2}$ the GVD parameter,
and $\gamma$ the nonlinear coefficient. The first term on the right-hand side
of (II) represents the chromatic dispersion. This effect on the waveform $A$
is larger when the bandwidth and/or the fiber length is increased. The last
term represents the Kerr nonlinearity, which has a cubic dependence on the
instantaneous signal power and also increases with the fiber length. The
solution $A$ of (II) can be numerically estimated by the SSFM [1, Ch. 2], [2].
Other effects, such as third-order dispersion (TOD) [1, Ch. 3], are not taken
into account in (II). TOD becomes significant for large bandwidths or when
$\beta_{2}$ is low. This effect will be taken into account for our simulations
in the O-band, although it is not used in the derivation of the models. For
the C-band, we consider solely the effects in (II).

In what follows, we first review three models available in the literature.
Sec. II-A and Sec. II-B describe the RP on $\gamma$ and on $\beta_{2}$,
respectively, while Sec. II-C presents the LP on $\gamma$. Finally, Sec. II-D
introduces the FLP on $\beta_{2}$.

###  II-A Regular Perturbation on the Nonlinear Coefficient

The RP on $\gamma$ was first derived in [5, 4]. To approximate the solution of
(II), the RP method represents the signal by a power series of a certain
coefficient. For the RP on $\gamma$, the RP solution can be written as

$\displaystyle A(t,z)$
$\displaystyle=\sum_{k=0}^{\infty}\gamma^{k}A_{k}^{(\gamma)}(t,z).$ (2)

To obtain the functions $A_{k}^{(\gamma)}$, (2) is substituted into (II) and
the terms multiplied by the $k$-th power of $\gamma$ are equated. An
approximate solution for (II) can be obtained by considering only the
functions $A_{0}^{(\gamma)}$ and $A_{1}^{(\gamma)}$. This approximation is the
first-order RP on $\gamma$ [5, Eqs. (7), (9)], [4, Eq. (12)]

$A(t,z)\approx A_{\text{RP}}^{(\gamma)}(t,z)=A_{0}^{(\gamma)}(t,z)+\gamma
A_{1}^{(\gamma)}(t,z),$ (3)

where

$\displaystyle A_{0}^{(\gamma)}(t,z)$
$\displaystyle=\mathcal{D}_{z}\\{A(\cdot,0)\\}(t),$ (4) $\displaystyle
A_{1}^{(\gamma)}(t,z)$ $\displaystyle=j\int_{0}^{z}\\!\mathrm{e}^{-\alpha
u}\mathcal{D}_{z-u}\left\\{|A_{0}(\cdot,u)|^{2}A_{0}(\cdot,u)\right\\}\\!(t)\text{d}u,$
(5)

and the dispersion operator $\mathcal{D}_{z}$ is

$\displaystyle\mathcal{D}_{z}\\{f\\}(t)$ $\displaystyle\triangleq\left(f\ast
h(\cdot,z)\right)(t),$ (6) $\displaystyle h(t,z)$
$\displaystyle=\dfrac{1}{\sqrt{j2\pi\beta_{2}z}}\mathrm{e}^{-\frac{j}{2\beta_{2}z}t^{2}}.$
(7)

The function $A_{0}^{(\gamma)}$ in (4) is called the dispersion-only solution
of (II). This solution can be seen as a model that is accurate only when the
nonlinear effect is negligible. The first-order RP on $\gamma$ in (3) is
accurate for low nonlinear effects and is illustrated as the yellow curve in
Fig. 1.

###  II-B Regular Perturbation on the GVD Parameter

We recently proposed the RP on $\beta_{2}$ in [7]. The same procedure as in
(2) can be applied by considering $A$ as a power series of $\beta_{2}$, i.e.,

$A(t,z)=\sum_{k=0}^{\infty}\beta_{2}^{k}A_{k}^{(\beta_{2})}(t,z).$ (8)

In analogy to RP on $\gamma$, the functions $A_{k}^{(\beta_{2})}$ are also
obtained by replacing (8) in (II) and equating the terms related to the $k$-th
power of $\beta_{2}$. For the first-order RP, involving the functions
$A_{0}^{(\beta_{2})}$ and $A_{1}^{(\beta_{2})}$, an approximate solution for
$A$ can be obtained by

$A(t,z)\approx
A_{\text{RP}}^{(\beta_{2})}(t,z)=A_{0}^{(\beta_{2})}(t,z)+\beta_{2}A_{1}^{(\beta_{2})}(t,z),$
(9)

where

$A_{0}^{(\beta_{2})}(t,z)=A(t,0)\mathrm{e}^{j\gamma|A(t,0)|^{2}G(z)},$ (10)

and

$A_{1}^{(\beta_{2})}(t,z)=B(t,z)\mathrm{e}^{j\gamma|A(t,0)|^{2}G(z)},$ (11)

with $B$ given by

$\displaystyle B(t,z)=$ $\displaystyle-M(t)z+G_{1}(z)R(t)+G_{2}(z)P(t)$
$\displaystyle-2j\gamma A(t,0)\Re\\{A^{*}(t,0)V(t,z)\\},$ (12) $\displaystyle
V(t,z)=$ $\displaystyle\ G(z)\left[M(t)z-G_{1}(z)R(t)-G_{2}(z)P(t)\right]$
$\displaystyle-G_{1}(z)M(t)+G_{2}(z)R(t)+G_{3}(z)P(t),$ (13) $\displaystyle
M(t)$ $\displaystyle=\dfrac{j}{2}\dfrac{\partial^{2}A(t,0)}{\partial t^{2}},$
(14) $\displaystyle R(t)$
$\displaystyle=\dfrac{\gamma}{2}A(t,0)\dfrac{\partial^{2}|A(t,0)|^{2}}{\partial
t^{2}}+\gamma\dfrac{\partial A(t,0)}{\partial
t}\dfrac{\partial|A(t,0)|^{2}}{\partial t},$ (15) $\displaystyle P(t)$
$\displaystyle=\dfrac{j\gamma^{2}}{2}A(t,0)\left(\dfrac{\partial|A(t,0)|^{2}}{\partial
t}\right)^{2},$ (16) $\displaystyle G_{1}(z)$ $\displaystyle=\dfrac{\alpha
z+\mathrm{e}^{-\alpha z}-1}{\alpha^{2}},$ (17) $\displaystyle G_{2}(z)$
$\displaystyle=\dfrac{2\alpha z+4\mathrm{e}^{-\alpha z}-\mathrm{e}^{-2\alpha
z}-3}{2\alpha^{3}},$ (18) $\displaystyle G_{3}(z)$
$\displaystyle=\dfrac{6\alpha z+18\mathrm{e}^{-\alpha z}-9\mathrm{e}^{-2\alpha
z}+2\mathrm{e}^{-3\alpha z}-11}{6\alpha^{4}}.$ (19)

Analogously to RP on $\gamma$, the function $A_{0}^{(\beta_{2})}$ in (10) for
RP on $\beta_{2}$ is an accurate model when dispersion is negligible, and is
called the NLPN model [3, 7]. The first-order RP on $\beta_{2}$ in (9) is
accurate for low accumulated dispersion and is illustrated as the red curve in
Fig. 1.

###  II-C Logarithmic Perturbation

LP is a mathematical technique similar to RP. LP on $\gamma$ was first
presented in [6, 4] and can be shown to converge faster to the true NLSE’s
solution than RP on $\gamma$. LP functions can be obtained directly by the RP
functions $A_{k}$. For example, following an approach similar to [4], consider
that the signal $A$ can be written as a power series of a coefficient $\theta$
(e.g., $\gamma$ or $\beta_{2}$ as done in (2) and (8)) as

$A(t,z)=\sum_{k=0}^{\infty}\theta^{k}A_{k}^{(\theta)}(t,z),$ (20)

where $A_{k}^{(\theta)}$ is the $k$-th RP function. We now want to express $A$
in its LP version, which takes the form

$A(t,z)=A_{0}^{(\theta)}(t,z)\exp\left(\sum_{k=1}^{\infty}\theta^{k}\psi_{k}^{(\theta)}(t,z)\right),$
(21)

where the function $\psi_{k}^{(\theta)}$ is the $k$-th LP function and
$A_{0}^{(\theta)}$ is the $0$-th order RP function. Representing the
exponential function in (21) by its Taylor expansion yields

$A(t,z)=A_{0}^{(\theta)}(t,z)\sum_{m=0}^{\infty}\dfrac{1}{m!}\left(\sum_{k=1}^{\infty}\theta^{k}\psi_{k}^{(\theta)}(t,z)\right)^{m}.$
(22)

The functions $\psi_{k}^{(\theta)}$ can now be obtained by equating (22) with
(20), and further equating the terms that have the same power of $\theta$. For
example, equating the terms multiplied by $\theta^{1}$, we obtain the first-
order LP function as

$\psi_{1}^{(\theta)}(t,z)=\dfrac{A_{1}^{(\theta)}(t,z)}{A_{0}^{(\theta)}(t,z)}.$
(23)

Equating the terms multiplied by $\theta^{2}$ we can also obtain the second-
order LP function

$\psi_{2}^{(\theta)}(t,z)=\dfrac{A_{2}^{(\theta)}(t,z)}{A_{0}^{(\theta)}(t,z)}-\dfrac{1}{2}\left(\dfrac{A_{1}^{(\theta)}(t,z)}{A_{0}^{(\theta)}(t,z)}\right)^{2}.$
(24)

The function $\psi_{2}^{(\theta)}$ in (24) depends on the RP term
$A_{2}^{(\theta)}$. For RP on $\gamma$, $A_{2}^{(\gamma)}$ is well defined [5,
Eq. (11)], [4, Eq. (12)]. However, for RP on $\beta_{2}$,
$A_{2}^{(\beta_{2})}$ is not known in the literature at the time this paper is
being written. Thus, we will restrict the analysis to first-order LP and RP
only.

TABLE I: Summary of the first-order perturbation methods discussed in this paper Coeff. $\theta$ | Perturbation method
---|---
1 order RP |  1 order LP |  1 order FLP
$A(t,z)=A_{0}^{(\theta)}(t,z)+\theta A_{1}^{(\theta)}(t,z)$ | $A(t,z)=A_{0}^{(\theta)}(t,z)\exp{\left(\theta\psi_{1}^{(\theta)}(t,z)\right)}$ | $\tilde{A}(\omega,z)=\tilde{A}_{0}^{(\theta)}(\omega,z)\exp{\left(\theta\tilde{\zeta}_{1}^{(\theta)}(\omega,z)\right)}$
$\gamma$ | $A_{0}^{(\gamma)}(t,z)$ | $\mathcal{D}_{z}\\{A(\cdot,0)\\}(t)$ | $\psi_{1}^{(\gamma)}(t,z)$ | $\dfrac{A_{1}^{(\gamma)}(t,z)}{A_{0}^{(\gamma)}(t,z)}$ | $\tilde{\zeta}_{1}^{(\gamma)}(\omega,z)$ | $\dfrac{\tilde{A}_{1}^{(\gamma)}(\omega,z)}{\tilde{A}_{0}^{(\gamma)}(\omega,z)}$
$A_{1}^{(\gamma)}(t,z)$ | (5)
$\beta_{2}$ | $A_{0}^{(\beta_{2})}(t,z)$ | $A(t,0)\mathrm{e}^{j\gamma|A(t,0)|^{2}G(z)}$ | $\psi_{1}^{(\beta_{2})}(t,z)$ | $\dfrac{A_{1}^{(\beta_{2})}(t,z)}{A_{0}^{(\beta_{2})}(t,z)}$ | $\tilde{\zeta}_{1}^{(\beta_{2})}(\omega,z)$ | $\dfrac{\tilde{A}_{1}^{(\beta_{2})}(\omega,z)}{\tilde{A}_{0}^{(\beta_{2})}(\omega,z)}$
$A_{1}^{(\beta_{2})}(t,z)$ | (11)

Setting $\theta=\gamma$ or $\theta=\beta_{2}$ in (21) and (23) and truncating
the sum in (21) at $k=1$, we obtain the first-order LP on $\gamma$ and on
$\beta_{2}$, respectively. The first-order LP on $\gamma$ is written as

$A(t,z)\approx
A_{\text{LP}}^{(\gamma)}(t,z)=A_{0}^{(\gamma)}(t,z)\exp\left(\gamma\dfrac{A_{1}^{(\gamma)}(t,z)}{A_{0}^{(\gamma)}(t,z)}\right),$
(25)

where $A_{0}^{(\gamma)}$ and $A_{1}^{(\gamma)}$ are given by (4) and (5),
respectively. The accuracy of $A_{\text{LP}}^{(\gamma)}$ is qualitatively
illustrated by the blue curve in Fig. 1. The first-order LP on $\beta_{2}$ is
similarly obtained as

$A(t,z)\approx
A_{\text{LP}}^{(\beta_{2})}(t,z)=A_{0}^{(\beta_{2})}(t,z)\exp\left(\beta_{2}\dfrac{A_{1}^{(\beta_{2})}(t,z)}{A_{0}^{(\beta_{2})}(t,z)}\right),$
(26)

where $A_{0}^{(\beta_{2})}$ and $A_{1}^{(\beta_{2})}$ are given by (10) and
(11), respectively.

###  II-D Frequency Logarithmic Perturbation

The linearity of (20) with respect to the functions $A_{k}^{(\theta)}$
suggests another approach to obtain a different LP solution. The new approach
consists of performing the same steps as in (20)–(23) in the frequency domain,
which we refer to as FLP. To obtain the FLP solution, we first express (20) in
the frequency domain, i.e.,

$\tilde{A}(\omega,z)=\sum_{k=0}^{\infty}\theta^{k}\tilde{A}_{k}^{(\theta)}(\omega,z),$
(27)

where $\tilde{A}$ represents the Fourier transform222We define the Fourier
transform of a function $A(\cdot,z)$ as
$\tilde{A}(\omega,z)\triangleq\int_{-\infty}^{\infty}\\!A(t,z)\mathrm{e}^{+j\omega
t}\text{d}t$, which depends on the angular frequency $\omega$ and is evaluated
at distance $z$. The inverse Fourier transform of $\tilde{A}(\cdot,z)$ is
$A(t,z)=[1/(2\pi)]\int_{-\infty}^{\infty}\\!\tilde{A}(\omega,z)\mathrm{e}^{-j\omega
t}\text{d}\omega$. of $A$ and $\omega$ is the angular frequency. Analogous to
(21), now $\tilde{A}$ is expressed in its FLP version as

$\tilde{A}(\omega,z)=\tilde{A}_{0}^{(\theta)}(\omega,z)\exp\left(\sum_{k=1}^{\infty}\theta^{k}\tilde{\zeta}_{k}^{(\theta)}(\omega,z)\right),$
(28)

where the function $\tilde{\zeta}_{k}^{(\theta)}$ is the $k$-th FLP function
and $\tilde{A}_{0}^{(\theta)}$ is the Fourier transform of the 0-th order RP
function.

Coherent Tx$\cdot$$\cdot$$\cdot$$\cdot$$\cdot$$\cdot$CoherentRx$20$ km SSMF$1$
km SSMFSplitter $1\\!:\\!64$$0\leq P_{\text{Tx}}\text{ [dBm]}\leq 20$C-band:
$-22.3\leq P_{\text{Rx}}\text{ [dBm]}\leq-2.3$O-band: $-26.5\leq
P_{\text{Rx}}\text{ [dBm]}\leq-6.5$$A(t,0)$$A(t,z)$Fiber Model$A_{M}(t,z)$NSD
(33)Symbol Mapping$x\in\mathcal{S}$01001Symbol Detector$y\in\mathbb{C}$Symbol
Demap.01001$\hat{x}\in\mathcal{S}$BERFiber Channel

Figure 2: PON system setup used for the simulations in this paper. This system
presents low accumulated dispersion and operates in the highly nonlinear
regime for the used range of input powers. The NSD is calculated using the
fiber output $A$, obtained via SSFM, and a fiber model output $A_{M}$. The NSD
exact formula given by (33). The BER is estimated using the bits corresponding
to the input symbols $x$ and the bits corresponding to the estimated input
symbols $\hat{x}$, where $x,\hat{x}\in\mathcal{S}$ and $\mathcal{S}$ is the
set of constellation points. Demap.: demapping.

In complete analogy with the procedure used to obtain (23), the first-order
FLP function is

$\tilde{\zeta}_{1}^{(\theta)}(\omega,z)=\dfrac{\tilde{A}_{1}^{(\theta)}(\omega,z)}{\tilde{A}_{0}^{(\theta)}(\omega,z)},$
(29)

which is used to obtain the first-order FLP on $\gamma$ and on $\beta_{2}$.
The first-order FLP on $\gamma$ is

$\tilde{A}(\omega,z)\approx\tilde{A}_{\text{FLP}}^{(\gamma)}(\omega,z)=\tilde{A}_{0}^{(\gamma)}(\omega,z)\exp\left(\gamma\dfrac{\tilde{A}_{1}^{(\gamma)}(\omega,z)}{\tilde{A}_{0}^{(\gamma)}(\omega,z)}\right),$
(30)

where $\tilde{A}_{0}^{(\gamma)}$ and $\tilde{A}_{1}^{(\gamma)}$ are the
Fourier transforms of (4) and (5), respectively. The first-order FLP on
$\beta_{2}$ is given by

$\tilde{A}(\omega,z)\approx\tilde{A}_{\text{FLP}}^{(\beta_{2})}(\omega,z)=\tilde{A}_{0}^{(\beta_{2})}(\omega,z)\exp\left(\beta_{2}\dfrac{\tilde{A}_{1}^{(\beta_{2})}(\omega,z)}{\tilde{A}_{0}^{(\beta_{2})}(\omega,z)}\right),$
(31)

where $\tilde{A}_{0}^{(\beta_{2})}$ and $\tilde{A}_{1}^{(\beta_{2})}$ are the
Fourier transforms of (10) and (11), respectively. The qualitative behaviour
of the accuracy of $\tilde{A}_{\text{FLP}}^{(\beta_{2})}$ is illustrated as
the green dashed curve in Fig. 1.

The functions $\tilde{\zeta}_{k}^{(\theta)}$ in (28) differ from the Fourier
transform of $\psi_{k}^{(\theta)}$ in (21), since the exponential of the LP
method was applied in the frequency domain. The expressions in (23) and (29)
are the simplest example of this fact, since they do not form, in general, a
Fourier transform pair. Therefore, we expect that these LP-based models result
in a different accuracy for each perturbation coefficient ($\gamma$ or
$\beta_{2}$). As it will be seen later in Sec. III, $A_{\text{LP}}^{(\gamma)}$
is more accurate than $A_{\text{FLP}}^{(\gamma)}$, while
$A_{\text{FLP}}^{(\beta_{2})}$ is more accurate than
$A_{\text{LP}}^{(\beta_{2})}$. We believe that the difference between
$\beta_{2}$ and $\gamma$ when comparing LP and FLP could be explained by the
solution of (II) for only the chromatic dispersion effect or only the Kerr
effect [1]. The solution for the chromatic dispersion effect only is an
exponential in the frequency domain, which resembles the FLP approach.
Similarly, the solution for the Kerr effect only is an exponential in the time
domain, which resembles the LP approach.

Together with the models in the previous sections, we obtained six
perturbation models: two RPs in (3) and (9); two LPs in (25) and (26); and two
FLPs in (30) and (31). Table I summarizes these six first-order perturbation
methods. As shown in Table I, LP and FLP can be obtained using the RP terms.

##  III Simulation Setup and Results

The model presented in this work is validated in a PON transmission scenario
where the accumulated dispersion is expected to be low. Fig. 2 shows the
coherent PON system setup under consideration. The fiber parameters are given
in Table II. As depicted in the figure, we consider a standard single mode
fiber (SSMF) of $20$ km, followed by a splitter of ratio $1\\!:\\!64$ and a
final fiber segment of $1$ km. With this split ratio, the total link loss is
$22.3$ dB for the C-band and $26.5$ dB for the O-band. The fiber input power
$P_{\text{Tx}}$ varies from $0$ to $20$ dBm, which leads to a received power
$P_{\text{Rx}}$ between $-22.3$ and $-2.3$ dBm in the C-band and between
$-26.5$ and $-6.5$ dBm in the O-band. The range of powers was chosen to cover
and go beyond launch powers for typical PON systems according to [13, 14]. All
the results were obtained using randomly generated bits which were mapped into
symbols drawn from a quadrature phase shift keying (QPSK) constellation
$\mathcal{S}=(\pm 1\pm j)/\sqrt{2}$. The coherent transmitter applies pulse
shaping and scales the waveform such that the average transmitted power is
$P_{\text{Tx}}$. The coherent receiver undo the waveform scaling, and then
applies matched filtering and sampling, without chromatic-dispersion
compensation. The symbol rate of the transmitted signal is $10$ Gbaud for both
C- and O-band systems. The considered pulse shape is a root-raised cosine
(RRC), with a roll-off factor of $0.1$. We consider a noiseless scenario since
for the considered bandwidth and received powers, nonlinear distortions
dominate over the shot noise [15, 16].

In the considered system setup, we want to evaluate the impact of
nonlinearities and dispersion on the models. For evaluating the impact of the
nonlinearities, the power was varied as specified before. The effect of the
dispersion in the models is evaluated by comparing the C-band and the O-band
scenarios, which have different $\beta_{2}$ values. We consider the effect of
the TOD in the O-band since in that regime values of $\beta_{2}$ are low.
Therefore, the equation used for the SSFM simulation is slightly different
from (II), i.e., it considers an additional term accounting for TOD [1, Eq.
(3.3.1)], resulting in

$\displaystyle\dfrac{\partial A(t,z)}{\partial z}=$
$\displaystyle-\dfrac{j\beta_{2}}{2}\dfrac{\partial^{2}A(t,z)}{\partial
t^{2}}+\dfrac{\beta_{3}}{6}\dfrac{\partial^{3}A(t,z)}{\partial t^{3}}$
$\displaystyle+j\gamma\mathrm{e}^{-\alpha z}|A(t,z)|^{2}A(t,z).$ (32)

We evaluate the model accuracy on two levels, comparing either the channel
output waveforms or detected symbols. For the former, no receiver is
considered. The accuracy is quantified using the normalized squared deviation
(NSD) metric [5, 7]

$\text{NSD}\triangleq\dfrac{\int_{-\infty}^{\infty}\\!|A_{M}(t,z)-A(t,z)|^{2}\text{d}t}{\int_{-\infty}^{\infty}|A(t,z)|^{2}\text{d}t},$
(33)

where $A_{M}$ is a model output (i.e., $A_{\text{LP}}^{(\gamma)}$ or
$A_{\text{FLP}}^{(\beta_{2})}$, for example) and $A$ is the true fiber output
obtained from the SSFM. The NSD integrates the absolute error squared over the
entire signal duration, and normalizes it with the energy of the signal $A$.
Therefore, the lower the NSD, the more accurate is the waveform predicted by
the model. The inputs for the NSD calculation are illustrated in Fig. (2).

TABLE II: Fiber parameters for C- and O-band transmission Parameter | C-band | O-band
---|---|---
Wavelength $\lambda$ [nm] | 1550 | 1310
$\alpha$ [dB/km] | $0.2$ | $0.4$
$\beta_{2}$ [ps2/km] | $-21.67$ | $-0.2$
$\gamma$ [$1$/W/km] | $1.2$ | $1.4$
$\beta_{3}$ [ps3/km] | $-$ | $0.0765$

Calculating the waveforms for both LP and FLP leads to a numerical issue
related to the ratio in (23) and (29). When the denominator in one of those
two equations tends to zero, the respective model becomes inaccurate. To
address this problem in the LP case, [6] proposed to replace
$A_{\text{LP}}^{(\theta)}(t,z)$ by $A_{\text{RP}}^{(\theta)}(t,z)$ whenever
$|A_{0}(t,z)|<\epsilon$, where $\epsilon>0$ is a fixed threshold. For the FLP
case, we use the same technique in the frequency domain, replacing
$\tilde{A}_{\text{FLP}}^{(\theta)}(\omega,z)$ by
$\tilde{A}_{\text{RP}}^{(\theta)}(\omega,z)$ whenever
$|\tilde{A}_{0}^{(\theta)}(\omega,z)|<\epsilon$ for a certain real
$\epsilon>0$. In addition, $\tilde{A}_{\text{FLP}}^{(\theta)}(\omega,z)$ is
also replaced with $\tilde{A}_{\text{RP}}^{(\theta)}(\omega,z)$ whenever
$|\tilde{A}_{\text{FLP}}^{(\theta)}(\omega,z)|>c\,|\tilde{A}_{\text{RP}}^{(\theta)}(\omega,z)|$
for a certain fixed real $c>1$. In order to increase the stability of the LP
solution, we determined $c$ similarly to the LP case. For the results
presented in this paper $\epsilon$ and $c$ are heuristically fixed to $c=1.15$
and $\epsilon=\text{max}_{t}\\{|A_{0}^{(\gamma)}(t,z)|\\}/10^{5}$ for LP on
$\gamma$ and $c=1.1$ and
$\epsilon=\text{max}_{\omega}\\{|\tilde{A}_{0}^{(\beta_{2})}(\omega,z)\\}/10^{6}$
for FLP on $\beta_{2}$.

For the symbol-level evaluation, decision regions are optimized according to
each model as described in Sec. III-B. In that section, the accuracy is
quantified in terms of BER (illustrated in Fig. 2), while in Sec. III-C, the
accuracy is quantified in terms of achievable information rate (AIR).

$0$$5$$10$$15$$20$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$4.2$
dB$1.5$ dB$\approx$ $\times 42$Input power [dBm]NSD ($\%$)RP on $\gamma$LP on
$\gamma$FLP on $\gamma$RP on $\beta_{2}$LP on $\beta_{2}$FLP on $\beta_{2}$
Figure 3: NSD for for RP, LP and FLP on $\gamma$ and on $\beta_{2}$ in the
C-band using the system in Fig. 2. The fiber parameters are given in Table II.

###  III-A Waveform Comparison

Fig. 3 shows the NSD for RP, LP, and FLP on $\gamma$ and on $\beta_{2}$. As
depicted in Fig. 3, FLP on $\beta_{2}$ (solid red line with squares) is the
most accurate at powers higher than $7.5$ dBm and NSD below $0.1\%$, while LP
on $\gamma$ (solid blue line) is the most accurate at powers lower than $7.5$
dBm. FLP on $\gamma$ (solid blue line with squares) and LP on $\beta_{2}$
(solid red line) have a slightly better performance than RP on $\gamma$
(dashed blue line) and RP on $\beta_{2}$ (dashed red line), respectively.
However, we do not consider FLP on $\gamma$ and LP on $\beta_{2}$ further in
this paper, since their performance is surpassed by LP on $\gamma$ and FLP on
$\beta_{2}$, respectively. RP on $\beta_{2}$ crosses the 0.1$\%$ NSD line at
an input power approximately 4.2 dB higher than RP on $\gamma$, at 14 and 9.8
dBm, respectively. This gap in favor of RP on $\beta_{2}$ was expected since
input powers greater than 10 dBm and small distances such as 20 km put the
fiber in the highly nonlinear regime with low accumulated dispersion. As
discussed in [7], RP on $\beta_{2}$ is accurate on this regime, while RP on
$\gamma$ loses accuracy at high powers. If we change from RP on $\gamma$ to LP
on $\gamma$, the latter outperforms RP on $\beta_{2}$ for powers below 16 dBm.
This gain in accuracy by changing from the LP on $\gamma$ to the RP on
$\gamma$ was previously shown in [6]. In addition, LP on $\gamma$ has its
performance increased due to low accumulated dispersion. By letting
$\beta_{2}\to 0$ in (25), LP on $\gamma$ tends to the NLPN solution, which is
accurate in very low dispersion scenarios [7].

$0.01$$0.1$$1$$10$$10^{-15}$$10^{-12}$$10^{-9}$$10^{-6}$$10^{-3}$$10^{0}$C-band$\beta_{2}$
forO-band$\approx$ $\times 42$$|\beta_{2}|$ [ps2/km]NSD ($\%$)RP on $\gamma$LP
on $\gamma$RP on $\beta_{2}$FLP on $\beta_{2}$ Figure 4: NSD versus
$|\beta_{2}|$ (negative $\beta_{2}$) for four models at an input power of 10
dBm. The system is represented in Fig. 2. All considered models except RP on
$\gamma$ get higher accuracy as $|\beta_{2}|$ decreases.

The dependence of the models on $|\beta_{2}|$ can be seen in Fig. 4, where
four models are compared at a fixed power of 10 dBm for different values of
$|\beta_{2}|$ and the other parameters for C-band transmission (with no TOD).
All the simulated $\beta_{2}$ values were negative. Among the four models, RP
on $\gamma$ is the only one that is virtually invariant to changes in
$\beta_{2}$ for the displayed values. Nevertheless, RP on $\gamma$ has worse
accuracy than LP on $\gamma$ for all displayed values of $\beta_{2}$. When
increasing $|\beta_{2}|$, LP on $\gamma$ increases its NSD at a rate of
approximately $10^{2}$ per decade. Although LP on $\gamma$ outperforms RP on
$\beta_{2}$ at 10 dBm for the C-band, the latter has an increasing rate of
approximately $10^{4}$ per decade, and surpasses the accuracy of LP on
$\gamma$ for $\beta_{2}$ values lower than $-6$ ps2/km. If we now also
consider FLP on $\beta_{2}$, we can gain approximately 42 times in NSD
accuracy with respect to RP on $\beta_{2}$ at 10 dBm. This gap can also be
seen in Fig. 3 and remains approximately constant for different values of
$\beta_{2}$, since both RP and FLP on $\beta_{2}$ have the same increasing
rate of approximately $10^{4}$ per decade. The higher increasing rate for RP
and FLP on $\beta_{2}$ when compared with LP on $\gamma$ shows that the two
former models converge to the true solution of (II) faster than the latter
model when decreasing the accumulated dispersion. For $\beta_{2}=-21.67$
(C-band), FLP on $\beta_{2}$ already outperforms LP on $\gamma$. In addition,
we see back in Fig. 3 that FLP on $\beta_{2}$ crosses the line for an NSD of
0.1$\%$ at an input power 1.5 dB higher than LP on $\gamma$.

$0$$5$$10$$15$$20$$10^{-13}$$10^{-12}$$10^{-11}$$10^{-10}$$10^{-9}$$10^{-8}$$10^{-7}$$10^{-6}$$\beta_{3}=0$$\beta_{3}=0.0765$$\approx$
$\times 91$$14.7$ dBm$18.2$ dBmInput Power [dBm]NSD ($\%$)W/ TODRP on
$\beta_{2}$LP on $\gamma$FLP on $\beta_{2}$W/o TODRP on $\beta_{2}$LP on
$\gamma$FLP on $\beta_{2}$ Figure 5: NSD for the four models in the O-band for
the system in Fig. 2. Results without TOD are also shown.

The results for LP on $\gamma$, RP on $\beta_{2}$, and FLP on $\beta_{2}$
using the O-band parameters in Table II are shown in Fig. 5, where the NSD is
displayed as a function of the input power. First, we consider a system
without TOD (dotted and dashed dotted lines). As shown in Fig. 5, the NSD for
LP on $\gamma$ and FLP on $\beta_{2}$ significantly decay when reducing the
input power in the absence of TOD. On the other hand, RP on $\beta_{2}$
converges to an NSD of approximately $4.6\cdot 10^{-11}\%$ for powers lower
than 2 dBm. This convergence to a non-zero NSD value reflects the mismatch
between (8) and (9) in the absence of nonlinearities. FLP on $\beta_{2}$
outperforms LP on $\gamma$ for all the displayed input powers, while RP on
$\beta_{2}$ outperforms LP on $\gamma$ for input powers higher than 2.5 dBm.
At 10 dBm, the difference in NSD between RP and FLP on $\beta_{2}$ is
approximately 91 times, which is different from the 42 factor for the C-band
results at the same input power (see Fig. 4). This discrepancy is due to the
new set of $\gamma$ and $\alpha$ values, which boost the difference between
the two models at that input power.

The results for the O-band in Fig. 5 shows that when TOD is considered, the
three models now converge to a constant NSD of $10^{-9}$ for input powers
lower than 3 dBm. This behavior can be explained by the absence of TOD in the
model derivations. The error introduced by not accounting for TOD becomes
approximately constant when input powers are lower than 3 dBm, and dominates
the error introduced by incorrectly modeling the other fiber effects. From 0
to 20 dBm, the performance of RP and FLP on $\beta_{2}$ in the system with TOD
is worse than without TOD. This behavior is expected, since for the system
without TOD, RP and FLP on $\beta_{2}$ cross the constant NSD of $10^{-9}$
(TOD error floor for low powers) only at 14.7 dBm and 18.2 dBm, respectively.
For LP on $\gamma$, the NSD for the system without TOD already crosses the
$10^{-9}$ line at 6 dBm. Therefore, results for the systems with and without
TOD converge for powers greater than 10 dBm.

###  III-B Decision Region Optimization

As discussed in [7], comparing models in discrete time can lead to slightly
different conclusions than on a waveform level. For this reason, this section
compares LP on $\gamma$, FLP on $\beta_{2}$, and SSFM results at the symbol
level, measured by BER. The results are shown for the C-band system with
parameters given in Table II. To obtain the symbols, the output waveform from
these three models is filtered by a matched filter and sampled as done for the
signal $A(\cdot,z)$ in Fig. 2. The resulting complex samples are used to
optimize decision regions for each model, originating a symbol detector.
Finally, SSFM simulations are performed to validate the accuracy of each
symbol detector when receiving the true (SSFM) output waveform. These results
are also compared with the PW detector in [12].

$-1$$0$$1$$-1$$0$$1$$\Re(y)$$\Im(y)$$P_{\text{Tx}}=17$
dBm$-1$$0$$1$$-1$$0$$1$$\Re(y)$$\Im(y)$$P_{\text{Tx}}=16$
dBm$-1$$0$$1$$-1$$0$$1$$\Re(y)$$\Im(y)$$P_{\text{Tx}}=15$ dBm (a) SSFM.

$-1$$0$$1$$-1$$0$$1$$\Re(y)$$\Im(y)$$P_{\text{Tx}}=17$
dBm$-1$$0$$1$$-1$$0$$1$$\Re(y)$$\Im(y)$$P_{\text{Tx}}=16$
dBm$-1$$0$$1$$-1$$0$$1$$\Re(y)$$\Im(y)$$P_{\text{Tx}}=15$ dBm (b) FLP on
$\beta_{2}$.

$-1$$0$$1$$-1$$0$$1$$\Re(y)$$\Im(y)$$P_{\text{Tx}}=17$
dBm$-1$$0$$1$$-1$$0$$1$$\Re(y)$$\Im(y)$$P_{\text{Tx}}=16$
dBm$-1$$0$$1$$-1$$0$$1$$\Re(y)$$\Im(y)$$P_{\text{Tx}}=15$ dBm (c) LP on
$\gamma$.

Figure 6: Optimized decision regions obtained with HB detection after the
propagation for $P_{\text{Tx}}\in\\{15,16,17\\}$ dBm of a QPSK constellation
based on: (a) SSFM; (b) FLP on $\beta_{2}$; (c) LP on $\gamma$. Each of the
four differently colored areas represents the decision region for a specific
constellation point. The respective constellation points are illustrated with
red crosses. A contour plot of the histogram for the received samples for the
transmitted symbol $(1-j)/\sqrt{2}$ is plotted on top of each subfigure.

The decision region optimization for the symbol detector is based on the
histogram of the received complex samples. This histogram-based (HB) detector
follows the principle of choosing the most probable transmitted symbol, given
that the corresponding received sample falls in certain area of the complex
plane. This principle corresponds to the maximum a-posteriori (MAP) rule

$\hat{x}=\operatorname*{\text{argmax}}_{s_{m}\in\mathcal{S}}\text{Pr}\\{X=s_{m}\
|\ Y=y\\}.$ (34)

where $X$ and $Y$ are random variables associated with the transmitted symbols
and received samples (see Fig. 2), respectively, and $s_{m}\in\mathcal{S}$
where $\mathcal{S}=\\{s_{1},\cdots,s_{M}\\}$ is the set of constellation
points with cardinality $M$ and $1\leq m\leq M$. The MAP rule, as stated in
(34), is optimal for a memoryless channel, and thus, suboptimal for the
optical fiber channel, which includes memory.

To build the HB detector and numerically approximate the rule in (34),
training symbols are transmitted through the fiber in order to obtain an
estimation of the probability distribution. Their respective received samples
will fall in a specific bin, which is an small area in the complex plane. For
each bin, we count the number $n_{m,b}$ of received samples that fall inside
the $b$-th bin and were drawn from the $m$-th constellation point, where
$1\leq b\leq B$, $1\leq m\leq M$, and $B$ is the total number of bins. The
value of $m=\hat{m}_{b}$ that maximizes $n_{m,b}$ is considered to be the most
probable transmitted constellation point for that bin. If
$\max_{m}\\{n_{m,b}\\}=0$, we replace $b$ by the closest bin $b^{\prime}$ in
Euclidean distance such that $\max_{m}\\{n_{m,b^{\prime}}\\}>0$, and decide
$\hat{m}_{b}=\text{argmax}_{m}\\{n_{m,b^{\prime}}\\}$. After obtaining
$\hat{m}_{b}$ for all possible $b$, the decision regions are defined and every
time a received sample $y$ is received in the $b$-th bin, we assume that
received symbol $\hat{x}=s_{\hat{m}_{b}}$ was transmitted. For computational
reasons, we only consider a subset
$\mathcal{A}=\\{z\in\mathbb{C}\,:\,|\Re\\{z\\}|<2,|\Im\\{z\\}|<2\\}$ of the
complex plane, divided into $B=400\times 400$ square bins of side 0.01. This
region is sufficient to contain virtually all the received samples in the
simulation. We simulate approximately $5.4\cdot 10^{8}$ symbols to obtain
histograms that define $n_{m,b}$. The decision regions were optimized for each
transmitted power separately.

Fig. 6 depicts the decision regions obtained by using HB detectors for SSFM,
FLP on $\beta_{2}$, and LP on $\gamma$ at input powers
$P_{\text{Tx}}\in\\{15,16,17\\}$ dBm. The transmitted constellation is
illustrated with red crosses. Fig. 6 also includes contour plots of the
histogram for the received samples when transmitting the constellation point
$(1-j)/\sqrt{2}$. The decision regions associated to this constellation point
are represented in yellow and contain most of the received samples shown by
the contour plots. The samples that fall outside the yellow regions are not
classified as $(1-j)/\sqrt{2}$ since there are more received samples
originating from another constellation point in those specific bins.

$14$$15$$16$$17$$18$$19$$20$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$\approx$
$\times 17$$\approx$ $\times 5.4$$\approx$ $\times 1.5$FEC
limit$\text{RS}(255,239)$FEC limit$\text{RS}(255,223)$$16.5$ dBmInput power
[dBm]BERMin. dist. HBSSFM HBFLP on $\beta_{2}$ HBLP on $\gamma$ HBPW detec.
$T=2^{11}$PW detec. $T=2^{10}$$\approx$ $\times 1.4$ Figure 7: BER versus
input power for receivers using HB detectors and PW detectors. The HB dectors
were obtained via SSFM, LP on $\gamma$, and FLP on $\beta_{2}$. In all cases,
the channel was simulated in the C-band using the SSFM. Min. dist.: Minimum
distance.

The optimum decision regions are considered as the ones obtained via SSFM in
Fig. LABEL:sub@fig:constssfm. At 15 dBm, the decision regions differ
significantly from the four quadrants of the complex plane that represent the
minimum distance decision regions. This difference originates from the
nonlinearities, which creates a non-Gaussian distribution of the received
samples. Due to the high nonlinear effect present at 17 dBm, the SSFM decision
regions and the contour plots follow a spiral shape. The decision regions
obtained by FLP on $\beta_{2}$ are shown in Fig. LABEL:sub@fig:constb2. At 15
dBm, these decision regions are similar to the SSFM ones for the same
$P_{\text{Tx}}$. However, at 17 dBm, these regions do not present a spiral
shape, as in Fig. LABEL:sub@fig:constssfm. This new behavior can be explained
by the contour plots. The curvature of the contour plots in Fig.
LABEL:sub@fig:constb2 is lower than the curvature in Fig.
LABEL:sub@fig:constssfm, which results in the symbols from the tail of the
contour plot falling in bins with a high number of points transmitted from
another constellation point. The curvature for the contour plot of symbols
obtained via LP on $\gamma$ in Fig. LABEL:sub@fig:constgamma is even lower
than for FLP on $\beta_{2}$ in Fig. LABEL:sub@fig:constb2. This fact results
in decision regions with less accentuated curves in center of Fig.
LABEL:sub@fig:constgamma for each $P_{\text{Tx}}$. At $15$ dBm, the LP on
$\gamma$ decision regions already differ significantly from the SSFM ones. For
example, at both $15$ and $16$ dBm, a small lobe (red dashed circle) present
in both SSFM and FLP on $\beta_{2}$ decision regions is not present in the LP
on $\gamma$ ones. The shape of the contour plots indicate that the decision
regions obtained with FLP on $\beta_{2}$ might perform closer to the SSFM
decision regions than the ones obtained with LP on $\gamma$, as will be
discussed next.

The HB detectors obtained using LP on $\gamma$, FLP on $\beta_{2}$, and SSFM
are compared in a system whose fiber propagation is modeled by the SSFM. The
results are shown in Fig. 7, where the BER is evaluated for different launch
powers using the obtained HB detectors and PW detectors. The latter are
obtained using the SSFM for fiber propagation. As shown in Fig. 7, the SSFM HB
detector (black curve) shows the lowest BER for the displayed input powers.
The minimum distance HB detector gives the worst performance since it assumes
a Gaussian distribution of the received samples. This detector is obtained by
minimum Euclidean distance from the received samples to the possible
transmitted symbols. Replacing the minimum distance HB detector by the SSFM
one at 16 dBm reduces the BER approximately 17 times (from $6.0\cdot 10^{-3}$
to $3.5\cdot 10^{-4}$). At 15 dBm, the BER for FLP on $\beta_{2}$ is 5.4 times
lower than for LP on $\gamma$. As expected, results for FLP on $\beta_{2}$ are
closer to the SSFM results than LP on $\gamma$. The highest gap between the
SSFM and FLP on $\beta_{2}$ occurs at 17 dBm, where the BER for FLP on
$\beta_{2}$ is 1.4 times higher than the BER for SSFM.

For the PW detector [12], instead of defining bins in the complex plane, a
certain amount ($T$) of the $N$ received samples is used as training symbols,
denoted by $y_{t}$, $1\leq t\leq T$. The remaining received samples are used
for testing, denoted by $y_{k}$, where $T+1\leq k\leq N$. The received symbol
$\hat{x}$ depends on the distance between $y_{k}$ and all the training symbols
$y_{t}$ in the following way: for every received sample $y_{k}$, we identify a
set of nearby training samples by calculating the set of time indices
$\mathcal{T}_{s}=\\{t=1,\ldots,T:x_{t}=s,|y_{t}-y_{k}|\leq R\\}$ for every
$s\in\mathcal{S}$ and for a given radius $R>0$. For each transmitted symbol
$s$, $\sum_{t\in\mathcal{T}_{s}}(1/|y_{t}-y_{k}|)$ is calculated, and the
symbol with the highest sum is taken as the received symbol $\hat{x}_{k}$. In
our simulations, we considered two scenarios: $T=2^{10}$ and $T=2^{11}$, both
with $N=2^{16}$. Each scenario was repeated 256 times and the resulting BER
averaged over the repetitions. The radius $R$ was optimized for each
transmitted sequence with a grid search, as illustrated in [12, Fig. 3].

As depicted in Fig. 7, the PW detector BER for $2^{11}$ training symbols
converges to the BER of FLP on $\beta_{2}$ at input powers higher than 17 dBm.
At input powers lower than 16 dBm, the PW BER for $2^{11}$ training symbols is
even higher than the BER for LP on $\gamma$. This result shows that the HB
detector can perform close to recently proposed methods in the literature. At
15 dBm, the penalty obtained for reducing the number of training symbols from
$2^{11}$ to $2^{10}$ is an increase in BER of approximately 1.5 times. The PW
detector performance could be further improved by increasing $T$, at the cost
of increased complexity at the receiver.

###  III-C Achievable Information Rates

FEC is present in modern PON systems to improve system performance [17]. In
Fig. 7, pre-FEC BER thresholds are shown for two Reed-Solomon (RS) codes [18].
The considered codes in Fig. 7 are $\text{RS}(n,k)$ with $k=239,223$ and
$n=255$, where $k$ and $n$ are the information and codeword lengths, resp.
These two codes are typical low-complexity RS codes used in PONs [19, 17] and
have a code rate of $0.93$ and $0.87$, respectively.

Along with RS codes, also stronger FEC codes such as low-density parity-check
codes or staircase codes have been proposed for PONs in the literature [20,
21, 22]. In this section, we evaluate the models for PON systems with hard-
decision (HD) FEC in terms of achievable information rates (AIRs) [23]. We
consider the simple RS codes described above and a theoretical HD limit for
the AIR. The latter is close to the performance of strong HD FEC codes [24,
Fig. 8]. The AIRs are obtained by closed-form expressions based on pre-FEC
BER.

Bits toRS Symb.RSEnc.Inter-leaverRS Symb.to Bits$010$FiberChannelRS Symb.to
BitsRSDec.Bits toRS Symb.Deinter-leaver$110$$110$$010$pre-FEC BER $p$post-FEC
BER $p_{\text{pos}}$BSC Figure 8: Encoding/decoding procedures used in this
paper for RS codes. The fiber channel block include all components from Fig.
2. The post-FEC BER $p_{\text{pos}}$ is estimated from the pre-FEC BER $p$
using (35). Enc.: encoder; Dec.: decoder; BSC: binary symmetric channel.
$15$$15.5$$16$$16.5$$17$$17.5$$18$$18.5$$19$$19.5$$20$$1.2$$1.4$$1.6$$1.8$$2$RS
system $p_{\text{pos}}\\!\approx\\!10^{-12}$Theoretical
bounds$\text{AIR}_{\text{RS}}$ for$\text{RS}(255,239)$$\text{AIR}_{\text{RS}}$
for$\text{RS}(255,223)$$1.8$ dB$0.4$ dB$1.9$ dBMin. dist.
$p_{\text{pos}}\\!\approx\\!10^{-5}$SSFM
$p_{\text{pos}}\\!\approx\\!10^{-5}$$\text{RS}(255,177)$$\text{RS}(255,199)$Input
power [dBm]AIR [bits/symbol]Min. dist. HB - RSMin. dist. HB - Th.SSFM HB -
RSSSFM HB - Th.FLP on $\beta_{2}$ HB - RSFLP on $\beta_{2}$ HB - Th.LP on
$\gamma$ HB - RSLP on $\gamma$ HB - Th. Figure 9: Comparison of AIRs when
using the decision regions obtained in Sec. III-B. The AIRs were obtained
using (36) and (35) at $p_{\text{pos}}\approx 10^{-12}$ for the RS system and
(37) for the theoretical bound. The dotted lines represent the comparison
between $\text{AIR}_{\text{RS}}$ when approximating $p_{\text{pos}}$ by (35)
(white triangles) and $\text{AIR}_{\text{RS}}$ when simulating the system in
Fig. 8 (filled diamonds), both for $p_{\text{pos}}\approx 10^{-5}$. Th.:
theoretical.

We consider a family of RS codes $\text{RS}(n,k)$ with multiple coding rates,
where $k$ is varied to obtain different code rates. We use a fixed codeword
length ($n=255$ symbols) in order to constrain the code complexity. For every
launch power, $k$ is determined by finding the highest $k$ such that the post-
FEC BER falls below a certain threshold. Following [13, Table IV.2], we use
$10^{-12}$ as post-FEC BER threshold.

The post-FEC BER $p_{\text{pos}}$ can be approximated by substituting the pre-
FEC BER $p$ from the system in Fig. 2 in the analytical expression for binary
symmetric channels (BSCs) and bounded-distance decoders [25]

$p_{\text{pos}}\approx\frac{1}{n}\sum_{r=t+1}^{n}\left(\dfrac{p}{p_{s}}r+\dfrac{1}{2(t-1)!}\right){n\choose
r}p_{s}^{r}(1-p_{s})^{n-r},$ (35)

where $t=\lfloor(n-k)/2\rfloor$ is the RS error-correction capability,
$p_{s}=1-(1-p)^{m}$ is the (RS) symbol error probability, and
$m=\lceil\log_{2}(n+1)\rceil$ is the number of bits per symbol. To improve the
total computation time, a binary search algorithm on $k$ over all integers
between 1 and 253 is performed. After finding $k$, the AIR for the RS system
($\text{AIR}_{\text{RS}}$) is determined by

$\text{AIR}_{\text{RS}}=\log_{2}(M)\dfrac{k}{n}=\dfrac{2}{255}k,$ (36)

since $M=4$ for QPSK. We call (36) an achievable information rate since we
assume that $p_{\text{pos}}<10^{-12}$ can be considered virtually error-free
for the system considered in this paper.

Since the expression in (35) is valid for BSCs, we need to modify the system
in Fig. 2 to fulfill that property. The fiber channel in Fig. 2 presents
memory due to the interaction of dispersion and nonlinearities. Therefore, in
our simulations we included a bit interleaver and a bit deinterleaver so that
the fiber channel in Fig. 2 is well-approximated by a BSC. The resulting
system, together with the RS encoding and decoding blocks, can be seen in Fig.
8. The system in Fig. 8 was only simulated for $p_{\text{pos}}\approx 10^{-5}$
and was used to validate the results of (35). After this validation, (35) was
used instead of simulating the system in Fig. 8.

Fig. 9 depicts the AIRs for the RS system using (35). Before analyzing the
results for $p_{\text{pos}}\approx 10^{-12}$, we validate (35) at
$p_{\text{pos}}\approx 10^{-5}$ to determine if the system in Fig. 8 can be
well-approximated by a BSC. This validation is done assuming a minimum
distance detector and is shown by the dotted lines in Fig. 9. The results for
the RS simulations of the exact system in Fig. 8 at $p_{\text{pos}}\approx
10^{-5}$ are shown with diamonds. An almost perfect overlap with the results
from (35) with the same $p_{\text{pos}}$ (shown with triangles) is observed.
The agreement between these curves suggest that the system in Fig. 8 can be
approximated by a BSC. Therefore, from now on only (35) is used for the RS
systems at $p_{\text{pos}}\approx 10^{-12}$.

We start by comparing the dashed lines in Fig. 9, which represent the
$\text{AIR}_{\text{RS}}$ results in (36). As shown in Fig. 9, using the
decision regions obtained by SSFM in a RS system can provide a gain of
approximately $1.8$ dB for a rate of $1.6$ bits/symbol in terms of nonlinear
tolerance. The crossing point between the rate of $\text{RS}(255,223)$ with
$\text{AIR}_{\text{RS}}$ for the SSFM decision regions is at $16.5$ dBm, which
closely matches with the crossing point with BER in Fig. 7. When comparing the
models, the histogram-based detector obtained using FLP on $\beta_{2}$
outperforms the LP on $\gamma$ one throughout the considered power range,
analogously to Fig. 7. Specifically for $\text{RS}(255,239)$, FLP on
$\beta_{2}$ outperforms LP on $\gamma$ by approximately $0.4$ dB.

The $\text{AIR}_{\text{RS}}$ is also compared with a theoretical bound on
hard-decision bit-wise AIRs for independent, identically distributed bit
errors. The theoretical AIR ($\text{AIR}_{\text{TH}}$) used in this paper is
defined as [26]

$\text{AIR}_{\text{TH}}=\log_{2}(M)\left(1-\text{H}_{\text{b}}(p)\right),$
(37)

where $\text{H}_{\text{b}}(p)=-p\log_{2}(p)-(1-p)\log_{2}(1-p)$ is the binary
entropy function for the given pre-FEC BER $p$ in the system of Fig. 2. The
$\text{AIR}_{\text{TH}}$ from (37) can be approached by strong FEC codes such
as staircase codes, as reported in [24, Fig. 8].

As shown in Fig. 9, the theoretical bounds (solid lines) from (37) show
significant gains over the $\text{AIR}_{\text{RS}}$ results. These gains show
that, by using codes more complex than $\text{RS}(255,k)$, higher code rates
can be achieved or the input power can be improved for a specific rate. For
the same code rate as in $\text{RS}(255,223)$, the input power for the
theoretical bound on the SSFM decision regions is approximately $1.9$ dB
higher than the one for the RS system. However, the complexity and latency of
codes that perform close to the theoretical bound should be carefully analysed
for a PON system design. The results for the theoretical bounds in in Fig. 9
also show that the histogram-based detector based on the FLP on $\beta_{2}$
outperforms the one based on LP on $\gamma$.

##  IV Conclusions

In this paper we presented a novel model for optical fiber transmission and
evaluated its performance for a passive optical network system. The proposed
model was derived as an improved version of the regular perturbation on the
GVD parameter model. The improvement was obtained by applying frequency
logarithmic perturbation on the GVD parameter. Both regular and frequency
logarithmic perturbation on the GVD parameter models are suitable in the
weakly dispersive and highly nonlinear regime, whereas the frequency
logarithmic perturbation is able to surpass the limitations of the regular
perturbation.

Apart from the regular perturbation on the GVD parameter, the proposed model
was compared with two other models present in the literature: regular and
logarithmic perturbation on the Kerr nonlinear coefficient. For a fixed
normalized squared deviation of $0.1\%$, the proposed model was accurate at
$1.5$ dB higher input powers compared to logarithmic perturbation on the Kerr
nonlinear coefficient. Both frequency logarithmic and regular perturbation on
the GVD parameter exhibit the highest convergence rate to the split-step
Fourier method results when reducing the dispersion effect. The proposed model
also proved to be more suitable for symbol and bit detection, with and without
FEC.

Possible extensions of this work are higher-order logarithmic perturbation
models, perturbation on the GVD parameter for dual-polarization systems, and
nonlinearity-compensation techniques based on the proposed model.

## Acknowledgments

The authors would like to thank Dr. Tobias Fehenberger (ADVA Optical
Networking) for fruitful discussions about the regular perturbation on
$\gamma$. We would also like to thank Dr. Nicola Calabretta (Eindhoven
University of Technology) and Dr. Domaniç Lavery (Infinera) for insightful
discussions about passive optical networks.

## References

  * [1] G. Agrawal, _Nonlinear Fiber Optics_ , 5th ed., ser. Optics and Photonics. Boston: Academic Press, 2013\.
  * [2] O. V. Sinkin, R. Holzlohner, J. Zweck, and C. R. Menyuk, “Optimization of the split-step Fourier method in modeling optical-fiber communications systems,” _J. Lightwave Technol._ , vol. 21, no. 1, pp. 61–68, Jan. 2003\.
  * [3] L. Beygi, E. Agrell, M. Karlsson, and P. Johannisson, “Signal statistics in fiber-optical channels with polarization multiplexing and self-phase modulation,” _J. Lightwave Technol._ , vol. 29, no. 16, pp. 2379–2386, Aug. 2011.
  * [4] E. Forestieri and M. Secondini, “Solving the nonlinear Schrödinger equation,” in _Optical Communication Theory and Techniques_ , 1st ed., E. Forestieri, Ed. Boston: Springer US, 2005, pp. 3–11.
  * [5] A. Vannucci, P. Serena, and A. Bononi, “The RP method: a new tool for the iterative solution of the nonlinear Schödinger equation,” _J. Lightwave Technol._ , vol. 20, no. 7, pp. 1102–1112, July 2002.
  * [6] E. Ciaramella and E. Forestieri, “Analytical approximation of nonlinear distortions,” _IEEE Photonics Technol. Lett._ , vol. 17, no. 1, pp. 91–93, Jan. 2005.
  * [7] V. Oliari, E. Agrell, and A. Alvarado, “Regular perturbation on the group-velocity dispersion parameter for nonlinear fibre-optical communications,” _Nature Communications_ , vol. 11, no. 933, pp. 1–11, Feb. 2020.
  * [8] V. Oliari, E. Agrell, and A. Alvarado, “Logarithmic perturbation models in the weak-dispersion regime with applications to passive optical networks,” in _Proc. European Conf. Optical Communication (ECOC)_ , Brussels, Belgium, Dec. 2020.
  * [9] M. Secondini and E. Forestieri, “Analytical fiber-optic channel model in the presence of cross-phase modulation,” _IEEE Photonics Technol. Lett._ , vol. 24, no. 22, pp. 2016–2019, Sep. 2012.
  * [10] M. Secondini, E. Forestieri, and G. Prati, “Achievable information rate in nonlinear wdm fiber-optic systems with arbitrary modulation formats and dispersion maps,” _J. Lightwave Technol._ , vol. 31, no. 23, pp. 3839–3852, Nov. 2013.
  * [11] M. Secondini, E. Agrell, E. Forestieri, D. Marsella, and M. R. Camara, “Nonlinearity mitigation in WDM systems: models, strategies, and achievable rates,” _J. Lightwave Technol._ , vol. 37, no. 10, pp. 2270–2283, Feb. 2019.
  * [12] A. Amari, X. Lin, O. Dobre, R. Venkatesan, and A. Alvarado, “A machine learning-based detection technique for optical fiber nonlinearity mitigation,” _IEEE Photonics Technol. Lett._ , vol. 31, no. 8, pp. 627–630, Mar. 2019.
  * [13] _10-Gigabit-capable passive optical networks (XG-PON): physical media dependent (PMD) layer specification_ , ITU-T G.987.2, 2016.
  * [14] _PON transmission technologies above 10 Gb/s per wavelength_ , ITU-T G.Sup64, 2018.
  * [15] D. Lavery, M. Paskov, R. Maher, B. C. Thomsen, S. J. Savory, and P. Bayvel, “Low complexity multichannel nonlinear predistortion for passive optical networks,” in _Advanced Photonics 2015_. Optical Society of America, Jun. 2015, paper SpS2C.5.
  * [16] D. Lavery, R. Maher, D. Millar, A. Alvarado, S. J. Savory, and P. Bayvel, “Why compensating fibre nonlinearity will never meet capacity demands,” _arXiv:1512.03426_ , Dec. 2019.
  * [17] _10-Gigabit-capable passive optical networks (XG-PON): transmission convergence (TC) layer specification_ , ITU-T G.987.3, 2014.
  * [18] I. Reed and G. Solomon, “Polynomial codes over certain finite fields,” _J. Soc. Indust. Appl. Math._ , vol. 8, no. 2, pp. 300–304, Jun. 1960.
  * [19] L. Schmalen, A. J. de Lind van Wijngaarden, and S. Ten Brink, “Forward error correction in optical core and optical access networks,” _Bell Labs Technical Journal_ , vol. 18, no. 3, pp. 39–66, 2013.
  * [20] B. Powell and K. Droskiewicz, “Latency & complexity for various 25/50/100G EPON FEC code proposals,” Sep. 2017 (accessed on August 20, 2020). [Online]. Available: https://www.ieee802.org/3/ca/public/meeting˙archive/2017/09/powell˙3ca˙1a˙0917.pdf.
  * [21] D. van Veen, V. Houtsma, A. de Lind van Wijngaarden, B. Powell, and E. Harstead, “FEC code for 25/50/100G EPON,” Mar. 2017 (accessed on August 20, 2020). [Online]. Available: https://www.ieee802.org/3/ca/email/pdfUGjlG0uq7U.pdf.
  * [22] A. Teixeira, _et al._ , “DSP enabled optical detection techniques for PON,” _J. Lightwave Technol._ , vol. 38, no. 3, pp. 684–695, Feb. 2020\.
  * [23] A. Alvarado, T. Fehenberger, B. Chen, and F. M. J. Willems, “Achievable information rates for fiber optics: applications and computations,” _J. Lightwave Technol._ , vol. 36, no. 2, pp. 424–439, Jan. 2018.
  * [24] B. P. Smith, A. Farhood, A. Hunt, F. R. Kschischang, and J. Lodge, “Staircase codes: FEC for 100 Gb/s OTN,” _J. Lightwave Technol._ , vol. 30, no. 1, pp. 110–117, Jan. 2012.
  * [25] _Forward error correction for high bit-rate DWDM submarine systems_ , ITU-T G.975.1, 2004.
  * [26] T. Fehenberger, A. Alvarado, P. Bayvel, and N. Hanik, “On achievable rates for long-haul fiber-optic communications,” _Opt. Express_ , vol. 23, no. 7, pp. 9183–9191, Apr. 2015.


2 | 1808 |  | 15796 |  | 91175 |  | 310203 |  | 1942 |  |  |
3 |  | 135 |  | 1943 |  | 15797 |  | 75779 |  | 135 |  |
4 |  |  | 1 |  | 135 |  | 1943 |  | 13856 |  | 1 |
5 |  |  |  |  |  | 1 |  | 135 |  | 1808 |  |
6 |  |  |  |  |  |  |  |  | 1 |  | 134 |
7 |  |  |  |  |  |  |  |  |  |  |  | 1

$d=5$

Table 24: BPS spectrum of the Wilson loop in the fundamental representation
for $E_{7}$ del Pezzo.

$2j_{L}\backslash 2j_{R}$ | 0
---|---
0 | 1
$d$ = $-1$

$2j_{L}\backslash 2j_{R}$ | 0 | 1 | 2
---|---|---|---
0 | 4125 |  |
1 |  | 249 |
2 |  |  | 1
$d$ = $1$

$2j_{L}\backslash 2j_{R}$ | 0 | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---|---
0 |  | 186126 |  | 249 |  |
1 | 4124 |  | 38877 |  | 1 |
2 |  | 249 |  | 4373 |  |
3 |  |  | 1 |  | 249 |
4 |  |  |  |  |  | 1
$d$ = $2$

$2j_{L}\backslash 2j_{R}$ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
---|---|---|---|---|---|---|---|---|---|---
0 | 3694119 |  | 11393622 |  | 252004 |  | 249 |  |  |
1 |  | 1434130 |  | 4880618 |  | 43498 |  | 1 |  |
2 | 39125 |  | 295005 |  | 1286881 |  | 4623 |  |  |
3 |  | 4622 |  | 43747 |  | 256377 |  | 250 |  |
4 | 1 |  | 250 |  | 4623 |  | 39374 |  | 1 |
5 |  |  |  | 1 |  | 250 |  | 4374 |  |
6 |  |  |  |  |  |  | 1 |  | 249 |
7 |  |  |  |  |  |  |  |  |  | 1

$d=3$ $2j_{L}\backslash 2j_{R}$ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 0 485875765
799689237 72308485 1357379 4374 1 72035619 292038369 545982777 28547872 300373
250 2 29699120 106691237 234719994 7613632 44247 1 3 1387759 7873883 30967129
80571370 1649255 4625 4 300621 1692754 7688504 23680624 300623 250 5 4375
44247 305246 1649753 6079371 43998 1 6 250 4625 44248 300624 1353506 4624 7 1
250 4625 43998 256875 250 8 1 250 4624 39375 1 9 1 250 4374 10 1 249 11 1
$d=4$

Table 25: BPS spectrum of the Wilson loop in the fundamental representation
for $E_{8}$ del Pezzo.

## References

  * [1] S. H. Katz, A. Klemm and C. Vafa,  _Geometric engineering of quantum field theories_ ,  Nucl. Phys. B 497, 173 (1997),  10.1016/S0550-3213(97)00282-4,  hep-th/9609239.
  * [2] E. Witten,  _Chern-Simons gauge theory as a string theory_ ,  Prog. Math. 133, 637 (1995),  hep-th/9207094.
  * [3] M. Aganagic, A. Klemm, M. Marino and C. Vafa,  _Matrix model as a mirror of Chern-Simons theory_ ,  JHEP 02, 010 (2004),  10.1088/1126-6708/2004/02/010,  hep-th/0211098.
  * [4] M. Aganagic, R. Dijkgraaf, A. Klemm, M. Marino and C. Vafa,  _Topological strings and integrable hierarchies_ ,  Commun. Math. Phys. 261, 451 (2006),  10.1007/s00220-005-1448-9,  hep-th/0312085.
  * [5] S. Katz, P. Mayr and C. Vafa,  _Mirror symmetry and exact solution of 4-D N=2 gauge theories: 1._ ,  Adv. Theor. Math. Phys. 1, 53 (1998),  10.4310/ATMP.1997.v1.n1.a2,  hep-th/9706110.
  * [6] T. Dimofte, S. Gukov and L. Hollands,  _Vortex Counting and Lagrangian 3-manifolds_ ,  Lett. Math. Phys. 98, 225 (2011),  10.1007/s11005-011-0531-8,  1006.0977.
  * [7] M. Aganagic, A. Klemm and C. Vafa,  _Disk instantons, mirror symmetry and the duality web_ ,  Z. Naturforsch. A 57, 1 (2002),  10.1515/zna-2002-1-201,  hep-th/0105045.
  * [8] M. Aganagic, A. Klemm, M. Marino and C. Vafa,  _The Topological vertex_ ,  Commun. Math. Phys. 254, 425 (2005),  10.1007/s00220-004-1162-z,  hep-th/0305132.
  * [9] A. Iqbal, C. Kozcaz and C. Vafa,  _The Refined topological vertex_ ,  JHEP 10, 069 (2009),  10.1088/1126-6708/2009/10/069,  hep-th/0701156.
  * [10] C. Kozçaz, S. Shakirov, C. Vafa and W. Yan,  _Refined Topological Branes_ ,  Commun. Math. Phys. 385(2), 937 (2021),  10.1007/s00220-020-03883-1,  1805.00993.
  * [11] H.-C. Kim, M. Kim and S.-S. Kim,  _5d/6d Wilson loops from blowups_ ,  JHEP 08, 131 (2021),  10.1007/JHEP08(2021)131,  2106.04731.
  * [12] N. A. Nekrasov,  _Seiberg-Witten prepotential from instanton counting_ ,  Adv. Theor. Math. Phys. 7(5), 831 (2003),  10.4310/ATMP.2003.v7.n5.a4,  hep-th/0206161.
  * [13] N. Nekrasov and A. Okounkov,  _Seiberg-Witten theory and random partitions_ ,  Prog. Math. 244, 525 (2006),  10.1007/0-8176-4467-9_15,  hep-th/0306238.
  * [14] H.-C. Kim, S. Kim, E. Koh, K. Lee and S. Lee,  _On instantons as Kaluza-Klein modes of M5-branes_ ,  JHEP 12, 031 (2011),  10.1007/JHEP12(2011)031,  1110.2175.
  * [15] D. Rodríguez-Gómez and G. Zafrir,  _On the 5d instanton index as a Hilbert series_ ,  Nucl. Phys. B 878, 1 (2014),  10.1016/j.nuclphysb.2013.11.006,  1305.5684.
  * [16] O. Bergman, D. Rodríguez-Gómez and G. Zafrir,  _Discrete $\theta$ and the 5d superconformal index_,  JHEP 01, 079 (2014),  10.1007/JHEP01(2014)079,  1310.2150.
  * [17] C. Hwang, J. Kim, S. Kim and J. Park,  _General instanton counting and 5d SCFT_ ,  JHEP 07, 063 (2015),  10.1007/JHEP07(2015)063,  [Addendum: JHEP 04, 094 (2016)],  1406.6793.
  * [18] R. Gopakumar and C. Vafa,  _M theory and topological strings. 1._ (1998),  hep-th/9809187.
  * [19] R. Gopakumar and C. Vafa,  _M theory and topological strings. 2._ (1998),  hep-th/9812127.
  * [20] T. W. Grimm, A. Klemm, M. Marino and M. Weiss,  _Direct Integration of the Topological String_ ,  JHEP 08, 058 (2007),  10.1088/1126-6708/2007/08/058,  hep-th/0702187.
  * [21] B. Haghighat, A. Klemm and M. Rauch,  _Integrability of the holomorphic anomaly equations_ ,  JHEP 10, 097 (2008),  10.1088/1126-6708/2008/10/097,  0809.1674.
  * [22] M.-x. Huang and A. Klemm,  _Direct integration for general $\Omega$ backgrounds_,  Adv. Theor. Math. Phys. 16(3), 805 (2012),  10.4310/ATMP.2012.v16.n3.a2,  1009.1126.
  * [23] D. Tong and K. Wong,  _Instantons, Wilson lines, and D-branes_ ,  Phys. Rev. D 91(2), 026007 (2015),  10.1103/PhysRevD.91.026007,  1410.8523.
  * [24] N. Nekrasov,  _BPS/CFT correspondence: non-perturbative Dyson-Schwinger equations and qq-characters_ ,  JHEP 03, 181 (2016),  10.1007/JHEP03(2016)181,  1512.05388.
  * [25] H.-C. Kim,  _Line defects and 5d instanton partition functions_ ,  JHEP 03, 199 (2016),  10.1007/JHEP03(2016)199,  1601.06841.
  * [26] T. Kimura, H. Mori and Y. Sugimoto,  _Refined geometric transition and $qq$-characters_,  JHEP 01, 025 (2018),  10.1007/JHEP01(2018)025,  1705.03467.
  * [27] P. Agarwal, J. Kim, S. Kim and A. Sciarappa,  _Wilson surfaces in M5-branes_ ,  JHEP 08, 119 (2018),  10.1007/JHEP08(2018)119,  1804.09932.
  * [28] B. Assel and A. Sciarappa,  _Wilson loops in 5d $\mathcal{N}=1$ theories and S-duality_,  JHEP 10, 082 (2018),  10.1007/JHEP10(2018)082,  1806.09636.
  * [29] N. Haouzi and J. Oh,  _On the Quantization of Seiberg-Witten Geometry_ ,  JHEP 01, 184 (2021),  10.1007/JHEP01(2021)184,  2004.00654.
  * [30] J. Chen, B. Haghighat, H.-C. Kim and M. Sperling,  _Elliptic quantum curves of class ${\mathcal{S}}_{k}$_,  JHEP 03, 028 (2021),  10.1007/JHEP03(2021)028,  2008.05155.
  * [31] J. Chen, B. Haghighat, H.-C. Kim, M. Sperling and X. Wang,  _E-string quantum curve_ ,  Nucl. Phys. B 973, 115602 (2021),  10.1016/j.nuclphysb.2021.115602,  2103.16996.
  * [32] J. Chen, B. Haghighat, H.-C. Kim, K. Lee, M. Sperling and X. Wang,  _Elliptic Quantum Curves of 6d SO(N) theories_ (2021),  2110.13487.
  * [33] S. Nawata, K. Zhang and R.-D. Zhu,  _ABCD of qq-characters_ (2023),  2302.00525.
  * [34] M.-x. Huang, K. Lee and X. Wang,  _Topological strings and Wilson loops_ ,  JHEP 08, 207 (2022),  10.1007/JHEP08(2022)207,  2205.02366.
  * [35] M.-x. Huang, A.-K. Kashani-Poor and A. Klemm,  _The $\Omega$ deformed B-model for rigid $\mathcal{N}=2$ theories_,  Annales Henri Poincare 14, 425 (2013),  10.1007/s00023-012-0192-x,  1109.5728.
  * [36] M.-x. Huang and A. Klemm,  _Holomorphic Anomaly in Gauge Theories and Matrix Models_ ,  JHEP 09, 054 (2007),  10.1088/1126-6708/2007/09/054,  hep-th/0605195.
  * [37] M.-x. Huang, A. Klemm and S. Quackenbush,  _Topological string theory on compact Calabi-Yau: Modularity and boundary conditions_ ,  Lect. Notes Phys. 757, 45 (2009),  10.1007/978-3-540-68030-7_3,  hep-th/0612125.
  * [38] H. Nakajima and K. Yoshioka,  _Instanton counting on blowup. 1._ ,  Invent. Math. 162, 313 (2005),  10.1007/s00222-005-0444-1,  math/0306198.
  * [39] H. Nakajima and K. Yoshioka,  _Instanton counting on blowup. II. K-theoretic partition function_ (2005),  math/0505553.
  * [40] H. Nakajima and K. Yoshioka,  _Perverse coherent sheaves on blowup, III: Blow-up formula from wall-crossing_ ,  Kyoto J. Math. 51(2), 263 (2011),  10.1215/21562261-1214366,  0911.1773.
  * [41] M.-x. Huang, K. Sun and X. Wang,  _Blowup Equations for Refined Topological Strings_ ,  JHEP 10, 196 (2018),  10.1007/JHEP10(2018)196,  1711.09884.
  * [42] A. Klemm,  _The B-model approach to topological string theory on Calabi-Yau n-folds_ ,  In _B-model Gromov-Witten theory_ , Trends Math., pp. 79–397. Birkhäuser/Springer, Cham (2018).
  * [43] H.-C. Kim, M. Kim, S.-S. Kim and K.-H. Lee,  _Bootstrapping BPS spectra of 5d/6d field theories_ ,  JHEP 04, 161 (2021),  10.1007/JHEP04(2021)161,  2101.00023.
  * [44] M. Dedushenko and E. Witten,  _Some Details On The Gopakumar-Vafa and Ooguri-Vafa Formulas_ ,  Adv. Theor. Math. Phys. 20, 1 (2016),  10.4310/ATMP.2016.v20.n1.a1,  1411.7108.
  * [45] D. Krefl and J. Walcher,  _Extended Holomorphic Anomaly in Gauge Theory_ ,  Lett. Math. Phys. 95, 67 (2011),  10.1007/s11005-010-0432-2,  1007.0263.
  * [46] M. Bershadsky, S. Cecotti, H. Ooguri and C. Vafa,  _Kodaira-Spencer theory of gravity and exact results for quantum string amplitudes_ ,  Commun. Math. Phys. 165, 311 (1994),  10.1007/BF02099774,  hep-th/9309140.
  * [47] M.-X. Huang, A. Klemm and M. Poretschkin,  _Refined stable pair invariants for E-, M- and $[p,q]$-strings_,  JHEP 11, 112 (2013),  10.1007/JHEP11(2013)112,  1308.0619.
  * [48] K. A. Intriligator, D. R. Morrison and N. Seiberg,  _Five-dimensional supersymmetric gauge theories and degenerations of Calabi-Yau spaces_ ,  Nucl. Phys. B 497, 56 (1997),  10.1016/S0550-3213(97)00279-4,  hep-th/9702198.
  * [49] M. Aganagic and S. Shakirov,  _Refined Chern-Simons Theory and Topological String_ (2012),  1210.2733.
  * [50] J. Choi, S. Katz and A. Klemm,  _The refined BPS index from stable pair invariants_ ,  Commun. Math. Phys. 328, 903 (2014),  10.1007/s00220-014-1978-0,  1210.4403.
  * [51] A. Klemm, J. Manschot and T. Wotschke,  _Quantum geometry of elliptic Calabi-Yau manifolds_ (2012),  1205.1795.
  * [52] C. F. Cota, A. Klemm and T. Schimannek,  _Topological strings on genus one fibered Calabi-Yau 3-folds and string dualities_ ,  JHEP 11, 170 (2019),  10.1007/JHEP11(2019)170,  1910.01988.
  * [53] S.-S. Kim and X.-Y. Wei,  _Refined topological vertex with ON-planes_ ,  JHEP 08, 006 (2022),  10.1007/JHEP08(2022)006,  2201.12264.
  * [54] T. M. Chiang, A. Klemm, S.-T. Yau and E. Zaslow,  _Local mirror symmetry: Calculations and interpretations_ ,  Adv. Theor. Math. Phys. 3, 495 (1999),  10.4310/ATMP.1999.v3.n3.a3,  hep-th/9903053.
  * [55] _https://github.com/wangxin07895123/Wilsonloop5D_ .
  * [56] M. Aganagic, V. Bouchard and A. Klemm,  _Topological Strings and (Almost) Modular Forms_ ,  Commun. Math. Phys. 277, 771 (2008),  10.1007/s00220-007-0383-3,  hep-th/0607100.
  * [57] M.-x. Huang and A. Klemm,  _Holomorphicity and Modularity in Seiberg-Witten Theories with Matter_ ,  JHEP 07, 083 (2010),  10.1007/JHEP07(2010)083,  0902.1325.
  * [58] M.-x. Huang,  _Modular anomaly from holomorphic anomaly in mass deformed $\mathcal{N}=2$ superconformal field theories_,  Phys. Rev. D 87(10), 105010 (2013),  10.1103/PhysRevD.87.105010,  1302.6095.
  * [59] N. Seiberg and E. Witten,  _Monopoles, duality and chiral symmetry breaking in N=2 supersymmetric QCD_ ,  Nucl. Phys. B 431, 484 (1994),  10.1016/0550-3213(94)90214-3,  hep-th/9408099.
  * [60] N. Seiberg and E. Witten,  _Electric - magnetic duality, monopole condensation, and confinement in N=2 supersymmetric Yang-Mills theory_ ,  Nucl. Phys. B 426, 19 (1994),  10.1016/0550-3213(94)90124-4,  hep-th/9407087.
  * [61] N. A. Nekrasov and S. L. Shatashvili,  _Quantization of Integrable Systems and Four Dimensional Gauge Theories_ ,  In _16th International Congress on Mathematical Physics_ , pp. 265–289,  10.1142/9789814304634_0015 (2009), 0908.4052.
  * [62] A. B. Goncharov and R. Kenyon,  _Dimers and cluster integrable systems_ ,  Ann. Sci. Éc. Norm. Supér. (4) 46(5), 747 (2013),  10.24033/asens.2201,  1107.5588.
  * [63] R. Eager, S. Franco and K. Schaeffer,  _Dimer Models and Integrable Systems_ ,  JHEP 06, 106 (2012),  10.1007/JHEP06(2012)106,  1107.1244.
  * [64] S. Moriyama,  _Spectral Theories and Topological Strings on del Pezzo Geometries_ ,  JHEP 10, 154 (2020),  10.1007/JHEP10(2020)154,  2007.05148.
  * [65] T. Furukawa, S. Moriyama and Y. Sugimoto,  _Quantum Mirror Map for Del Pezzo Geometries_ ,  J. Phys. A 53(38), 38 (2020),  10.1088/1751-8121/ab93fe,  1908.11396.
  * [66] D. Gaiotto and H.-C. Kim,  _Duality walls and defects in 5d $\mathcal{N}=1$ theories_,  JHEP 01, 019 (2017),  10.1007/JHEP01(2017)019,  1506.03871.
  * [67] A. Sciarappa,  _Exact relativistic Toda chain eigenfunctions from Separation of Variables and gauge theory_ ,  JHEP 10, 116 (2017),  10.1007/JHEP10(2017)116,  1706.05142.
  * [68] A. Grassi and M. Marino,  _The complex side of the TS/ST correspondence_ ,  J. Phys. A 52(5), 055402 (2019),  10.1088/1751-8121/aaec4b,  1708.08642.
  * [69] A. Grassi and M. Mariño,  _A Solvable Deformation of Quantum Mechanics_ ,  SIGMA 15, 025 (2019),  10.3842/SIGMA.2019.025,  1806.01407.
  * [70] A. Mironov and A. Morozov,  _Nekrasov Functions and Exact Bohr-Zommerfeld Integrals_ ,  JHEP 04, 040 (2010),  10.1007/JHEP04(2010)040,  0910.5670.
  * [71] M. Aganagic, M. C. N. Cheng, R. Dijkgraaf, D. Krefl and C. Vafa,  _Quantum Geometry of Refined Topological Strings_ ,  JHEP 11, 019 (2012),  10.1007/JHEP11(2012)019,  1105.0630.
  * [72] M.-x. Huang and X.-f. Wang,  _Topological Strings and Quantum Spectral Problems_ ,  JHEP 09, 150 (2014),  10.1007/JHEP09(2014)150,  1406.6178.
  * [73] J. D. Edelstein, M. Gomez-Reino and M. Marino,  _Blowup formulae in Donaldson-Witten theory and integrable hierarchies_ ,  Adv. Theor. Math. Phys. 4, 503 (2000),  10.4310/ATMP.2000.v4.n3.a1,  hep-th/0006113.
  * [74] N. Nekrasov,  _Blowups in BPS/CFT correspondence, and Painlevé VI_ (2020),  2007.03646.
  * [75] S. Jeong and N. Nekrasov,  _Riemann-Hilbert correspondence and blown up surface defects_ ,  JHEP 12, 006 (2020),  10.1007/JHEP12(2020)006,  2007.03660.
  * [76] C. A. Keller and J. Song,  _Counting Exceptional Instantons_ ,  JHEP 07, 085 (2012),  10.1007/JHEP07(2012)085,  1205.4722.
  * [77] J. Gu, B. Haghighat, K. Sun and X. Wang,  _Blowup Equations for 6d SCFTs. I_ ,  JHEP 03, 002 (2019),  10.1007/JHEP03(2019)002,  1811.02577.
  * [78] J. Gu, A. Klemm, K. Sun and X. Wang,  _Elliptic blowup equations for 6d SCFTs. Part II. Exceptional cases_ ,  JHEP 12, 039 (2019),  10.1007/JHEP12(2019)039,  1905.00864.
  * [79] J. Gu, B. Haghighat, A. Klemm, K. Sun and X. Wang,  _Elliptic blowup equations for 6d SCFTs. Part III. E-strings, M-strings and chains_ ,  JHEP 07, 135 (2020),  10.1007/JHEP07(2020)135,  1911.11724.
  * [80] J. Kim, S.-S. Kim, K.-H. Lee, K. Lee and J. Song,  _Instantons from Blow-up_ ,  JHEP 11, 092 (2019),  10.1007/JHEP11(2019)092,  [Erratum: JHEP 06, 124 (2020)],  1908.11276.
  * [81] J. Gu, B. Haghighat, A. Klemm, K. Sun and X. Wang,  _Elliptic blowup equations for 6d SCFTs. Part IV. Matters_ ,  JHEP 11, 090 (2021),  10.1007/JHEP11(2021)090,  2006.03030.
  * [82] Z. Duan, K. Lee, J. Nahmgoong and X. Wang,  _Twisted 6d (2, 0) SCFTs on a circle_ ,  JHEP 07, 179 (2021),  10.1007/JHEP07(2021)179,  2103.06044.
  * [83] K. Lee, K. Sun and X. Wang,  _Twisted Elliptic Genera_ (2022),  2212.07341.
  * [84] H.-C. Kim, M. Kim and Y. Sugimoto,  _Blowup equations for little strings_ ,  JHEP 05, 029 (2023),  10.1007/JHEP05(2023)029,  2301.04151.
  * [85] K. Sun,  _Blowup Equations and Holomorphic Anomaly Equations_ (2021),  2112.14753.
  * [86] P. Jefferson, S. Katz, H.-C. Kim and C. Vafa,  _On Geometric Classification of 5d SCFTs_ ,  JHEP 04, 103 (2018),  10.1007/JHEP04(2018)103,  1801.04036.
  * [87] J. Chen, Y. Lü and X. Wang,  _D-type Minimal Conformal Matter: Quantum Curves, Elliptic Garnier Systems, and the 5d Descendants_ (2023),  2304.04383.
  * [88] X. Wang, G. Zhang and M.-x. Huang,  _New Exact Quantization Condition for Toric Calabi-Yau Geometries_ ,  Phys. Rev. Lett. 115, 121601 (2015),  10.1103/PhysRevLett.115.121601,  1505.05360.
  * [89] A. Grassi, Y. Hatsuda and M. Marino,  _Topological Strings from Quantum Mechanics_ ,  Annales Henri Poincare 17(11), 3177 (2016),  10.1007/s00023-016-0479-4,  1410.3382.