# End-to-end Optimization of Optical Communication Systems based on Directly
Modulated Lasers

Sergio Hernandez  Department of Electrical and Photonics Engineering,
Technical University of Denmark, 2800 Kongens Lyngby, Denmark<EMAIL_ADDRESS>Christophe Peucheret  Univ. Rennes, CNRS, FOTON - UMR6082, 22305 Lannion,
France  Francesco Da Ros  Department of Electrical and Photonics Engineering,
Technical University of Denmark, 2800 Kongens Lyngby, Denmark  Darko Zibar
Department of Electrical and Photonics Engineering, Technical University of
Denmark, 2800 Kongens Lyngby, Denmark

###### Abstract

The use of directly modulated lasers (DMLs) is attractive in low-power, cost-
constrained short-reach optical links. However, their limited modulation
bandwidth can induce waveform distortion, undermining their data throughput.
Traditional distortion mitigation techniques have relied mainly on the
separate training of transmitter-side pre-distortion and receiver-side
equalization. This approach overlooks the potential gains obtained by
simultaneous optimization of transmitter (constellation and pulse shaping) and
receiver (equalization and symbol demapping). Moreover, in the context of DML
operation, the choice of laser-driving configuration parameters such as the
bias current and peak-to-peak modulation current has a significant impact on
system performance. We propose a novel end-to-end optimization approach for
DML systems, incorporating the learning of bias and peak-to-peak modulation
current to the optimization of constellation points, pulse shaping and
equalization. The simulation of the DML dynamics is based on the use of the
laser rate equations at symbol rates between 15 and 25 Gbaud. The resulting
output sequences from the rate equations are used to build a differentiable
data-driven model, simplifying the calculation of gradients needed for end-to-
end optimization. The proposed end-to-end approach is compared to 3 additional
benchmark approaches: the uncompensated system without equalization, a
receiver-side finite impulse response equalization approach and an end-to-end
approach with learnable pulse shape and nonlinear Volterra equalization but
fixed bias and peak-to-peak modulation current. The numerical simulations on
the four approaches show that the joint optimization of bias, peak-to-peak
current, constellation points, pulse shaping and equalization outperforms all
other approaches throughout the tested symbol rates.

††journal: jocn

##  1 Introduction

Directly modulated lasers (DMLs) are a compelling option for short-reach
intensity modulation/direct detection (IM/DD) systems, thanks to their low
energy consumption, small form factor and reduced cost [1, 2]. The aim in such
systems is to maximize the symbol rate $R_{s}$ while maintaining sufficient
optical received power $P_{\mathrm{rec}}$ in order to increase net data
throughput. This requires the laser to operate in the large-signal regime,
where the high modulation index allows high extinction ratio and peak-to-peak
output power. However, the large-signal DML dynamics introduce significant
waveform distortions, potentially resulting in nonlinear intersymbol
interference as the symbol rate is increased [3, 4]. This is due to the
modulation-induced changes in carrier and photon concentration within the
laser active region, that cause nonlinear memory effects in the output optical
field. Consequently, the DML’s response time sets a limit on its modulation
bandwidth, restricting the laser’s throughput [5].

Figure 1: Proposed numerical simulation of a back-to-back DML system. The DSP
pipeline at TX (green) and RX (red) aims to encode and decode information
subject to the constraints imposed by the simulated analog domain (yellow).
The elements optimized in this work are underlined

Although increasing the bias current to the laser can enhance modulation
bandwidth, it comes at the expense of increased energy consumption and lower
extinction ratio, resulting in a receiver power penalty. Another alternative
to partially overcome bandwidth limitations is tuning transmitter (TX),
receiver (RX), digital signal processing (DSP) and the laser-driving
configurations (bias and peak-to-peak current to the laser) separately. Yet,
the large amount of parameters to be optimized within the DSP pipeline (pulse
shaping, pre-distortion, receiver-side equalization) may require a high amount
of simulation/system evaluations to converge towards optimal configurations,
making such approach intractable in cases where such evaluations are time
consuming.

Recent advances in DML development have focused on delivering >100 Gbps bit
rates with low energy consumption (<1 pJ/bit) through the use of coupled-
cavity laser structures [6, 7, 8]. The additional optical feedback enabled by
such structures enables higher modulation bandwidths through the use of
photon-photon resonance (PPR) and detuned loading. Impairment compensation in
DML-based systems [1, 9] has however relied mainly on separate transmitter-
side pre-distortion and receiver-side optimization, omitting the potential
gains of optimizing both TX and RX jointly. This arises from the fact that the
large-signal DML dynamics are governed by nonlinear differential equations,
hindering the calculation of analytical gradients and therefore the
simultaneous optimization of transmitter and receiver using standard gradient-
base optimization techiques. An accurate, differentiable model of the DML
dynamics can prove useful in this scenario, allowing the propagation of
gradients between TX and RX while accounting for the DML-induced signal
impairments. Data-driven modeling is a viable option in this context provided
that sufficient data is available, yet the choice of model structure (Volterra
filters, neural networks) can pose a challenge. In [10], we conduct a
comparison between model structures using the laser rate equations as source
of DML waveform data. The results show the potential of transformer-based
neural networks [11] in the prediction of the DML dynamics.

DSP algorithms have extensively been used to compensate for transmission
impairments within optical communication systems. End-to-end (E2E) learning
has attracted a special interest in this scope, as it enables the simultaneous
optimization of TX (constellation and pulse shaping) and RX, given that a
differentiable model of the system under test is available [12, 13]. Through
the use of alternative approaches like gradient-free optimization and
reinforcement learning, it is also possible to avoid modeling the
communication channel and use local gradient approximations [14, 15]. The
focus in E2E learning is to substitute one or several functions of the DSP
pipeline at both ends of the channel by an adaptive optimizable function [16].
Autoencoders (AEs) are especially popular in this context, as they allow the
compression and decompression of data [17] in a similar fashion to how
communication systems map symbols to optical waveforms and viceversa [16].
Approaches such as geometric constellation shaping (GCS), that have become
standard in state-of-the-art in long-haul coherent systems, can be implemented
based on this principle [18, 19, 20, 21, 15].

In this paper, we demonstrate E2E optimization of TX and RX DSP, together with
bias and peak-to-peak modulation current for impairment compensation of DML
systems. The proposed simulation approach is shown in Fig. 1. The modeling of
the DML dynamics is performed using a data-driven surrogate model, based on
the numerical solution of the rate equations as source of data. The waveform
data generation uses 4-level pulse amplitude modulation (4PAM) symbols in a
back-to-back (B2B) simulation. The proposed AE-based approach is used to
optimize DSP configurations at TX (GCS, pulse shaping) and RX (equalization,
symbol detection) simultaneously. In addition, the input current offset
$I_{\mathrm{bias}}$ and the peak-to-peak modulation current $I_{\mathrm{pp}}$
are used as learnable parameters. This provides insight on the optimal
compromise between extinction ratio and distortion of the optical waveform.
Additionally to the AE approach, three more approaches are included as a
benchmark: the uncompensated system without equalization, a receiver-side
linear feed-forward equalizer (FFE) and a second E2E approach with learnable
pulse shape (LPS) and RX-side Volterra nonlinear equalization (VNLE), but
excluding $I_{\mathrm{bias}}$ and $I_{\mathrm{pp}}$ from the optimization. The
results show the advantage of the joint optimization of bias, modulation
current and DSP for DML systems, with the AE yielding performance gains over
the RX-only equalization and a considerable advantage over the VNLE setup in
terms of symbol error rate (SER) and mutual information (MI).

This paper is structured as follows: Section 2 describes the dynamic behaviour
of DMLs and the state-of-the-art in optical communication laser development.
Section 3 analyzes the available algorithms for the modeling of DMLs in terms
of complexity, interpretability and compatibility with gradient-based
optimization. The fundamentals of E2E learning and its application to DMLs in
the literature is reviewed in Section 4. The motivation and structure of the
proposed DML model and the optimization approach around it can be found in
Section 5. The results of the E2E-optimized B2B DML system simulation are
described in Section 6. The conclusions are summarized in Section 7.

##  2 DML dynamics

The rate equations governing the photon density S(t), carrier density N(t) and
phase $\phi(t)$ in DMLs are given by [5]:

$\displaystyle\frac{dS(t)}{dt}=\Gamma g_{0}(N(t)-N_{0})\frac{1}{1+\epsilon
S(t)}S(t)-\frac{S(t)}{\tau_{p}}+\frac{\Gamma\beta N(t)}{\tau_{n}}\,,$ (1)
$\displaystyle\frac{dN(t)}{dt}=\frac{I(t)}{qV}-\frac{N(t)}{\tau_{n}}-g_{0}(N(t)-N_{0})\frac{1}{1+\epsilon
S(t)}S(t)\,,$ (2)
$\displaystyle\frac{d\phi(t)}{dt}=\frac{1}{2}\alpha\left[\Gamma
g_{0}(N(t)-N_{0})-\frac{1}{\tau_{p}}\right]\,,$ (3)

where $\Gamma$ is the mode confinement factor, $g_{0}$ is the gain slope
constant, $N_{0}$ is the carrier density at transparency, $\epsilon$ is the
gain compression factor, $\tau_{p}$ is the photon lifetime, $\beta$ is the
fraction of spontaneous emission coupled into the lasing mode, $\tau_{n}$ is
the electron lifetime, $I(t)$ is the injected current, $q$ is the electron
charge, $V$ is the active layer volume and $\alpha$ is the linewidth
enhancement factor. The output optical power is given by:

$P(t)=\frac{S(t)V\eta_{0}h\nu}{2\Gamma\tau_{p}}\,,$ (4)

where $\eta_{0}$ is the differential quantum efficiency, $h$ is Planck’s
constant, and $\nu$ is the unmodulated optical frequency. Thus, the modulated
optical power $P(t)$ is proportional to the photon density $S(t)$, while the
rate of change of the optical phase (instantaneous angular frequency) is
proportional to the carrier density over threshold.

###  2.1 Small-signal regime

The term $N(t)-N_{0}$ in Eq. 1 and 2 hints the reservoir-like behaviour of the
carrier density fluctuation in the laser active region. When $I$ reaches the
laser’s threshold current level $I_{\mathrm{th}}$ (optical gain overcomes
cavity losses) the carrier density reaches its maximum steady-state level
$N_{0}$ (the reservoir becomes full). The carrier excess generated by
$I-I_{th}$ "overflows" the reservoir, generating net stimulated emission of
photons. The higher $I-I_{th}$ is, the more photons are emitted, and through
Eq. 4 the larger $P(t)$ will be. This also entails that the stimulated
recombination of carriers increases, bringing $N(t)$ back to its threshold
level $N_{0}$ after some transient period. The interaction between carriers
and photons creates a damped oscillatory behaviour in $S(t)$ and $N(t)$
(carrier-photon resonance) when $I(t)$ is modulated, as the increase of one
drives the decrease of the other until steady state is reached. The optical
gain generated by the laser, dependent of $N(t)$ and $S(t)$, is given by Eq.
5:

$g(N,S)=\frac{g_{0}}{1+\epsilon
S(t)}\ln\left(\frac{N(t)+N_{s}}{N_{0}+N_{s}}\right)\,,$ (5)

where $N_{s}$ is a fitting parameter used to ensure the logarithm is finite
and defined for $N(t)>0$. Both $g_{0}$ and $N_{s}$ are usually fitted from the
measured response of the laser [22]. Given that the optical gain $g$ is
monotonically (although non-linearly) related to $N$ and $S$, it can be
assumed that for small variations of the carrier density $\Delta N$ the
optical gain variation is proportional to $\Delta N,\Delta S$. This leads to
$g=g_{th}+a\Delta N-a_{p}\Delta S$, using the local slopes $a=\partial
g/\partial N$ (differential gain) and $a_{p}=\partial g/\partial S$. The
approximation yields accurate results as long as
$I_{\mathrm{pp}}<<I_{\mathrm{bias}}$ for $I>I_{th}$, the small-signal regime
conditions. Although such conditions are usually impractical in a
communication setting, the small-signal analysis provides valuable insight
into the characteristics of a DML. It can be shown that, under the small
signal approximation and under sinusoidal excitation at angular frequency
$\omega$, the frequency response of the laser can be expressed according to:

$H(\omega)=\frac{\omega_{R}^{2}}{\omega_{R}^{2}-\omega^{2}+j\omega\gamma}\,,$
(6)

where $f_{R}$ is the natural resonant frequency of the system, $j$ is the
imaginary unit and $\gamma$ is the laser’s damping factor. Assuming operation
over threshold, $f_{R}$ can be approximated as [22]:

$f_{R}\approx\frac{1}{2\pi}\sqrt{\frac{v_{g}a\bar{S}}{\tau_{p}}}=\frac{1}{2\pi}\sqrt{\frac{\Gamma
v_{g}a}{qV}\eta_{i}(I_{\mathrm{bias}}-I_{th})}\,,$ (7)

where $\bar{S}$ denotes the average photon density. $f_{R}$ is also called
relaxation-oscillation or carrier-photon resonance frequency, as it determines
the frequency of the transient damped oscillations created in $S(t)$ and
$N(t)$ due to the aforementioned carrier-photon interaction under current
modulation. $\gamma$ can be expressed as:

$\gamma=Kf_{R}^{2}+\gamma_{0}\,,$ (8)

where $K$ is the "K-factor" governing the laser response at high modulation
frequencies and $\gamma_{0}$ is the damping offset. Both parameters are
usually obtained through fitting from the modulation response curve (Eq. 6).

Figure 2: Small-signal modulation response of a simulated DML. The dotted line
shows the 3 dB bandwidth level

The linear relation between $\omega_{R}^{2}$ and $(I_{\mathrm{bias}}-I_{th})$
describes the impact of $I_{\mathrm{bias}}$ in the modulation response of the
laser. Such relation is shown in Fig. 2, where various frequency responses
$|H(\omega)|$ are represented for different values of $I_{\mathrm{bias}}$. It
becomes apparent how the response shows a peak (that coincides approximately
with $f_{R}$) and a rapid damping for frequencies beyond it. A 3-dB bandwidth
$f_{3dB}$ can be defined as the frequency for which the magnitude of the
modulation response decreases to half of its DC value. For small $\gamma$, the
approximation $f_{3dB}\approx 1.55f_{R}$ is often used. The small-signal
intensity modulation (IM) of DMLs is therefore dictated by the
$I_{\mathrm{bias}}$ level, with higher levels associated with a higher 3dB-
bandwidth, although with diminishing returns due to the
$(I_{\mathrm{bias}}-I_{th})^{1/2}$ factor in Eq. 7.

###  2.2 Large-signal regime

The cost and power constraints in short-reach IM/DD systems favor amplifier-
free operation. The large-signal regime describes the DML laser behaviour in
most communication settings, where a high modulation index is desirable to
overcome receiver noise while avoiding the use of amplifiers. This is due to
the analytical solutions to the rate equations being unavailable in this
regime (where $I_{\mathrm{pp}}$ and $I_{\mathrm{bias}}$ are comparable in
magnitude). The most direct effect of a large $I_{\mathrm{pp}}$ is the
introduction of signal intensity distortion due to the effect of relaxation
oscillations. This effect is depicted in Fig. 3a, where overshoot and
undershoot can be observed right after changes of the modulating current
value. The larger instantaneous variation of $N,S$ enhances the amplitude of
the relaxation oscillations, introducing signal components uncorrelated to
$I$. The instantaneous emission frequency $\nu$ can also show significant
differences with respect to the small signal regime, as shown in Fig. 3b. The
figure shows the variation of $\nu$ resulting from a current modulation. This
is due to the signal chirp and relaxation oscillations resulting from the
higher modulation index. Frequency chirping can induce a critical limitation
of the transmission distance over a dispersive optical fiber. This is due to
the interaction between chirp and CD in optical fiber, that induces pulse
broadening on the envelope of the received optical signal. The effect of
chromatic dispersion on a signal with nonlinear chirp is depicted in Fig. 3c,
d, where the eye diagram of the simulated optical signal before and after
transmission over 2 km of standard single-mode fiber (SSMF) is shown. The
detrimental impact of pulse broadening is twofold: it decreases the peak
optical power of the pulses while introducing ISI between neighboring symbols.
It therefore entails substantial worsening of the effective SNR on top of the
linear fiber-induced attenuation. The instantaneous chirp expression [23, 24]
is obtained from the derivative of Eq. 3:

Figure 3: Waveform distortion introduced by DMLs: a) relaxation oscillations,
b) instantaneous chirp, c) simulated back-to-back eye diagram under 4PAM
modulation and d) simulated eye diagram after 2km of SSMF, $R_{s}=30$ Gbaud,
central wavelength $\lambda=1500$ nm. The power levels
$P_{0},P_{1},P_{2},P_{3}$ and instantaneous lasing frequency levels
$\nu_{0},\nu_{1},\nu_{2},\nu_{3}$ correspond to the 4 symbols in the 4PAM
intensity modulation of $I$.

$\Delta\nu=\frac{\alpha}{4\pi}\left[\frac{1}{P(t)}\frac{dP(t)}{dt}+\kappa
P(t)\,\right],$ (9)

and the factor $\kappa$ is defined by:

$\kappa=\frac{2\Gamma}{\eta h\nu_{0}V}\epsilon\,,$ (10)

where $h$ is Planck’s constant. The first factor in the sum of Eq. 9 is known
as transient chirp, while the second is called adiabatic chirp. The prevalence
of each chirp factor depends on the dynamics of $P$, and therefore on the
variation of the input current $I$. The derivative term of the transient chirp
makes it sensitive to fast variations of the optical power. This makes it
especially prevalent in high $R_{s}$ conditions, where the rapid variation of
power makes its derivative large. The use of sharp pulses, like square pulses,
induce high instantaneous derivative and high relaxation oscillation
amplitude, leading to additional transient chirp in the output waveform. The
adiabatic chirp becomes prevalent under high $I_{\mathrm{bias}}$ and optical
output powers, especially at lower $R_{s}$. Chirp in DMLs is therefore
unavoidable, but the careful selection of certain configuration parameters
($I_{\mathrm{pp}},I_{\mathrm{bias}}$, pulse shape) could help mitigate it.

The limited modulation bandwidth of DMLs is associated with several factors,
including the aforementioned carrier-photon interaction. Such factors
introduce linear and nonlinear memory effects, that become apparent as $R_{s}$
increases. Although Eq. 6 gives meaningful insight on the characteristics of
the DML, the DML response depends on the laser configuration, including
$I_{\mathrm{bias}},I_{\mathrm{pp}}$, and pulse shaping, and can only be
obtained through numerical simulations. Additionally, the susceptibility of
the system to timing and amplitude impairments varies depending on the channel
characteristics. Thus, the use of eye diagrams (Fig. 3c and d) is widely used
in the evaluation of experimental DML-based systems in the large signal
regime, as it allows to qualitatively assess the signal degradation at the
receiver.

Figure 4: Reflectivity vs detuning of a DBR. The modulation-induced chirp
drives the lasing wavelength towards the maximum reflectivity peak

###  2.3 Coupled-cavity effects

Through the use of additional passive and active laser sections in the laser
structure, it is possible to overcome the physical limitations that hinder the
DML modulation bandwidth beyond 40 GHz [25, 26, 27, 28, 22]. Among the various
techniques used, detuned loading and photon-photon resonance (PPR) have
delivered promising results in several implementations. Detuned loading is
based on the combined effect of distributed reflectors, mostly distributed
Bragg reflectors (DBRs) and modulation-induced chirping in order to enhance
the dynamic response of the DML [29, 30]. This is achieved through the design
of the gain and reflective sections, so the range of main lasing wavelengths
is detuned from the main lobe of the DBR reflectivity spectrum. Fig. 4 shows
the instantaneous lasing frequencies $\nu_{0},\nu_{1}$, associated with 2
different pump current levels $I_{0},I_{1}$. Both frequencies fall on
significantly different regions of the DBR reflectivity spectrum, leading to a
difference in mirror loss and the effective differential gain $a_{eff}$. The
approximation of $f_{R}$ in Eq. 7 must therefore be modified to account for
the detuning $\Delta\lambda$ from the DBR central wavelength:

$f_{R,DL}\approx\frac{1}{2\pi}\sqrt{\frac{v_{g}S}{\tau_{p}}\Re\left[\tilde{\Gamma}_{z}(\Delta\lambda)g_{0}(1+j\alpha)\right]}\,,$
(11)

where $\tilde{\Gamma}_{z}(\Delta\lambda)$ is the complex reactive confinement
factor in the active section [30]. The larger the imaginary part of this
factor becomes, the larger the impact of $\alpha$ (proportional to chirp)
becomes, enhancing $f_{R}$. When $\Delta\lambda=0$,
$\tilde{\Gamma}_{z}(\Delta\lambda)$ becomes real, and therefore no $f_{R}$
advantage is obtained from the passive section. Detuned loading is therefore
an effective method to exploit the FM derived from direct modulation by
enhancing the IM, resulting in higher modulation bandwidth. This is not only
desirable from the IM bandwidth perspective, but it can also mitigate chirp,
leading to reduced distortion after fiber transmission [31].

PPR bandwidth enhancement is based on the interaction between longitudinal
modes in the laser cavity. When such modes are close in wavelength, they can
still be coupled, [32, 33] leading to a potentially faster modulation response
compared to the single-mode configuration. This is usually done through
amplification of a secondary side mode using a coupled reflective cavity, as
in the case of detuned loading. Both approaches are mutually compatible,
resulting in further enhancements in modulation bandwidth. In the case of PPR,
a key design goal is to tune the grating phase to match the round-trip phase
of the main mode and the side mode, forcing constructive interference between
themselves. In this fashion, $\Gamma$ shows also a time dependency, due to the
time-varying interaction between the modes in the cavity. The combined gain of
the two modes makes possible to force a second resonance peak, of higher
frequency that the carrier-photon resonance, to appear in the modulation
response of the DML.

###  2.4 High-speed DMLs

The modulation performance of DMLs, based on both in-plane lasers [2, 26, 34,
35, 36, 37] and vertical cavity surface-emitting lasers (VCSELs) [38, 39, 40,
41, 42] structures, has been extensively developed through cavity design. It
must be noted that the commercial interest of DMLs often resides in their
energy efficiency and cost, and considerations like reliability, output power
and thermal performance play a role in their development. Their comparison
must therefore go beyond modulation bandwidth, and a more holistic assessment
must be made [2]. Within the cavity design choice, several structures and
substrates have been employed in the development of high-bandwidth ($f_{3dB}$)
DML lasers. >50 GHz bandwidth was achieved in [43] using push-pull modulation,
where the current is injected on two locations within the cavity. The driving
of the laser is designed so an increase in one of the currents leads to a
decrease in the other, and viceversa. Through the combined use of a MQW active
region and a single DBR section enhancing detuned loading and PPR, uncooled
>100 GHz bandwidth was obtained in [8], maintaining >70GHz at 85°C. The
combination of detuned loading and PPR has been used in several works,
obtaining $f_{3dB}$ between 40 and 65GHz [6, 7, 1, 44]. As a general overview,
several configurations have achieved bandwidths over 50GHz, but 100 GHz have
been reached using more complex and costly fabrication processes [2].

##  3 System Modeling

Modeling a dynamic system is a fundamental first step to its optimization. In
the case of DMLs, several alternatives are available depending on the
computational resources available, time constraints and desired accuracy and
exhaustiveness of the modeling. In this section we introduce some of the most
common approaches, sorting them from the more physics-intensive ones to the
purely data-driven.

Figure 5: P-I curve associated with the simulated DML. Depending on the
physical structure of the laser and thermal effects, the power response of the
laser might show nonlinearities

###  3.1 Parameter extraction and ODE solvers

Given that the previously introduced rate equations (Eq. 2, 1 and 3) yield
accurate prediction of the DML response in line with experimental verification
[45, 46, 47], they may seem as the most direct way of estimating the system
response. However, they entail 3 main challenges:

  * •

Most of their parameters are not easily measurable

  * •

They require the use of numerical ordinary differential equation (ODE) solvers
in the large signal regime

  * •

Analytical gradients are unavailable due to the computational structure of ODE
solvers

The problem of parameter extraction has been extensively treated in the
literature [45, 46, 47, 48]. Iterative techniques allow to obtain all the
relevant parameters based on experimentally obtained data, like the lasing
spectrum, the modulation response (Fig. 2), the static light-current (L-I)
curve (Fig. 5) or the spontaneous emission spectrum of the laser. Through the
use of machine learning, it is possible to automatize this process, simply by
providing the algorithm with the necessary figures of merit, avoiding manual
calculation of parameters [49]. Yet, quantities like $N_{0}$ are not directly
measurable, and must be estimated through related parameters, leading to
potential inaccuracy.

Even when all the necessary parameters are available, the use of ODE solvers
may prove impractical in some scenarios. The large values of $S$ and $N$ can
lead to instability in the numerical calculation, leading to divergence and
failure of the method. This makes convergence analysis a must have when using
ODE solvers, and re-initialization of the calculation may be needed [50]. The
sequential nature of the widely used Euler-based methods [51] makes
parallelization challenging, leading to computational bottlenecks depending on
the analyzed system. The largest advantage of such methods is their
configurable precision, but the computational toll associated with high
precision must also be considered. Lastly, ODE solvers make use of local
gradient approximations in order to solve differential equations. This
collides with the concept of automatic differentiation that many machine
learning and optimization frameworks use, where systems are usually built on
functions with known analytical derivatives [52, 53], making gradient
calculations faster. This drawback makes an automatically differentiable model
of DMLs desirable, as it may allow the simultaneous optimization of large
parameter spaces in a relatively reduced computation time. The next
subsections will describe some of the available alternatives for the modeling
of DMLs.

###  3.2 Circuit-equivalent models

Circuit-equivalent models [50, 54, 55, 56, 57, 58, 59] represent a less
complex alternative to rate equation solvers, while sharing the
interpretability of physics-based models. They make use of electrical
components (amplifiers, capacitors, resistors…) to model the dynamic behaviour
of the laser by relating their characteristics to the laser’s rate equation
parameters. Thus, they are able to reproduce the behaviour of a wide variety
of lasers, as long as the corresponding parameters can be measured. The use of
current and voltage signals as sole dynamic variables makes them an ideal
candidate for circuit simulators, thus facilitating the task of building a
numerical simulation. Some literature has included the Langevin phase and
intensity noise sources in the simulation, [60, 61] thus introducing
significant acceleration with respect to noisy rate equations, where
stochastic local gradients can worsen the performance of adaptive step
solvers.

On the other hand, circuit-based models share some of the pitfalls of ODE
solvers. Firstly, although circuit simulators are optimized for the task,
exhaustive models can be relatively complex, introducing significant
computational burden to the simulation [59]. This becomes especially true in
optimization scenarios, where circuit-specialized software cannot be utilized,
and the numerous differential relations between variables must be obtained
numerically, i.e. with ODE solvers. In conclusion, although they can be useful
in the design and simulation of DMLs, circuit-equivalent models are so far
impractical as part of a gradient-based optimization pipeline, and
automatically differentiable alternatives should be developed for this purpose
[12].

Figure 6: Relation between complexity and interpretability of the proposed
modeling techniques

###  3.3 Interpretable data-driven modeling

The mathematical modeling of dynamic systems (based on underlying physics,
data, or a mixture of both) has been studied extensively in the control theory
community, in what is usually called system identification [62, 63]. A
fundamental advantage of data-driven modeling is its compatibility with
numerical optimization techniques: they are mostly based on continuous, easily
differentiable functions that provide numerical stability in the search of
optimal configurations. Some of the most popular approaches for discrete
linear systems are variations of autoregressive such as integrated moving
average models (ARIMA) [64, 65, 66]. These models combine the use of 3
techniques: the evaluation of past outputs to predict future ones (AR), the
discrete differentiation of the time series in order to force its stationarity
(I) and a moving average (MA) of the past prediction errors. However, due to
the nonlinear nature of the DML behaviour, many classical system
identification models cannot be applied. Additionally, ARIMA models are based
on endogenous time series, i.e. the sequence they make predictions on is
determined solely by its past values. In order to model nonlinear dynamics,
higher-order temporal dependencies and/or nonlinear activations are often
used. Nonlinear autoregressive exogenous models (NARX) [67, 68] are a common
choice in this case, as they combine linear operations with a static
nonlinearity that allows to model a larger function space. Defining an input
$u(n)$, an estimated output $\hat{y}(n)$, a input feature vector $x(n)$, a
learnable parameter space $\theta$, a and a static nonlinear function $h$, a
generic NARX model expression can be obtained:

$x(n)=[\hat{y}(n-1),...\hat{y}(n-N_{y}),u(n),...,u(n-N_{u})]\,,\\\ $ (12)
$\hat{y}(n,\theta)=h(x(n),\theta)\,,$ (13)

where $N_{u},N_{y}$ are the input and output memory length, respectively. This
general definition includes a large variety of functions, with different
nonlinearities, recursivity schemes and complexity levels.

Volterra filters are often used in communication-oriented DSP to model and/or
compensate for the dynamics of nonlinear system [69, 70, 71, 72]. They are
among the simplest NARX models, as they do not employ recursion, and
nonlinearity is achieved through multiplication of past input samples with
each other. They are based on Volterra series [73], a variation of the Taylor
series where instead of evaluating the analyzed function around a single
point, it is evaluated over an infinite range of past samples. Given the
previous definitions of $\hat{y}(n),u(n)$ the discrete-time Volterra series
can be evaluated as:

$\displaystyle\hat{y}(n)=h_{0}+\sum_{k_{1}=0}^{\infty}h_{1}(k_{1})u(n-k_{1})+$
(14)
$\displaystyle\sum_{k_{1}=0}^{\infty}\sum_{k_{2}=0}^{\infty}h_{2}(k_{1},k_{2})u(n-k_{1})u(n-k_{2})+...$

where $h_{0},h_{1},...,h_{n}$ are the Volterra kernels and
$k_{1},k_{2},...,k_{n}$ are the delays associated with every kernel order
$1,2,...n$. As in the case of the Taylor series, increasing the order of the
series delivers potentially higher accuracy in the representation, at the cost
of increasing complexity [74, 75]. This forces real implementations to
truncate the Volterra series in both kernel order and delay, leading to
Volterra filters. It must be noted that the Volterra filter of order 1
corresponds to the linear convolution operator, and therefore the Volterra
kernels can be interpreted as a higher-order impulse response of the analyzed
system. Fig. 6 shows a comparison in terms of interpretability and potential
complexity between data-driven and physics-driven modeling methods. The main
challenge when using Volterra filters is the calculation of optimal kernels.
There are several techniques developed for this purpose [76, 77, 78], but
numerical multivariate optimization techniques can also be used at the expense
of higher computation time.

Hammerstein-Wiener (H-W) models [67] are an alternative to NARX-based system
identification, and provide design flexibility while including a larger number
of hyperparameters to be tuned. They are based on the combination of an input
static nonlinear function $h$, a linear infinite impulse response (IRR) filter
$z(n)$ defined by the real or complex coefficients $p_{i},q_{i}$ and an output
static nonlinear function $g$. Its mathematical expression is:

$w(n)=h\left[u(n)\right]\,,$ (15)
$z(n)=\sum_{i=1}^{N_{z}}q_{i}z(n-i)+\sum_{j=0}^{N_{u}}p_{i}w(m-j)\,,$ (16)

$\hat{y}(n)=g\left[z(w(n))\right]\,,$ (17)

where $N_{w},N_{z}$ are the number of input and output delays, respectively.
The choice of $h,g$ and the delay orders $N_{w},N_{z}$ must be tailored to the
specific system to be modeled. Although the nonlinearity choice is usually
hyperparameter-optimized among a pool of usual functions (sigmoidal, piece-
wise continuous, wavelet), the linear transfer functions can be estimated
through gray-box modeling (combination of a-priori physical knowledge with
data-driven approaches), numerical optimization or stochastic models like the
expectation-maximization (EM) algorithm can also be used [79]. One of the main
advantages of H-W models resides in their interpretability, given that the
linear transfer function governing them can be analyzed with traditional
spectral analysis as any other linear filter (Laplace/Z transforms, etc). When
the physics of the nonlinear system to be modeled are unknown or too complex
to be represented in such structure, purely data-based approaches, like neural
networks (NNs), are a popular approach [63].

###  3.4 Neural networks

NNs are mathematical structures where complex calculations would be performed
through the combination of large numbers of simpler operations [17]. Their
potential resides in the use of relatively simple linear functions followed by
a nonlinear activation. Although, as in the case of Volterra filters, they can
be considered a NARX model, their advantage resides in their structural
flexibility, that allows them to specialize depending on the task to be
performed. The concatenation of simple mathematical operations within NNs
leads to large, usually non-orthogonal parameter spaces. Such spaces must be
optimized to meet a certain objective, abstracted into a loss or cost
function. This is usually done through gradient-based numerical optimization
techniques, where instead of calculating the overall function gradient
analytically, it is approximated based on reduced subsets of the training data
(mini-batches). This procedure accelerates the overall gradient calculation,
leading to high performance models without explicit physical knowledge of the
system under investigation. The main inconvenient of this procedure is however
their lack of interpretability: the concatenation of nonlinearities within NNs
them makes it difficult to understand the interaction between the different
parameters in the network.

NNs have been widely used in the modeling of optical communication systems and
subsystems [80, 81, 82, 83]. Many of the NN architectures in the field have
tried to embed some temporal context into the network, as many of the systems
found in communications have memory elements. Time-delay neural networks
(TDNNs) [84, 85] are a special case of feedforward neural networks (FFNNs),
where no recursive element is added to network. Instead, the input of the
network is built as a sliding window of past samples, thus providing the
network with temporal context while maintaining relatively low complexity.
Single-dimensional convolutional neural networks (CNNs) [86] work on the same
principle as TDNNs, although they sometimes inherit pooling layers from its
higher-dimensional counterparts. Another possible approach to address this
problem is using recurrent neural networks (RNNs), where the network output
depends not only on the input features, but also on the state tensor in each
neuron. This tensor gives the network encoded temporal context on the past
inputs and outputs of the model, leading to potentially better prediction of
models with temporal dependencies. Gated recurrent units (GRUs) [87] and long-
short term memory (LSTM) [88] are two of the main architectures in this
paradigm. One major limitation of such networks is the modeling of long
temporal dependencies, due to the exploding and vanishing gradient problem
[89]. The transformer architecture has been extensively used due to its
success overcoming these problems in the natural language processing community
[90]. This has encouraged its use in the modeling of optic-fiber channels [91,
92]. Our work in [10] demonstrated the use of a variation of the transformer
architecture, the convolutional attention transformer (CAT) [11] in the
modeling of DMLs. Based on the structural similarities between neural ODEs and
the residual connections in transformers [93], we use the outputs of the DML
rate equations to obtain a differentiable DML model. The model compared
favourably to hyperparameter-optimized Volterra filters, TDNNs and LSTMs [10].

##  4 End-to-end Optimization

The throughput limitations of DMLs have also been studied from the DSP
perspective. In this scope, we can distinguish 2 different trends: one of them
aims to exploit the available bandwidth to increase spectral efficiency, while
the other aims to compensate the DML-induced waveform distortion. Discrete
multi-tone (DMT) [94], probabilistic constellation shaping (PCS) [95] and GCS
are the dominant technologies in the former category. DMT (the baseband,
guided-system equivalent of orthogonal frequency-division multiplexing, OFDM)
aims to exploit the available bandwidth through the use of multiple digital
narrow-band sub-carriers, instead of using a unique one spanning all the
available spectrum. In non-flat spectrum conditions, like the induced by PPR
in DMLs, the SNR varies along different regions of the spectrum. The use of
DMT allows higher granularity on the configuration for each sub-channel,
optimizing it to the local conditions of each narrow band. PCS aims to replace
uniform probability mass function (PMF) of the transmitted symbols by a
different PMF in order to optimize a certain metric (energy per bit, SNR). GCS
[20] follows similar approach, but it modifies the energy allocated to each
symbol instead of tweaking their probabilities. The combination of DMT and
constellation shaping is usually called entropy loading [96], and it allows to
adapt the symbol distribution to the spectral channel response, saving power
in the high-SNR narrow bands while maximizing it in the low-SNR bands. Many
DML implementations have included a combinations of these technologies to
increase the throughput of DMLs [97, 98, 28, 99]. Pulse shaping is another
source of throughput optimization. Faster-than-Nyquist (FTN) [100] signaling
aims to exploit this by introducing a controlled amount of ISI in each symbol,
correcting their introduced correlation on the RX side. This allows to
increase the effective throughput, as long as the receiver is able to
compensate both the TX- and the channel-induced ISI. [28] combines FTN with
entropy loading for increased spectral density of the DML bandwidth. Within
distortion compensation schemes, most approaches are focused on either
tailoring the peak-to-peak current to invert the nonlinear DML response or
using pre-distortion and equalization techniques. [101] provides a
linearization method to correct the DML nonlinearities, while [102] aims to
suppress relaxation-oscillations using ML approaches. Even optical filters can
be designed through semi-analytical approaches to find optimal cavity
configurations [4]. The development of DML-specific compensation approaches is
also extensive, using different equalizer structures like Volterra [103, 104],
TDNN [9], RNN [105] or deep belief networks [106].

Despite the remarkable performance leap obtained in some of the aforementioned
data-driven approaches, their optimization is based on single-sided or
sequential learning approaches: they do not optimize the TX and RX sides
jointly. In sequential learning, each side of the communication link (TX and
RX) is optimized keeping the configuration of the other one fixed, thus
avoiding the need for gradient propagation across the transmission channel.
E2E learning [16] aims to condense several functions within TX (constellation
shaping, pulse shaping, bit labelling) and RX (EQ, symbol decoding) into a
single neural network, where the link configuration is optimized
simultaneously at both ends based on a certain metric (loss function). Thus,
the performance of the configuration obtained through E2E learning is limited
mainly by the amount of DSP functions optimized and the accuracy of the
modeling of the elements between TX and RX. This can pose a considerable
challenge in the case of DML-based systems, where the building a
differentiable model of the DML is not straightforward, as discussed in the
previous section. Gradient-free approaches, based on derivative-free
optimizers or reinforcement learning, have been proposed in the training of
E2E [14, 107, 21]. Nonetheless, they introduce severe computational overhead
and their gradient approximations can lead to numerical optimization issues
[15]. Several gradient-based have been proposed for DML (VCSEL) systems,
implementing GCS [108] and pre-distortion + EQ [109] based on data-driven DML
models. The present work aims to jointly optimize transmitter GCS and LPS and
receiver EQ with the driving configuration of the DML
($I_{\mathrm{bias}},I_{\mathrm{pp}}$), thus tailoring E2E learning to the
specific characteristics of DML-based systems. This entails extending our DML
modeling work in [10], obtaining a model that is able to predict the dynamics
of the laser regardless of its biasing.

Figure 7: Architecture of the proposed AE

The usual approach to E2E learning in communication systems is using
autoencoders (AEs) [17] as substitutes of the DSP pipeline [16, 110, 111].
This is due to their functional resemblance to a communication link: they
compress information subject to physical constraints (encoding) to then
retrieve it with the minimum loss of information possible (decoding). The AE
was created to compress and decompress information with minimal information
loss, making it an ideal candidate for this task. It consists on an encoder,
that performs the dimensionality reduction, and a decoder, that aims to
retrieve the information compressed by the encoder (as shown in Fig. 7). AE-
based optimization of communication systems [110, 21] usually relies on an
encoder to map an alphabet of one-hot encoded symbols $i_{n}$ to a vector of
samples $h_{n}$ representing a symbol (assuming 2 sps in Fig. 7) to be
transmitted. The decoder then maps the received samples $g_{n}$ back to
scalars $o_{n}$ representing output probabilities, usually after a series of
hidden layers, represented by $l_{n}$ in Fig. 7. It must be noted that the
matrix multiplication in AEs introduces a significant complexity overhead to
the TX and RX operation. Therefore, AEs are often used only to find optimal
link configurations, that are then implemented in the form of less complex
look-up tables and decision circuits. This allows to take advantage of the
performance gains yielded by autoencoders while adjusting complexity to the
requirements of optical transceivers [112]. The standard AE [17] is an FFNN,
and therefore it lacks any memory mechanism in its structure. This explains
the use of a finite impulse response (FIR) filter on the decoder side, that
allows the AE to capture and compensate for ISI and undesired memory effects.
The $o_{n}$ are then masked with a normalized exponential function (softmax),
depicted in Eq. 18.

$\sigma(o_{n})=\frac{e^{o_{n}}}{\sum_{j=1}^{N}e^{o_{j}}}$ (18)

where $N$ is the number of neurons in the decoder output layer. The softmax
function converts the single-dimensional output tensor in an array of
normalized probabilities, thus giving insight on the certainty of the symbol
prediction. As in any supervised learning scheme, the AE needs to be trained
based on a certain loss function. Given that the goal of the link optimization
is minimizing the loss of information over the transmission channel, a metric
that captures this information loss is desirable. This leads to the use of
categorical cross-entropy (CE) as an indirect estimator of mutual information
between transmitted and received symbols [21]. The expression of CE is shown
in Eq. 19:

$J_{CE}(\theta)=\frac{1}{N}\sum_{n=1}^{N}\left[-\sum_{m=1}^{M}i_{m}^{(n)}\log
o_{m}^{(n)}(\theta)\right]\,,$ (19)

where $N$ represents the symbol batch size and $M$ represents the modulation
order utilized. In the backpropagation stage of the training process, the
gradient over trainable parameters $\theta$ is calculated to minimize CE. This
leads to higher mutual information between both ends of the link, due to the
approximation of its lower bound in Eq. 20:

$I(X;Y)\geq H(X)-\hat{H}(X|Y)\,,$ (20)

where $H(X)$ is the entropy of the probability distribution of the transmitted
symbols $X$ and $\hat{H}(X|Y)$ is the upper bound on the conditional entropy
of $X$ given the probability distribution of the received symbols $Y$.
$\hat{H}(X|Y)$ can be approximated by the CE between transmitted and received
symbols, given that the expression of the true channel transition is unknown
[21].

Parameter | Symbol | Value | Units
---|---|---|---
Confinement factor | $\Gamma$ | 0.24 | -
Photon lifetime | $\tau_{p}$ | 2.60 | ps
Carrier lifetime | $\tau_{c}$ | 3.17 | ns
Transparent carr. dens. | $N_{0}$ | $2.00\cdot 10^{24}$ | $m^{-3}$
Active cross section | $\sigma_{g}$ | $3.34\cdot 10^{-20}$ | $m^{2}$
Eff. refractive index | $n_{g}$ | 4 | -
Active volume | V | $3.60\cdot 10^{-17}$ | $m^{3}$
Gain compression | $\epsilon$ | $2.00\cdot 10^{23}$ | $m^{3}$
Spont. emiss. factor | $\beta$ | $1.00\cdot 10^{-3}$ | -
Diff. quant. eff. | $\eta_{0}$ | 0.20 | -
Table 1: Utilized parameters for the rate equation laser simulation Figure 8:
Block diagram of the proposed DML optimization approach. The trainable blocks
are highlighted in red

##  5 Simulation setup

The proposed DML system optimization is based on two main components: the
laser surrogate model and the E2E optimization approach. The surrogate model
aims solely to reproduce the laser behaviour as accurately as possible. The
E2E approach finds the optimal system configurations subject to the
constraints imposed by the surrogate model dynamics and its loss function,
once the surrogate model has been trained.

###  5.1 Surrogate model

The structure of the surrogate CAT model, depicted in Fig. 9, is based on 3
building blocks: learned positional embeddings, convolutional attention and
dense layers with ReLU activation. The positional embeddings aim to capture
the position of each sample within the sequence, giving temporal context to
the network. The convolutional attention block aims to discern the most
relevant input samples for the calculation of each output sample through the
use of matrix multiplication and learned convolutional filters. Lastly, the
dense layers provide depth to the network, as they contain most of the
trainable parameters within it. Although the model architecture, shown in Fig.
9, is identical to [10], in this work we introduce a new data generation
process. This allows to maintain model complexity with respect to our previous
work while achieving a DML model agnostic to the waveform shape, the bias, and
the peak-to-peak amplitude of the input sequence.

The surrogate modeling is performed in a data-driven fashion, using the laser
rate equations (Eq. 1, 2 and 3) as source of data. The dataset of input
waveforms to the rate equations is obtained by generating 4PAM-modulated
symbols, pulse shaped at 2 samples per symbol (sps). The rate equation
parameters utilized are specified in Table 1. The utilized parameters
correspond to a generic DFB laser with $f_{R}=10.6$ GHz and a
$f_{3\mathrm{dB}}$ of 25.5 GHz at $I_{\mathrm{bias}}=75$ mA. Although the
parameters are not chosen to match the specifications of a real device, they
are of the same order of magnitude as those reported in the literature for
single-mode quantum-well DFB lasers [113, 114]. It must be noted that
$\phi(t)$ is not considered, given that the simulated setup is back-to-back.
Therefore the values of $\alpha$ and $\kappa$ do not have an impact in the
behaviour of the surrogate model. The region of the modulation response beyond
$f_{R}$ is of special interest due to the increased data throughput and
significant waveform distortion introduced. The CAT model was therefore
trained separately on 3 relatively high $R_{s}$, namely $\\{15,20,25\\}$
Gbaud.

The generalizability in terms of waveform shaping is addressed by using 2
different types of pulse: square and stochastic pulses. The steep slope of the
square pulse-shaped symbols allows to capture the transient dynamics of the
DML response, while the stochastic pulses are used to prevent overfitting. The
stochastic pulses are implemented through a 2-tap finite impulse response
(FIR) filter, where the filter coefficients are drawn from a uniform
distribution in the interval $[-0.5,0.5]$. The pulses are normalized in
amplitude to avoid distorting the 4PAM symbols. In order to maintain a high
modelling accuracy regardless of the combination of $I_{\mathrm{bias}}$ and
$I_{\mathrm{pp}}$ utilized, both quantities are randomized for the generation
of the input waveform data set. This is achieved by randomizing
$I_{\mathrm{bias}}$, drawing it from a uniform distribution in the interval
$I_{\mathrm{bias}}=[50,100]$ mA. $I_{\mathrm{pp}}$ is randomized in a similar
fashion but constrained to the range $I_{\mathrm{pp}}=[0,80]$ mA, in order to
maintain a driving current significantly higher than the threshold current
$I_{th}=4.06$ mA for all $I_{\mathrm{bias}}$ and $I_{\mathrm{pp}}$
combinations. This will allow the optimization of the E2E learning approaches
in a wide variety of different input currents, allowing the use of
$I_{\mathrm{bias}}$ and $I_{\mathrm{pp}}$ as trainable parameters.

Figure 9: Architecture of the CAT. The network sub-blocks are described in
Section 5

Once the driving waveforms are generated, the sequences are then oversampled
to 32 sps and low-pass filtered (LPF) to constrain their bandwidth to
$0.9R_{s}$. This is done to ensure the accuracy and convergence of the
numerical solution to the rate equations, even if the surrogate model will
operate on a 2-sps basis. The target (ground truth) output sequences for the
CAT surrogate models are obtained by inputting the previously described
driving waveforms to a numerical solver of the laser rate equations at a
certain (randomized) $I_{\mathrm{bias}}$. The output photon density sequences
from the rate equations are then downsampled and converted to optical power
waveforms through Eq. 4. The loss function of the model is defined comparing
the sequences from the rate equation solver to the predictions of the
surrogate model using normalized root mean squared error (NRMSE). The use of
NRMSE provides better loss interpretability, as the normalization applied
makes the calculation of the loss relative to the amplitude of each generated
sequence. Once the surrogate model is trained, its weights are fixed and it is
used as part of an E2E optimization approach.

###  5.2 E2E approach

| BL/FFE | VNLE | AE
---|---|---|---
$I_{\mathrm{bias}}$ (mA) | 75 | 75 | Learnable
$I_{\mathrm{pp}}$ (mA) | $[8,80]$ | $[8,80]$ | Learnable
Pulse shap. | Square | Learnable | Learnable
GCS | Equispaced | Equispaced | Learnable
Filt. taps | 0/21 | 272 | 11
DML model | Rate eqs. | CAT | CAT
Loss | MSE | MSE | CE
Table 2: Parameters of the compensation approaches

The E2E optimization approach aims to find the optimal combination of TX, RX
and laser-driving parameters in order to minimize the probability of decision
errors after detection. The approach should therefore reach optimal
performance regardless of the targeted laser cavity structure, as long as the
surrogate DML model is able to reproduce the laser behavior accurately. Thus,
our approach paves the way for the E2E optimization of optical communication
system based on external modulators, like electro-absorption modulators
(EAMs), as long as data availability allows for the training of a sufficiently
accurate surrogate model. The investigated back-to-back IM/DD system is
represented in Fig. 8. The modulation is based on equiprobable 4PAM symbols,
upsampled to 2 sps. The TX establishes the GCS constellation intensity levels
and the 2-tap pulse shaping in order to generate the modulation current
sequences. The hardware limitations of the digital-to-analog and analog-to-
digital converters (DAC/ADC) are modeled as FIR LPFs, with a bandwidth
$B_{LPF}=0.9\cdot R_{s}$. The impact of impairments associated with DAC and
ADC is not considered in this paper. However, the AE is trained to take into
account the limited time resolution and BW of DAC and ADC. After filtering,
the sequences are constrained to the range $[-0.5,0.5]$ and then amplified to
$[-40,40]$ mA for a maximum peak-to-peak $I_{\mathrm{pp}}$ of 80 mA, in order
to match the surrogate data generation. $I_{\mathrm{bias}}$ is again
constrained to the range $[50,100]$ mA. During training, the emulation of the
DML response is carried out by the surrogate model, as its architecture allows
for automatic differentiation. The testing is however conducted on the laser
rate equations, using the same parameters as for data generation. This allows
to obtained more reliable performance metrics, as the surrogate may distort
them due to modeling inaccuracies.

Figure 10: Training and testing NRMSE scores of the CAT surrogate, using the
laser rate equations as ground truth Figure 11:  Symbol error rate (SER, top)
and mutual information (MI, bottom) results at symbol rates
$R_{s}=\\{15,20,25\\}$ Gbaud (left, middle, right) for the uncompensated
baseline (BL), finite impulse response equalizer (FIR EQ), Volterra nonlinear
equalizer (VNLE) and autoencoder (AE) setups

After generating the output sequences from the DML models, additive white
Gaussian noise (AWGN) samples are added to the signal after square-law
detection. The noise variance is fixed to yield 22 dB electrical signal-to-
noise ratio (SNR) at the highest average $P_{\mathrm{rec}}$ (corresponding to
$I_{\mathrm{pp}}=80$ mA). The value of $I_{\mathrm{bias}}$ does not impact
$P_{\mathrm{rec}}$, as only the peak-to-peak power is considered in it.
Lastly, the adaptive RX DSP performs 2-sps-wise EQ, downsampling and symbol
decision. The SER and MI performance metrics are obtained by comparing the
estimated symbol probabilities at RX with the originally sent sequence at TX,
using maximum likelihood detection (MLD) or softmax activation depending on
the approach. The compensation of the system is based on four different
approaches: the baseline uncompensated system (BL), an RX-side FIR FFE, a
second-order VNLE with a transmitter LPS and an AE. Table 2 lists the specific
parameters for each approach. The FFE is meant to emulate the conventional
approach to DML-based system optimization, while the VNLE and AE give insight
on the performance advantage of E2E models. It must be noted that only the AE
handles $I_{\mathrm{pp}}$ and $I_{\mathrm{bias}}$ as trainable parameters,
while all the other models are kept at fixed bias and swept through various
$I_{\mathrm{pp}}$ (and therefore $P_{\mathrm{rec}}$) values. Fig. 7 shows the
encoder and decoder AE architecture from input symbol $i_{n}$ to output symbol
likelihoods $o_{n}$. The encoder, or TX side, is based on a single linear
layer to map the input one-hot encoded vectors into 2-sample pulses. All the
obtained pulses are then serialized into 1024-sample time sequences to be used
as driving current to the DML. The $I_{\mathrm{bias}}$ level is added to the
signal after filtering and amplification, as explained previously. On the
decoder (RX) side, the received sequences are first FIR filtered in order to
provide a memory mechanism to the AE. After deserialization, 3 leaky-ReLU-
activated feedforward layers with softmax activation at the output convert the
filtered samples into symbol probabilities for decision. The AE loss is then
calculated using cross-entropy (CE) between the originally transmitted symbols
and their assigned symbol probability at the receiver.

##  6 Results

###  6.1 Surrogate model

The surrogate model training compromises $2^{23}$ samples in sequences of 1024
samples, while the testing data set includes $2^{17}$ samples. The
randomization of pulse shape and laser driving configurations is performed on
a per-sequence basis. The training and testing NRMSE scores obtained by the
CAT model are depicted in Fig. 10. Throughout the three analyzed $R_{s}$, the
testing NRMSE loss is slightly lower than its training counterpart, giving no
sign of overfitting. Another interesting trend is the slight increase in NRMSE
as $R_{s}$ increases. This trend might be related to the predominance of the
nonlinear effects at higher $R_{s}$, leading to potentially more complex DML
dynamics. In any case, the three models perform well below the $1\%$ mark.

###  6.2 E2E approach

The symbol dataset of the E2E approaches compromises $2^{20}$ symbols, with an
80/20 partition between training and validation, respectively. The dataset is
split in mini-batches of 512 symbols for exploiting parallelization to reduce
the training time. The $I_{\mathrm{pp}}$ to all the models except the AE is
swept in the range $[8,80]$ with a step of 8 mA between consecutive levels.

Even though all the approaches are subject to the same $I_{\mathrm{pp}}$
constraints, the E2E approaches are able to exploit the DML transient response
through LPS (and GCS in the case of the AE), yielding higher
$P_{\mathrm{rec}}$. The AE is able to optimize its $I_{\mathrm{pp}}$
dynamically within the constrained range, therefore only one AE power level
was analyzed. The training is iterated over three different symbol rates: 15,
20 and 25 Gbaud, matching the used in the surrogate training. The MI and SER
results tested on the laser rate equations are shown as a function of
$P_{\mathrm{rec}}$ in Fig. 11. Based on the figure, the AE delivers the best
SER and MI performance overall, hinting that the optimization of
$I_{\mathrm{bias}}$ level could have a higher impact on performance than the
equalization in certain cases. This is more accentuated on the lower $R_{s}$
than the higher ones, where the waveform distortion worsens the SNR, giving
equalization a higher relative impact. Another interesting metric is the
optimal $I_{\mathrm{bias}}$ obtained by the AE, resulting in
$\\{62.31,69.31,70.30\\}$ mA at $\\{15,20,25\\}$ Gbaud, respectively. This
trend provides interesting insight on the relation between the $R_{s}$ and the
optimal $I_{\mathrm{bias}}$ conditions. For the $R_{s}$ analyzed, the E2E
approaches show a clear performance advantage over the BL and RX-only
optimization. This is expected, given that the lack of nonlinearity in the
latter makes the complete compensation of the DML-induced distortion
infeasible at high symbol rates. The better performance of the E2E approaches
serves as a further validation of the surrogate, as a poor model of the DML
could lead to low performance when tested on the rate equations.

##  7 Conclusion

We propose a novel end-to-end, directly modulated laser optimization approach
using a differentiable data-driven model as laser surrogate, allowing the
propagation of gradients between transmitter and receiver. The surrogate is
built based on the laser rate equations using various bias and peak-to-peak
current values, in order to make it robust to such values. We compare
different system architectures with conventional receiver-side optimization,
varying the received optical power and the symbol rate of the simulation. The
proposed autoencoder approach including bias and peak-to-peak current
optimization shows significant performance gain compared to its receiver-side
counterpart, showcasing the potential of end-to-end approaches in the
optimization of directly modulated laser systems.

## Funding

The Villum Fonden (VI-POPCOM VIL54486) and Villum YIP OPTIC-AI (no. VIL29334)
projects are acknowledged.

## References

  * [1] N.-P. Diamantopoulos, T. Fujii, S. Yamaoka, H. Nishi, K. Takeda, T. Tsuchizawa, T. Segawa, T. Kakitsuka, and S. Matsuo, “60 GHz Bandwidth Directly Modulated Membrane III-V Lasers on SiO2/Si,” Journal of Lightwave Technology 40, 3299–3306 (2022).
  * [2] J. Huang, C. Li, R. Lu, L. Li, and Z. Cao, “Beyond the 100 Gbaud Directly Modulated Laser for Short Reach Applications,” Journal of Semiconductors 42, 041306 (2021).
  * [3] D. Che and X. Chen, “Modulation Format and Digital Signal Processing for IM-DD Optics at Post-200G Era,” Journal of Lightwave Technology 42, 588–605 (2024).
  * [4] M. S. Alam, R. Maram, K. A. Shahriar, P. Ricciardi, and D. V. Plant, “Chirped Managed Laser for Multilevel Modulation Formats: A Semi-analytical Approach for Efficient Filter Design,” Journal of Lightwave Technology pp. 1–11 (2023).
  * [5] L. A. Coldren, S. W. Corzine, and M. L. Mashanovitch, _Diode Lasers and Photonic Integrated Circuits_ (John Wiley & Sons, 2012).
  * [6] X. Pang, T. Salgals, H. Louchet, D. Che, M. Gruen, Y. Matsui, T. Dippon, R. Schatz, M. Joharifar, B. Krüger, F. Pittala, Y. Fan, A. Udalcovs, L. Zhang, X. Yu, S. Spolitis, V. Bobrovs, S. Popov, and O. Ozolins, “200 Gb/s Optical-Amplifier-Free IM/DD Transmissions Using a Directly Modulated O-Band DFB+R Laser Targeting LR Applications,” Journal of Lightwave Technology 41, 3635–3641 (2023).
  * [7] G. V. Rajeswari, M. Moehrle, F. Ehrensack, U. Troppenz, A. Sigmund, and M. Schell, “Novel >57 GHz bandwidth O-band InGaAlAs MQW RW DFB,” in _Novel In-Plane Semiconductor Lasers XXII,_ vol. 12440 A. A. Belyanin and P. M. Smowton, eds., International Society for Optics and Photonics (SPIE, 2023), p. 1244007.
  * [8] S. Yamaoka, N.-P. P. Diamantopoulos, H. Nishi, T. Fujii, K. Takeda, T. Hiraki, S. Kanazawa, T. Kakitsuka, and S. Matsuo, “Uncooled 100-GBaud Directly Modulated Membrane Lasers on SiC Substrate,” Journal of Lightwave Technology 41, 3389–3396 (2023).
  * [9] A. G. Reza, M. T. Costas, C. Browning, F. D. Otero, and L. Barry, “67.5 Gbit/s PAM-8 Signal Transmissions Over 25-km SMF With a 1550-nm 10G-Class DML Using Machine Learning,” in _Conference on Lasers and Electro-Optics,_ (Optica Publishing Group, 2022), p. SM3J.4.
  * [10] S. Hernandez, O. Jovanovic, C. Peucheret, F. D. Ros, and D. Zibar, “Differentiable Machine Learning-Based Modeling for Directly-Modulated Lasers,” IEEE Photonics Technology Letters 36, 266–269 (2024).
  * [11] S. Li, X. Jin, Y. Xuan, X. Zhou, W. Chen, Y.-X. Wang, and X. Yan, “Enhancing the Locality and Breaking the Memory Bottleneck of Transformer on Time Series Forecasting,” in _Proceedings of the 33rd International Conference on Neural Information Processing Systems,_ (2019).
  * [12] M. Srinivasan, J. Song, A. Grabowski, K. Szczerba, H. K. Iversen, M. N. Schmidt, D. Zibar, J. Schröder, A. Larsson, C. Häger, and H. Wymeersch, “End-to-End Learning for VCSEL-Based Optical Interconnects: State-of-the-Art, Challenges, and Opportunities,” Journal of Lightwave Technology 41, 3261–3277 (2023).
  * [13] L. Minelli, F. Forghieri, A. Shahpari, T. Shao, and R. Gaudino, “TDECQ optimization of VCSEL-MMF nonlinear Digital Pre-Distorters using End-to-end Learning,” in _European Conference on Optical Communication (ECOC),_ (2023), p. Th.B.5.6.
  * [14] F. A. Aoudia and J. Hoydis, “Model-Free Training of End-to-End Communication Systems,” IEEE Journal on Selected Areas in Communications 37, 2503–2516 (2019).
  * [15] M. P. Yankov, O. Jovanovic, D. Zibar, and F. D. Ros, “Recent Advances in Constellation Optimization for Fiber-optic Channels,” in _European Conference on Optical Communication (ECOC) 2022,_ (Optica Publishing Group, 2022), p. Mo3D.4.
  * [16] T. O’Shea and J. Hoydis, “An Introduction to Deep Learning for the Physical Layer,” IEEE Transactions on Cognitive Communications and Networking 3, 563–575 (2017).
  * [17] I. Goodfellow, Y. Bengio, and A. Courville, _Deep Learning_ (MIT Press, 2016).
  * [18] V. Oliari, B. Karanov, S. Goossens, G. Liga, O. Vassilieva, I. Kim, P. Palacharla, C. Okonkwo, and A. Alvarado, “Hybrid Geometric and Probabilistic Shaping; Is It Really Necessary?” in _Optica Advanced Photonics Congress 2022,_ (2022), p. SpTu1J.4.
  * [19] A. Soleimanzade, M. Ardakani, and H. Ebrahimzad, “Optimization of the Non-Linearity Tolerant 4D Geometric Shaped Constellations for Optical Fiber Communication Systems Using Neural Networks,” Journal of Lightwave Technology 42, 670–680 (2024).
  * [20] E. Sillekens, G. Liga, D. Lavery, P. Bayvel, and R. I. Killey, “High-Cardinality Geometrical Constellation Shaping for the Nonlinear Fibre Channel,” Journal of Lightwave Technology 40, 6374–6387 (2022).
  * [21] O. Jovanovic, F. Da Ros, D. Zibar, and M. P. Yankov, “Geometric Constellation Shaping for Fiber-Optic Channels via End-to-End Learning,” Journal of Lightwave Technology 41, 3726–3736 (2023).
  * [22] S. Yamaoka, N.-P. Diamantopoulos, H. Nishi, R. Nakao, T. Fujii, K. Takeda, T. Hiraki, T. Tsurugaya, S. Kanazawa, H. Tanobe, T. Kakitsuka, T. Tsuchizawa, F. Koyama, and S. Matsuo, “Directly modulated membrane lasers with 108 GHz bandwidth on a high-thermal-conductivity silicon carbide substrate,” Nature Photonics 15, 28–35 (2021).
  * [23] T. Koch and J. Bowers, “Nature of Wavelength Chirping in Directly Modulated Semiconductor Lasers,” Electronics Letters 20, 1038–1040 (1984).
  * [24] A. Villafranca, J. Lasobras, and I. Garces, “Precise Characterization of the Frequency Chirp in Directly Modulated DFB Lasers,” in _2007 Spanish Conference on Electron Devices,_ (2007), pp. 173–176.
  * [25] Y. Mao, Z. Ren, L. Guo, H. Wang, R. Zhang, Y. Huang, D. Lu, Q. Kan, C. Ji, and W. Wang, “Modulation Bandwidth Enhancement in Distributed Reflector Laser Based on Identical Active Layer Approach,” IEEE Photonics Journal 10 (2018).
  * [26] S. Sulikhah, S. L. Lee, and H. W. Tsao, “Improvement on Direct Modulation Responses and Stability by Partially Corrugated Gratings Based DFB Lasers With Passive Feedback,” IEEE Photonics Journal 13, 1–14 (2021).
  * [27] Y. Matsui, R. Schatz, T. Pham, W. A. Ling, G. Carey, H. M. Daghighian, D. Adams, T. Sudo, and C. Roxlo, “55 GHz Bandwidth Distributed Reflector Laser,” Journal of Lightwave Technology 35, 397–403 (2017).
  * [28] D. Che, Y. Matsui, R. Schatz, R. Rodes, F. Khan, M. Kwakernaak, and T. Sudo, “200-Gb/s Direct Modulation of a 50-GHz Class Laser With Advanced Digital Modulations,” Journal of Lightwave Technology 39, 845–852 (2021).
  * [29] U. Feiste, “Optimization of modulation bandwidth in DBR lasers with detuned Bragg reflectors,” IEEE Journal of Quantum Electronics 34, 2371–2379 (1998).
  * [30] M. Chaciński and R. Schatz, “Impact of Losses in the Bragg Section on the Dynamics of Detuned Loaded DBR Lasers,” IEEE Journal of Quantum Electronics 46, 1360–1367 (2010).
  * [31] Y. Matsui, R. Schatz, D. Che, F. Khan, M. Kwakernaak, and T. Sudo, “Low-chirp Isolator-free 65-GHz-bandwidth Directly Modulated Lasers,” Nature Photonics 15, 59–63 (2021).
  * [32] P. Bardella and I. Montrosset, “A New Design Procedure for DBR Lasers Exploiting the Photon–Photon Resonance to Achieve Extended Modulation Bandwidth,” IEEE Journal of Selected Topics in Quantum Electronics 19, 1502408–1502408 (2013).
  * [33] G. Morthier, R. Schatz, and O. Kjebon, “Extended Modulation Bandwidth of DBR and External Cavity Lasers by Utilizing a Cavity Resonance for Equalization,” IEEE Journal of Quantum Electronics 36, 1468–1475 (2000).
  * [34] S. Guan, Y. Zhang, J. Zheng, J. Su, Z. Sun, L. Lu, T. Fang, L. Li, R. Xiao, Y. Shi, and X. Chen, “Modulation Bandwidth Enhancement and Frequency Chirp Suppression in Two-Section DFB Laser,” Journal of Lightwave Technology 40, 7383–7389 (2022).
  * [35] K. Shinohara, R. Miyagoshi, Y. Suzuki, R. Suzuki, G. Sakaino, M. Shimada, and K. Matsumoto, “106-Gbps PAM4 Operation at an Extinction Ratio above 3.5 dB using a Conventional Buried-Heterostructure Directly Modulated Laser,” in _Optical Fiber Communication Conference (OFC) 2023,_ (2023), p. M2D.6.
  * [36] S. Ohno, M. Onga, T. Nakajima, A. Nakanishi, N. Sasada, S. Tanaka, R. Nakajima, and K. Naoe, “10-km Transmission of 106-Gb/s PAM4 with Directly Modulated DFB Lasers in the CWDM Range,” in _Optical Fiber Communication Conference (OFC) 2023,_ (2023), p. M2D.7.
  * [37] X. Zhu, J. Guo, H. Li, Z. Li, D. Zhou, L. Zhao, W. Wang, and S. Liang, “High Speed Directly Modulated DFB Lasers Having MQW Based Passive Reflectors,” IEEE Photonics Technology Letters 35, 333–336 (2023).
  * [38] H.-T. Cheng, Y.-C. Yang, T.-H. Liu, and C.-H. Wu, “Recent Advances in 850 nm VCSELs for High-Speed Interconnects,” Photonics 9 (2022).
  * [39] M. Hoser, W. Kaiser, D. Quandt, J. Bueno, S. Saintenoy, and S. Eitel, “Highly Reliable 106 Gb/s PAM-4 850 nm Multi-Mode VCSEL for 800G Ethernet Applications,” in _Optical Fiber Communication Conference (OFC) 2022,_ (2022), p. Tu2D.5.
  * [40] D. Wu, X. Yu, H. Wu, W. Fu, and M. Feng, “Single-mode 850nm VCSELs Demonstrate 96 Gb/s PAM4 OM4 Fiber Link for Extended Reach to 1km,” in _Optical Fiber Communication Conference (OFC) 2022,_ (2022), p. W2A.7.
  * [41] J. Wang, Y. Ji, Z. Yang, H. Hu, J. Chen, H. Li, F. Li, S. Li, J. Kapraun, C. Shen, and C. J. Chang-Hasnain, “300-m Multimode Fiber Transmission of 106Gbps PAM-4 Using 850nm High-Contrast-Grating Few-mode VCSELs,” in _European Conference on Optical Communication (ECOC),_ (2022), p. Tu3E.3.
  * [42] C. Ge, L. Dong, X. Gu, and F. Koyama, “1060nm Single-mode Intra-cavity Metal-aperture VCSEL for over 2km Standard 1300nm SMF Transmission,” in _CLEO 2023,_ (2023), p. STh4Q.2.
  * [43] J. Chi, X. Li, C. Niu, and J. Zhao, “Output Optical Power Enhancement of Push-Pull Modulated DFB Laser With Asymmetric Structure,” IEEE Photonics Journal 15, 1–8 (2023).
  * [44] D. Che, Y. Matsui, X. Chen, R. Schatz, and P. Iannone, “400-Gb/s Direct Modulation using a DFB$+$R Laser,” Opt. Lett. 45, 3337–3339 (2020).
  * [45] J. Chi, C. Niu, and J. Zhao, “Parameter Extraction For Quantum Well DFB Lasers Based on 1D Traveling Wave Model,” IEEE Photonics Journal 14, 1–8 (2022).
  * [46] L. Bjerkan, A. Royset, L. Hafskjaer, and D. Myhre, “Measurement of Laser Parameters for Simulation of High-Speed Fiberoptic Systems,” Journal of Lightwave Technology 14, 839–850 (1996).
  * [47] J. Cartledge and R. Srinivasan, “Extraction of DFB Laser Rate Equation Parameters for System Simulation Purposes,” in _Conference Proceedings LEOS’96 9th Annual Meeting IEEE Lasers and Electro-Optics Society,_ vol. 2 (1996), pp. 248–249.
  * [48] I. Tomkos, I. Roudas, R. Hesse, N. Antoniades, A. Boskovic, and R. Vodhanel, “Extraction of Laser Rate Equations Parameters for Representative Simulations of Metropolitan-area Transmission Systems and Networks,” Optics Communications 194, 109–129 (2001).
  * [49] A. Marchisio, E. Ghillino, V. Curri, A. Carena, and P. Bardella, “Particle Swarm Optimization-assisted Approach for the Extraction of VCSEL Model Parameters,” Opt. Lett. 49, 125–128 (2024).
  * [50] S. J. Zhang, N. H. Zhu, E. Y. B. Pun, and P. S. Chung, “Rate-equation-based Circuit Model of High-speed Semiconductor Lasers,” Microwave and Optical Technology Letters 49, 539–542 (2007).
  * [51] J. Butcher, _Numerical Methods for Ordinary Differential Equations_ (Wiley Blackwell, 2016).
  * [52] A. Griewank, “A mathematical view of automatic differentiation,” Acta Numerica 12, 321–398 (2003).
  * [53] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind, “Automatic differentiation in machine learning: a survey,” J. Mach. Learn. Res. 18, 5595–5637 (2017).
  * [54] W. W. Feng and N. H. Zhu, “Analysis of Chirp Characteristics of DFB Lasers and Integrated Laser-modulators,” Optical and Quantum Electronics 36, 1237–1245 (2004).
  * [55] Q. Ding, Q. Yang, J. Li, C. Gu, Z. Li, X. Tan, F. Gao, X. Chen, and Q. Niu, “Efficient and systematic parameter extraction based on rate equations by DFB equivalent circuit model,” Optics Express 31, 40604–40619 (2023).
  * [56] J. Katz, S. Margalit, C. Harder, D. Wilt, and A. Yariv, “The Intrinsic Electrical Equivalent Circuit of a Laser Diode,” IEEE Journal of Quantum Electronics 17, 4–7 (1981).
  * [57] I. Habermayer, “Nonlinear Circuit Model for Semiconductor Lasers,” Optical and Quantum Electronics 13, 461–468 (1981).
  * [58] R. S. Tucker, “Circuit model of double-heterojunction laser below threshold,” IEEE Proceedings I (Solid-State and Electron Devices) 128, 101–106(5) (1981).
  * [59] S. Ghoniemy, L. MacEachern, and S. Mahmoud, “Extended Robust Semiconductor Laser Modeling for Analog Optical Link Simulations,” IEEE Journal of Selected Topics in Quantum Electronics 9, 872–878 (2003).
  * [60] A. Horri, S. Z. Mirmoeini, and R. Faez, “The Noise Equivalent Circuit Model of Quantum-dot Lasers,” Journal of Russian Laser Research 33, 217–226 (2012).
  * [61] E. Mortazy, V. Ahmadi, and M. Moravvej-Farshi, “An Integrated Equivalent Circuit Model for Relative Intensity Noise and Frequency Noise Spectrum of a Multimode Semiconductor Laser,” IEEE Journal of Quantum Electronics 38, 1366–1371 (2002).
  * [62] L. Ljung, _System Identification_ (John Wiley and Sons, Ltd, 2017), chap. 5 - Models for Time-varying and Nonlinear Systems, pp. 140–161.
  * [63] O. Nelles, _Nonlinear System Identification From Classical Approaches to Neural Networks, Fuzzy Models, and Gaussian Processes_ (Springer Cham, 2020), 2nd ed.
  * [64] T. Jakaša, I. Andročec, and P. Sprčić, “Electricity Price Forecasting — ARIMA Model Approach,” in _2011 8th International Conference on the European Energy Market (EEM),_ (2011), pp. 222–225.
  * [65] D. R. Irmawati, R. M. Atok, and Suhartono, “Hybrid Singular Spectrum Analysis-ARIMA Modelling for Direct and Indirect Forecasting of Farmer’s Term of Trade in East Java,” in _2018 International Conference on Information and Communications Technology (ICOIACT),_ (2018), pp. 889–894.
  * [66] S. Zhang and Y. Liu, “Forecasting Model of Total Import and Export based on ARIMA Algorithm Optimized by BP Neural Network,” in _2023 IEEE 3rd International Conference on Data Science and Computer Application (ICDSCA),_ (2023), pp. 1534–1538.
  * [67] L. Ljung, C. Andersson, K. Tiels, and T. B. Schön, “Deep Learning and System Identification,” 21st IFAC World Congress 53, 1175–1181 (2020).
  * [68] J. Schoukens and L. Ljung, “Nonlinear System Identification: A User-Oriented Road Map,” IEEE Control Systems Magazine 39, 28–99 (2019).
  * [69] H.-Y. Chen, C.-C. Wei, C.-Y. Lin, L.-W. Chen, I.-C. Lu, and J. Chen, “Frequency- and Time-domain Nonlinear Distortion Compensation in High-speed OFDM-IMDD LR-PON with High Loss Budget,” Opt. Express 25, 5044–5056 (2017).
  * [70] Y. Xu, L. Huang, W. Jiang, L. Xue, W. Hu, and L. Yi, “Automatic Optimization of Volterra Equalizer With Deep Reinforcement Learning for Intensity-Modulated Direct-Detection Optical Communications,” Journal of Lightwave Technology 40, 5395–5406 (2022).
  * [71] N. Stojanovic, F. Karinou, Z. Qiang, and C. Prodaniuc, “Volterra and Wiener Equalizers for Short-Reach 100G PAM-4 Applications,” Journal of Lightwave Technology 35, 4583–4594 (2017).
  * [72] F. P. Guiomar, S. B. Amado, C. S. Martins, and A. N. Pinto, “Time-Domain Volterra-Based Digital Backpropagation for Coherent Optical Systems,” Journal of Lightwave Technology 33, 3170–3181 (2015).
  * [73] S. Elliott, “8 - Active Control of Nonlinear Systems,” in _Signal Processing for Active Control,_ (Academic Press, 2001), Signal Processing and its Applications, pp. 367–409.
  * [74] K. L. J. Tsimbinos, “Computational Complexity of Volterra based Nonlinear Compensators,” Electronics Letters 32, 852–854(2) (1996).
  * [75] J. Tsimbinos and K. Lever, “The Computational Complexity of Nonlinear Compensators based on the Volterra Inverse,” in _Proceedings of 8th Workshop on Statistical Signal and Array Processing,_ (1996), pp. 387–390.
  * [76] M. J. Korenberg, S. B. Bruder, and P. J. McLlroy, “Exact Orthogonal Kernel Estimation from Finite Data Records: Extending Wiener’s Identification of Nonlinear Systems,” Annals of Biomedical Engineering 16, 201–214 (1988).
  * [77] E. Thomas, J. van Hemmen, and W. Kistler, “Calculation of Volterra Kernels for Solutions of Nonlinear Differential Equations,” Siam Journal on Applied Mathematics 61 (2000).
  * [78] S. Orcioni, M. Pirani, and C. Turchetti, “Advances in Lee-Schetzen Method for Volterra Filter Identification,” Multidimensional Systems and Signal Processing 16, 265–284 (2005).
  * [79] A. Wills, T. B. Schön, L. Ljung, and B. Ninness, “Identification of Hammerstein–Wiener Models,” Automatica 49, 70–81 (2013).
  * [80] X. Liu, Y. Fan, Y. Zhang, M. Cai, L. Liu, L. Yi, W. Hu, and Q. Zhuge, “Fusing Physics to Fiber Nonlinearity Model for Optical Networks Based on Physics-Guided Neural Networks,” Journal of Lightwave Technology 40, 5793–5802 (2022).
  * [81] Y. Song, Y. Zhang, C. Zhang, J. Li, M. Zhang, and D. Wang, “PINN for Power Evolution Prediction and Raman Gain Spectrum Identification in C$+$L-Band Transmission System,” in _Optical Fiber Communication Conference (OFC) 2023,_ (2023), p. Th1F.5.
  * [82] U. C. de Moura, D. Zibar, A. Margareth Rosa Brusin, A. Carena, and F. Da Ros, “Fiber-Agnostic Machine Learning-Based Raman Amplifier Models,” Journal of Lightwave Technology 41, 83–95 (2023).
  * [83] D. Wang, Y. Song, J. Li, J. Qin, T. Yang, M. Zhang, X. Chen, and A. C. Boucouvalas, “Data-driven Optical Fiber Channel Modeling: A Deep Learning Approach,” Journal of Lightwave Technology 38, 4730–4743 (2020).
  * [84] A. Waibel, T. Hanazawa, G. Hinton, K. Shikano, and K. J. Lang, “Phoneme Recognition Using Time-Delay Neural Networks,” Readings in Speech Recognition pp. 393–404 (1990).
  * [85] A. G. Reza and J.-K. K. Rhee, “Nonlinear Equalizer Based on Neural Networks for PAM-4 Signal Transmission Using DML,” IEEE Photonics Technology Letters 30, 1416–1419 (2018).
  * [86] Y. LeCun and Y. Bengio, _Convolutional Networks for Images, Speech, and Time Series_ (MIT Press, 1998), p. 255–258.
  * [87] K. Cho, B. van Merrienboer, C. Gulcehre, F. Bougares, H. Schwenk, and Y. Bengio, “Learning phrase representations using RNN encoder-decoder for statistical machine translation,” in _Conference on Empirical Methods in Natural Language Processing (EMNLP 2014),_ (2014).
  * [88] S. Hochreiter and J. Schmidhuber, “Long Short-Term Memory,” Neural Computation 9, 1735–1780 (1997).
  * [89] T. K. Rusch and S. Mishra, “UnICORNN: A recurrent model for learning very long time dependencies,” in _Proceedings of the 38th International Conference on Machine Learning,_ vol. 139 M. Meila and T. Zhang, eds. (2021), pp. 9168–9178.
  * [90] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. u. Kaiser, and I. Polosukhin, “Attention is All you Need,” in _Advances in Neural Information Processing Systems,_ vol. 30 I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, eds. (2017).
  * [91] N. Zhang, H. Yang, Z. Niu, L. Zheng, C. Chen, S. Xiao, and L. Yi, “Transformer-Based Long Distance Fiber Channel Modeling for Optical OFDM Systems,” Journal of Lightwave Technology 40, 7779–7789 (2022).
  * [92] Y. Zhu, J. Ye, L. Yan, T. Zhou, P. Li, X. Zou, and W. Pan, “Transformer-Based High-Fidelity Modeling Method for Radio Over Fiber Link,” Journal of Lightwave Technology 41, 2657–2665 (2023).
  * [93] B. Li, Q. Du, T. Zhou, Y. Jing, S. Zhou, X. Zeng, T. Xiao, J. Zhu, X. Liu, and M. Zhang, “ODE Transformer: An Ordinary Differential Equation-Inspired Model for Sequence Generation,” in _Proceedings of the 60th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers),_ S. Muresan, P. Nakov, and A. Villavicencio, eds. (2022), pp. 8335–8351.
  * [94] A. Gatto, M. Rapisarda, P. Parolari, and P. Boffi, “Discrete Multitone Modulation for Short-Reach Mode Division Multiplexing Transmission,” Journal of Lightwave Technology 37, 5185–5192 (2019).
  * [95] G. Böcherer, P. Schulte, and F. Steiner, “Probabilistic Shaping and Forward Error Correction for Fiber-Optic Communication Systems,” Journal of Lightwave Technology 37, 230–244 (2019).
  * [96] D. Che and W. Shieh, “Approaching the Capacity of Colored-SNR Optical Channels by Multicarrier Entropy Loading,” Journal of Lightwave Technology 36, 68–78 (2018).
  * [97] D. Kim and H. Kim, “Capacity-achieving Symbol Distributions for Directly Modulated Laser and Direct Detection Systems,” Optics Express 31, 12609–12623 (2023).
  * [98] N.-P. Diamantopoulos, H. Yamazaki, S. Yamaoka, M. Nagatani, H. Nishi, H. Tanobe, R. Nakao, T. Fujii, K. Takeda, T. Kakitsuka, H. Wakita, M. Ida, H. Nosaka, F. Koyama, Y. Miyamoto, and S. Matsuo, “$>$100-GHz Bandwidth Directly-Modulated Lasers and Adaptive Entropy Loading for Energy-Efficient $>$300-Gbps/$\lambda$ IM/DD Systems,” Journal of Lightwave Technology 39, 771–778 (2021).
  * [99] C. Kottke, C. Caspar, V. Jungnickel, R. Freund, M. Agustin, and N. N. Ledentsov, “High Speed 160 Gb/s DMT VCSEL Transmission Using Pre-equalization,” in _Optical Fiber Communication Conference,_ (Optica Publishing Group, 2017), p. W4I.7.
  * [100] D. Dasalukunte, F. Rusek, and V. Owall, “Multicarrier Faster-Than-Nyquist Transceivers: Hardware Architecture and Performance Analysis,” IEEE Transactions on Circuits and Systems I: Regular Papers 58, 827–838 (2011).
  * [101] N. Bamiedakis, D. G. Cunningham, and R. V. Penty, “Linearisation Method of DML-Based Transmitters for Optical Communications Part I: Theory and Simulation Studies,” Journal of Lightwave Technology 39, 5815–5827 (2021).
  * [102] Q.-A. Ding, Huixin-Liu, X. Cheng, X.-J. Wang, W. Huang, L. Du, G. Zhu, L. Zheng, and Q. Yang, “Suppressing DFB Lasers Relaxation Oscillations by High-Order Continuous Shaping Current With CNN Prediction,” IEEE Journal of Quantum Electronics 59, 1–11 (2023).
  * [103] Z. Feng, N. Li, W. Li, P. He, M. Luo, Q. Hu, L. Huang, and Y. Jiang, “A Simplified Volterra Equalizer Based on System Characteristics for Direct Modulation Laser (DML)-Based Intensity Modulation and Direct Detection (IM/DD) Transmission Systems,” Photonics 10 (2023).
  * [104] N.-P. Diamantopoulos, H. Nishi, W. Kobayashi, K. Takeda, T. Kakitsuka, and S. Matsuo, “On the Complexity Reduction of the Second-Order Volterra Nonlinear Equalizer for IM/DD Systems,” Journal of Lightwave Technology 37, 1214–1224 (2019).
  * [105] Z. Xu, C. Sun, T. Ji, H. Ji, and W. Shieh, “Cascade Recurrent Neural Network Enabled 100-Gb/s PAM4 Short-Reach Optical Link Based on DML,” in _Optical Fiber Communication Conference (OFC) 2020,_ (2020), p. W2A.45.
  * [106] F. Tian and C. Yang, “Deep Belief Network-Hidden Markov Model-based Nonlinear Equalizer for VCSEL based Optical Interconnect,” Science China Information Sciences 63 (2020).
  * [107] F. A. Aoudia and J. Hoydis, “End-to-End Learning of Communications Systems Without a Channel Model,” in _2018 52nd Asilomar Conference on Signals, Systems, and Computers,_ (2018), pp. 298–303.
  * [108] M. Srinivasan, J. Song, C. Häger, K. Szczerba, H. Wymeersch, and J. Schröder, “Learning Optimal PAM Levels for VCSEL-based Optical Interconnects,” in _European Conference on Optical Communication (ECOC) 2022,_ (2022), p. We2C.3.
  * [109] L. Minelli, F. Forghieri, T. Shao, A. Shahpari, and R. Gaudino, “TDECQ-Based Optimization of Nonlinear Digital Pre-Distorters for VCSEL-MMF Optical Links Using End-to-End Learning,” J. Lightwave Technol. 42, 621–635 (2024).
  * [110] B. Karanov, M. Chagnon, F. Thouin, T. A. Eriksson, H. Bülow, D. Lavery, P. Bayvel, and L. Schmalen, “End-to-End Deep Learning of Optical Fiber Communications,” Journal of Lightwave Technology 36, 4843–4855 (2018).
  * [111] R. T. Jones, T. A. Eriksson, M. P. Yankov, and D. Zibar, “Deep Learning of Geometric Constellation Shaping Including Fiber Nonlinearities,” in _2018 European Conference on Optical Communication (ECOC),_ (2018).
  * [112] P. Freire, S. Srivallapanondh, B. Spinnler, A. Napoli, N. Costa, J. E. Prilepsky, and S. K. Turitsyn, “Computational Complexity Optimization of Neural Network-Based Equalizers in Digital Signal Processing: A Comprehensive Approach,” Journal of Lightwave Technology pp. 1–25 (2024).
  * [113] J. Cartledge and G. Burley, “The effect of laser chirping on lightwave system performance,” Journal of Lightwave Technology 7, 568–573 (1989).
  * [114] J. Cartledge and R. Srinivasan, “Extraction of dfb laser rate equation parameters for system simulation purposes,” Journal of Lightwave Technology 15, 852–860 (1997).