# Improving LSTM Neural Networks for Better Short-Term Wind Power Predictions

Maximilian Du  Fayetteville-Manlius High School
Manlius, NY, United States
<EMAIL_ADDRESS>
###### Abstract

This paper improves wind power prediction via weather forecast-contextualized
Long Short-Term Memory Neural Network (LSTM) models. Initially, only wind
power data was fed to a generic LSTM, but this model performed poorly, with
erratic and naive behavior observed on even low-variance data sections. To
address this issue, weather forecast data was added to better contextualize
the power data, and LSTM modifications were made to address specific model
shortcomings. These models were tested through both a Normalized Mean Absolute
Error and the Naive Ratio (NR), which is a score introduced by this paper to
quantify the unwanted presence of naive character in trained models. Results
showed an increased accuracy with the addition of weather forecast data on the
modified models, as well as a decrease in naive character. Key contributions
include making improved LSTM variants, usage of weather forecast data, and the
introduction of a new model performance index.

###### Index Terms:

Wind Power Prediction, LSTM, Renewable Energy Integration, Persistence
Algorithm Quantification

##  I Introduction

Wind energy is both a clean and renewable source of energy, and its usage
across the world is increasing. However, because the wind is uncontrollable,
wind turbine farms can cause grid instability, and the cost of cycling
conventional power plants to compensate for this variability is around 157
million dollars per year in areas with only 35% integration [1]. In the
future, this instability and subsequent costs will only worsen as wind power
penetration increase and conventional power capacity decrease [2].

Predicting wind power output reduces these adverse effects of wind power
penetration by playing important roles in load balancing [3] and reserve
optimization[4]. Furthermore, as conventional capacity decreases, more
accurate predictions over ultra-short intervals of time (5-10 minutes) are
needed to maintain power grid stability [5].

There are many techniques available for time series problems, e.g. wind power
prediction, such as mathematical and statistical modeling. However, machine
learning approaches can be superior in terms of robustness and accuracy
because they are very adaptable and do not rely on simulating a turbine
environment [6]. One such machine learning model is the Long Short-Term Memory
Neural Network (LSTM). The LSTM is a variant of the Recurrent Neural Network
(RNN) that uses a continuous cell state to carry temporal information. Such a
cell state is modified by the Forget, Input/Input Transform, and Output gates
to forget, add, and recall information, respectively. The LSTM was shown to be
powerful at time-series prediction, and once trained, they are relatively
lightweight [7].

Existing research has already approached the problem of wind power prediction
using LSTMs, e.g. Xiaoyun et al [6]. However, while such research considered
the importance of _current weather data_ from the Numerical Weather Prediction
(NWP) model as contextualization, the usage of weather _forecast data_ also
provided by the NWP remains as an additional and less-explored option.

This paper improves the accuracy of the LSTM model by using NWP forecast data.
Furthermore, this paper proposes a new way of numerically quantifying the
presence of Naive Character, an undesirable inclination to directly use
current states as predictions for future states. To decrease such character
and increase prediction accuracy, multiple new LSTM modifications were
designed, tested, and compared using a Normalized Mean Absolute Error (NMAE)
as well as a newly proposed Naive Ratio.

##  II Motivational Example: Problems Encountered with the Generic LSTM

As a starting point, unmodified (generic) LSTMs were trained and evaluated on
power data using the Tensorflow Python library. Although LSTMs have been
established as reliable models, there are several problems observed with the
trained model.

###  II-A Naivety and Step Plots

Naive character, or the _Persistence Algorithm_ , is the undesirable usage of
the current state as the prediction for the next state. This is mathematically
represented as $x^{\prime}_{t+1}=x_{t}$. On a step plot, which holds each
point on the plot constant until the next point, naive character is observed
as a one-step shift between truth and prediction. Other works such as Chang et
al [8] and Xiaoyun et al [6] have generated prediction vs. truth line plots,
but these line plots do not visibly show naive character, especially when
viewed over a long time range relative to their timestep.

Figure 1: Step plot of Generic LSTM output shows Naive behavior, which
indicates that this model has not learned any real trends.

In Figure 1, a step plot, this offset between the truth and prediction plots
is apparent. The presence of this naive behavior on the generic LSTM shows
that the trained model has found the persistence algorithm as a local minimum
instead of truly modeling wind power trends.

###  II-B Predication Oscillation and Negative Outputs

Another problem observed through step plots and cell state distributions was
prediction oscillation caused by cell state divergence.

Figure 2: A different step plot section of Generic LSTM output shows
prediction oscillation, which indicates that this model is unstable.

As can be seen in Figure 2, which shows a different section of the test set,
the prediction values oscillate from point to point. This problem can be
justified by the magnitudes of the cell states, whose outliers diverged up to
$\pm 10^{5}$ at the end of training. With this range, optimization favored
very small output weights to produce predictions of similar magnitude as the
truth. These diminishing weights attenuate more common mid-range cell values,
leading to under-sensitivity to specific states and over-sensitivity to
others, resulting in unstable and inconsistent outputs.

Negative outputs, which are undesirable for power predictions, were also
observed in this model due to the lack of output regulation.

##  III Data Acquisition and Processing

Data from two open-source databases was used for training and testing the
models.

###  III-A Wind Power Data

This database from the National Renewable Energy Laboratories (NREL) [9] was
chosen as the source of power data. For this paper, data from a 16 MW capacity
turbine farm off the coast of Long Island, New York, United States was used
due to its relatively greater power fluctuations.

###  III-B Weather Forecast Data

It was hypothesized that weather forecast data would enable the model to have
some future context while harnessing the power of existing weather models.
This forecast data was obtained from the Rapid Refresh (RAP) database from
NOAA [10]. While there existed a large amount of weather information (317
parameters over 18 hours of forecast), to prevent overloading the LSTM with
unnecessary amounts of data, specific weather parameters across an empirically
chosen 2-hour forecast interval were used. These included the following:

  * •

Surface air pressure

  * •

Ground temperature

  * •

2-meter air temperature

  * •

Relative humidity

  * •

Wind gust peak velocities

All of the above parameters were chosen based on their ability to represent
current and future weather conditions relevant to wind power generation.

###  III-C Data Concatenation and Implementation

The NREL wind power set was found in five-minute increments, while the RAP
dataset was found in one-hour increments. To combine them, the RAP dataset
values were kept constant in blocks of 12 wind power points.

The addition of weather forecast data to the power data meant increasing the
input vector size. As a starting point, a direct feed-in was used, where all
input values were made into input nodes. This was done to determine the
effectiveness of the forecast parameters in improving model performance.

Compression algorithms like Principal Component Analysis (PCA), used in
Xiaoyun et al [6], can reduce noise and complexity, so a model using PCA was
made to gauge performance impacts. It was observed that the first principle
component contained more than 98% of the variance, so it was the only
component used.

##  IV Proposed Structure Modifications

With the addition of weather data and the usage of PCA, there was a
significant increase in the prediction accuracy of the generic LSTM (Table I).
However, the forecast data did not address prediction oscillation, and naive
character was still observed, which indicated sub-optimal trend learning.
Thus, some structural modifications were needed to address specific problems
of the Generic LSTM and improve performance.

###  IV-A Main Modification: Cell Augmentation and Hyperbolic Tangent Addition

Figure 3: Model of LSTM with added hyperbolic tangent and cell augmentation,
which greatly increased accuracy.

To address prediction oscillation, a hyperbolic tangent activation function
(shown as an orange oval in Figure 3) was added to the cell propagation, which
keeps values within the interval $(-1,1)$. This removed the cell divergence
phenomenon found in the generic model, and oscillations were no longer
observed.

To further decrease naive character and to account for the increased amount of
information, the cell state dimension of the LSTM was made independent of
input size, allowing for the expansion of state sizes to keep more temporal
information. To create this cell dimension independence, a fully-connected
single-layer neural network with a rectified linear unit (ReLU) activation
function was added after output gating (shown as an orange square in Figure
3). This activation function eliminated negative outputs.

Initial testing showed that these modifications significantly increased
prediction accuracy (Table III), so the structure was kept as a backbone
model, denoted mLSTM, or the modified LSTM.

Because the addition of the hyperbolic tangent re-introduced the vanishing
gradient problem that prevents effective training on continuous data, _hybrid
backpropagation_ was used. This method uses discrete blocks of sequential data
for model training while initializing cell and hidden states with those from
the previous iteration, allowing for some long-term trend learning without the
use of an extended backpropagation.

###  IV-B Combined Input and Forget Gates (CIFG) with Peephole Connections

The CIFG modification, explored in Greff et al [11], uses a single gate to
modify both the Input Transform and Forget values, setting $f_{t}=1-i_{t}$.
Thus, the more attenuated a cell value becomes after forget gating, the more
it is modified by the input gate. This increases memory efficiency by
preventing excessive null (zero) values in the cell state.

The peephole, explored in Gers et al [12], gives the gates information about
the current cell state and allows the LSTM to learn precise timings by
enabling two-way data flow, which is beneficial for long-term trends and can
reduce naive behavior.

These modifications were proposed before, but there exists limited literature
on their usage together and none on the new _mLSTM_ , yet it has been reported
that both modifications have a positive impact on LSTM models.

###  IV-C Compression Layer

Data compression in the form of PCA was previously used on the weather data to
decrease input complexity. In order to gauge the effects of a trainable
compression algorithm on the model, a compression layer was proposed for the
_mLSTM_. This modification takes the raw input vector and compresses it
through a fully-connected layer. The construction is similar to that of an
autoencoder, but it is trained directly with the model. Thus, the compression
layer can learn to isolate application-relevant data components, which may or
may not be the components that are the most important for data preservation.
To keep results comparable with the PCA algorithm, which only used one
principle value, the compression layer was also trained to compress down to a
single scalar value.

##  V Cross-Model Comparison

###  V-A Hyperparameter Optimization

The _mLSTM_ and its variants had three hyperparameters that were adjusted:

  * •

Learning rate

  * •

Cell dimension

  * •

Instance iteration length

The learning rate tunes the speed of gradient descent, the cell dimension sets
the size of the cell and hidden states, and the instance iteration length
determines the number of past data points fed to the model per iteration. A
higher learning rate results in faster convergence, but it also increases the
chance of sub-optimal endpoints. Likewise, a larger cell dimension increases
the ability to store more information, but past a certain threshold, the extra
space will only serve to increase complexity. Lastly, more data points given
to the model leads to greater long-term trend recognition but it also
increases the effects of the vanishing gradient problem, even with hybrid
backpropagation.

There are no guarantees that each of the LSTM modifications shares the same
set of optimal hyperparameters, so a genetic hyperparameter optimization was
run on each modification. The genetic algorithm created 10 child models and
then partially trained each child to 3000 epochs. Each model was then
evaluated on the same test set, and the hyperparameters of the best two
children were kept. Ten new children were made from the two vectors, using
crossover and mutation to increase variance and decrease the chance of getting
stuck in a local minimum. This was done until stabilization of all values was
achieved, which was empirically set at 20 generations for all models.

###  V-B Training, Validation, and Testing

Each model, including the generic LSTM, _mLSTM_ , and its derivatives, was
trained for $160000$ epochs on the same partition ($80\%$) of the dataset made
from the combination of NREL and RAP data. Due to a large number of data
points ($84096$), a training session would cycle through the training set less
than 15 times, minimizing model overfitting. The validation set was also the
same across all models, and it served as a quantifier of training progression.
To ensure a fair comparison, the same test set was used across all models,
which consisted of $1000$ data points from a non-trained partition of the
dataset.

In order to measure model accuracy, the Normalized Mean Absolute Error (NMAE)
was used. The NMAE and its variant, the Mean Absolute Error (MAE), is defined
as follows, where $MAX(x)$ represents the wind farm capacity of 16 MW:

$e_{MAE}(x^{\prime},x)=\frac{1}{N}*\sum_{i=1}^{N}|x^{\prime}(i)-x(i)|$ (1)
$e_{NMAE}(x^{\prime},x)=\frac{100\%}{MAX(x)}*\frac{1}{N}*\sum_{i=1}^{N}|x^{\prime}(i)-x(i)|$
(2)

The NMAE allows test results to be better extrapolated to farms with different
capacities, while the MAE gives a better quantification of true power error.
Due to the non-convex nature of LSTM optimization, the initial weight states
proved to be non-trivial in the final performance of the trained models. Thus,
each model was trained 5 times and a margin of error (MOE) at $95\%$ was
calculated.

###  V-C Naive Ratio

Throughout experimentation, it was observed that some model modifications had
accuracies with non-statistically significant differences, yet step plots of
their test sets revealed a wide range of naive character. Thus, to better
compare models, a metric of naive character was created, which quantifies the
step plot observations. Under a pure persistence algorithm, the model
predictions would be one timestep behind the truth, meaning that shifting the
prediction backward one timestep would decrease the loss to zero. In the
opposite case, where the model is perfectly accurate, shifting the prediction
backward one timestep would increase the loss to a non-zero value. The models
were neither of the boundary cases, but comparing a loss where the predictions
are not time-shifted to a loss where the predictions are shifted one time-step
into the past can give insight into the magnitude of naive behavior exhibited
by the trained model.

The back-shifted loss function is defined as follows, where subscript $t+1$
represents the back shift (the NMAE subscript was removed for conciseness):

$e_{t+1}(x^{\prime},x)=\frac{1}{N}*\sum_{i=1}^{N-1}|x^{\prime}(i+1)-x(i)|$ (3)

If $e_{t+1}(x^{\prime},x)<e_{t}(x^{\prime},x)$ then the model is defined to
have _Naive Character_ , and the presence of such character is placed on a
continuous spectrum through the proposed _Naive Ratio_ (NR), which is defined
as the following:

$NR=\frac{e_{t}(x^{\prime},x)}{e_{t+1}(x^{\prime},x)},NR\in(0,\infty)$ (4)

A Naive Ratio score of anything greater than $1$ means that the model has some
naive character, and a score of $\infty$ indicates a perfectly naive model. A
score of anything less than $1$ indicates that naive character takes a
minority contribution to the model’s behavior, and a perfectly accurate model
would have a Naive Ratio score of $0$.

###  V-D Control Models

There were two control models used in experimentation that set benchmark
values for the _mLSTM_ and its variants. As a baseline machine learning model,
a generic LSTM with PCA was used, and as a baseline naive model, a persistence
algorithm was used.

##  VI Results and Analysis

###  VI-A Impact of Weather Data on Performance

TABLE I: Effect of Weather Data on Generic LSTM Model |  Power Data Only |  With Weather |  With Weather and PCA
---|---|---|---
NMAE (%) |  2.409 |  1.649 |  1.130
95% MOE |  0.450 |  0.231 |  0.212

The forecast data worked as intended, with the forecast-supplemented model
having a statistically significant improvement over the direct power model,
and a PCA compression also having a statistically significant improvement over
an uncompressed model (Table I).

TABLE II: Effect of Data Type on Naive Ratio Score (Generic) |  Weather |  Weather With PCA |  Power Data Only
---|---|---|---
Naive Ratio Score |  1.050 |  1.353 |  1.145

In Table II, it is observed that the addition of uncompressed weather forecast
data to the Generic LSTM models resulted in the decrease of naive character.
Combined with the previously observed increase in accuracy, it can be
concluded that weather forecast data is providing adequate context for wind
power.

###  VI-B Model Modification Performance

TABLE III: Evaluation of Model Performance |  Naive Ctrl. |  LSTM w/PCA Ctrl. |  mLSTM Direct Feed |  mLSTM PCA |  mLSTM Compr. |  mLSTM CFIG + Peep
---|---|---|---|---|---|---
NMAE (%) |  1.451 |  1.130 |  0.522 |  0.472 |  0.438 |  0.590
95% MOE |  0.262 |  0.212 |  0.038 |  0.038 |  0.031 |  0.038
Naive Ratio |  $\infty$ |  1.353 |  0.774 |  0.876 |  0.677 |  0.939

Shown in Table III is a compilation of test results on the new models and
controls. It is seen that the _mLSTM_ -based models had higher accuracies and
lower Naive Ratio scores. This increase in performance is attributed to the
bounding of cell state values and the expansion of cell state sizes, which
keep the weights at a reasonable level and allow for greater storage of past
information, respectively.

Among the variants, the compression layer addition showed the highest accuracy
and lowest naive character. This improvement is clearly shown in Figure 4,
which compares this version to the power-data-only generic LSTM on the same
section of the test set. This is consistent with the observation made in Table
I that compression increases accuracy, which is caused by the lowered
computational complexity and the isolation of trend-relevant data from noise.
However, in observing the differences between the compression layer and PCA on
the _mLSTM_ , the Naive Ratio score must be considered. The _mLSTM_ with PCA
showed an increase in naive character from the unmodified _mLSTM_ , while the
compression layer showed the opposite. Because an increase in naive character
corresponds to a decrease in trend learning, it can be reasonably concluded
that the PCA, while being able to capture data variation, was not able to
optimally preserve trend-relevant data. However, the compression layer,
trained with the model, was able to adapt its compression and learned to keep
such trend-relevant data while still reducing complexity. The differences in
naive character between the PCA and compression _mLSTM_ models shows that
optimizing for variation may attenuate more relevant yet subtler data. The
increased naive character with PCA is also observed on the Generic LSTM (Table
II). Taken holistically, it can be concluded that an adaptive compression
method was a better feature encoding method for this application.

It is worth noting that the NMAE of the persistence algorithm ($1.451\%$,
Table III) was lower than that of a generic LSTM on uncompressed forecast data
($1.649\%$, Table I). This shows that pure naive behavior has a deceptively
low NMAE, which justifies its presence in some trained models. Thus, the naive
ratio score is very important in disseminating false naive accuracy from real
model accuracy.

The results also show that combining CFIG and peephole on the _mLSTM_
decreased performance, both in terms of NMAE and Naive Ratio. However, the
addition of the CFIG allowed for a more uniform forget gate distribution after
training, which indicates a stable convergence. The worse performance can be
attributed to the Peephole addition. In an unlisted trial, a Peephole-only
_mLSTM_ had an NMAE of $2.328\%$ and a Naive Ratio of $1.040$, both of which
are significantly worse than an unmodified _mLSTM_. Peepholes connect the cell
states with the gates, which can cause unwanted gradient complications like
local minima, which inhibits optimization. The addition of the CFIG
modification decreased gate complexity and counteracted some negative effects
of the Peephole, but the impact of the Peephole was still present in the final
model.

Figure 4: By adding forecast data and making some LSTM structural
modifications, the naive character and prediction loss greatly decrease.

##  VII Conclusion, Limitations, and Future Development

It was shown that using weather forecast data to supplement the power data is
effective in increasing prediction accuracy, with an adaptive compression
algorithm, cell regulation, and ReLU units serving to further increase
accuracy and decrease naive character. The proposed _mLSTM_ models are able to
outperform traditional LSTM models at a NMAE of $0.438\%$ (Table III) with the
compression modification, or approximately 75 kW on a 16 MW farm. This
improvement is also qualitatively observed in Figure 4. Such improvements,
paired with future model developments, can be one step towards a more
renewable energy-integrated future.

However, it is worth noting that such findings are not guaranteed to
generalize to all wind farms, whose behaviors vary with location and climate.
As such, future developments include trials at multiple locations to measure
the location consistency of the new models. The usage of region-based
predictions will also be explored, which contain less noise than single-farm
predictions.

##  VIII Acknowledgements

Special thanks to Zhenhua Liu and Joshua Comden of Stony Brook University for
their mentorship and proofreading.

## References

  * [1] D Lew, G Brinkman, and E Ibanez.  Western wind and solar integration study phase 2.  2013\.
  * [2] American Wind Energy Association .  20% wind energy by 2030.  2009\.
  * [3] B. Khorramdel, C Chung, N Safari, and G Price.  A fuzzy adaptive probabilistic wind power prediction framework using diffusion kernel density estimators.  IEEE Trans. Power Syst., page 958–965, 2018.
  * [4] AY Wang, W Liao, and Y Chang.  Gated recurrent unit network-based short-term photovoltaic forecasting.  2018\.
  * [5] G.H. Riahy and M Abedi.  Short term wind speed forecasting for wind turbine applications using linear prediction method.  Renewable Energy, 33:35–41, 2008.
  * [6] Qu Xiaoyun, Kang Xiaoning, Zhang Chao, Jian Shuai, and Ma Xiuda.  Short-term prediction of wind power based on deep long short-term memory.  IEEE PES Asia-Pacific Power and Energy Conference, pages 1148–1152, 2016.
  * [7] S. Hochreiter and J. Schmidhuber.  Long short-term memory.  Neural Computation, 9:1735–1780, 1997.
  * [8] G Chang, J Lu, Y Chang, and Y Lee.  An improved neural network-based approach for short-term wind speed and power forecast.  Elsevier: Renewable Energy, 2017.
  * [9] NREL wind prospector database.  maps.nrel.gov/wind-prospector, 2010.
  * [10] Rapid refresh (RAP) weather database.  www.ncdc.noaa.gov/data-access/model-data/model-datasets/rapid-refresh-rap, 2010\.
  * [11] Klaus Greff, Rupesh Srivastava, Jan Koutnik, Bas Steunebrink, and Jurgen Schmidhuber.  Lstm: A search space odyssey.  2017\.
  * [12] Felix Gers and Jurgen Schmidhuber.  Recurrent nets that time and count.  2000\.


# Supplemental Material: Tunable Magnonic Chern Bands and Chiral Spin Currents
in Magnetic Multilayers

Zhongqiang Hu<EMAIL_ADDRESS>Department of Electrical Engineering and
Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139,
USA  Liang Fu  Department of Physics, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA  Luqiao Liu  Department of Electrical Engineering and
Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139,
USA

##  Appendix A  Interlayer dipolar field from propagating magnons

We first discuss the magnetic bilayers shown in Fig. S1(a), which consist of
two different ferromagnetic layers with in-plane magnetic moments. For
convenience, we set up a global coordinate system $(x,y,z)$, and denote the
angle between the equilibrium moment orientation and the $x$ axis as
$\theta_{j}$ [Fig. S1(b)]. Here, $j=1,2$ is the layer index. In the discussion
throughout the paper including the main text and Supplemental Material, we use
the dimensionless reduced moment $\bm{m}_{j}=\bm{M}_{j}/M_{{\rm s}j}$, where
$M_{{\rm s}j}$ is the saturated magnetization. Within each layer, we also set
up a local coordinate system ($x_{j},y_{j},z$), with the equilibrium moment
orientation chosen as the $y_{j}$ axis.

Figure 1: (a) Magnetic bilayers consisting of two different ferromagnetic
layers with in-plane magnetic moments. (b) The global and local coordinate
systems, as seen from the top view of the bilayers. The equilibrium moment
orientations in two layers have different angles with the $x$ axis. (c) When
calculating the interlayer dipolar energy, $\bm{r}$ and $\bm{r}^{\prime}$ are
in the different layers. (d) When calculating the intralayer dipolar energy,
$\bm{r}$ and $\bm{r}^{\prime}$ are in the same layer.

The dipolar energy between a pair of moments $\bm{m}_{1}$ (located at
$\bm{r}$) and $\bm{m}_{2}$ (located at $\bm{r}^{\prime}$) is

$\delta E_{\rm d}=-\frac{\mu_{0}M_{{\rm s1}}M_{{\rm
s2}}}{4\pi}\frac{3(\bm{m}_{1}\cdot\bm{R})(\bm{m}_{2}\cdot\bm{R})-(\bm{m}_{1}\cdot\bm{m}_{2})R^{2}}{R^{5}},$
(1)

where $\bm{R}=\bm{r}-\bm{r}^{\prime}=(X,Y,Z)$. If defining a matrix

$\begin{aligned} \bm{F}=-\frac{\mu_{0}M_{{\rm s1}}M_{{\rm s2}}}{4\pi
R^{5}}\begin{bmatrix}3X^{2}-R^{2}&3XY&3XZ\\\ 3XY&3Y^{2}-R^{2}&3YZ\\\
3XZ&3YZ&3Z^{2}-R^{2}\end{bmatrix}\end{aligned},$ (2)

we can write Eq. (1) as $\delta E_{\rm d}=\bm{m}_{2}^{T}\bm{F}\bm{m}_{1}$. We
note that here $\bm{m}_{j}=[m_{jx},m_{jy},m_{jz}]^{T}$ is defined as a column
vector. As shown in Fig. S1(c), the dipolar energy between two layers (i.e., a
pair of continuous bodies) is

$\begin{aligned} E_{{\rm
d}}=\iint\bm{m}_{2}^{T}(\bm{r^{\prime}})\bm{F}(\bm{r-r^{\prime}})\bm{m}_{1}(\bm{r})d^{3}\bm{r}d^{3}\bm{r^{\prime}}\end{aligned}.$
(3)

After the Fourier transformation
$\tilde{\bm{m}}_{j}=\mathcal{F}\\{\bm{m}_{j}\\}$ and
$\tilde{\bm{F}}=\mathcal{F}\\{\bm{F}\\}$, it can be written as

$\begin{aligned} E_{{\rm
d}}=\int\tilde{\bm{m}}_{2}^{\dagger}(\bm{k})\tilde{\bm{F}}(\bm{k})\tilde{\bm{m}}_{1}(\bm{k})\frac{d^{3}\bm{k}}{(2\pi)^{3}}\end{aligned}.$
(4)

According to [1], $\tilde{\bm{F}}(\bm{k})$ is approximated as

$\begin{aligned} \tilde{\bm{F}}(\bm{k})\approx-\frac{\mu_{0}M_{{\rm
s1}}M_{{\rm
s2}}}{3k^{2}}\begin{bmatrix}3k_{x}^{2}-k^{2}&3k_{x}k_{y}&3k_{x}k_{z}\\\
3k_{x}k_{y}&3k_{y}^{2}-k^{2}&3k_{y}k_{z}\\\
3k_{x}k_{z}&3k_{y}k_{z}&3k_{z}^{2}-k^{2}\end{bmatrix}\end{aligned},$ (5)

when $k=|\bm{k}|$ is not large compared with the Brillion zone boundary.
Considering the thin-film configuration with the lowest standing mode along
the $z$ direction, i.e.,
$\bm{m}_{j}(\bm{r})=\bm{m}_{j}(\bm{r}_{\parallel})S_{0}(z;d_{j})$, where
$d_{j}$ is the layer thickness and

$\begin{gathered}S_{0}(z;d_{1})=\left\\{\begin{aligned} 1\quad&(-d_{1}<z<0)\\\
0\quad&({\rm otherwise}),\end{aligned}\right.\\\
S_{0}(z;d_{2})=\left\\{\begin{aligned} 1\quad&(0<z<d_{2})\\\ 0\quad&({\rm
otherwise}),\end{aligned}\right.\end{gathered}$ (6)

we can rewrite Eq. (4) as

$\displaystyle E_{{\rm
d}}=\int\tilde{\bm{m}}_{2}^{\dagger}(\bm{k}_{\parallel})\tilde{S}_{0}^{*}(k_{z};d_{2})\tilde{\bm{F}}(\bm{k})\tilde{S}_{0}(k_{z};d_{1})\tilde{\bm{m}}_{1}(\bm{k}_{\parallel})\frac{d^{3}\bm{k}}{(2\pi)^{3}}=\int\tilde{\bm{m}}_{2}^{\dagger}(\bm{k}_{\parallel})\bm{G}(\bm{k}_{\parallel})\tilde{\bm{m}}_{1}(\bm{k}_{\parallel})\frac{d^{2}\bm{k}_{\parallel}}{(2\pi)^{2}},$
(7)

where

$\begin{gathered}\tilde{S}_{0}(k_{z};d_{1})=\mathcal{F}\\{S_{0}(z;d_{1})\\}=d_{1}{e^{ik_{z}d_{1}/2}}{\rm
sinc}(k_{z}d_{1}/2),\\\
\tilde{S}_{0}(k_{z};d_{2})=\mathcal{F}\\{S_{0}(z;d_{2})\\}=d_{2}e^{-ik_{z}d_{2}/2}{\rm
sinc}(k_{z}d_{2}/2),\\\ \end{gathered}$ (8)

and

$\displaystyle\bm{G}({\bm{k}_{\parallel}})=\int\tilde{S}_{0}^{*}(k_{z};d_{2})\tilde{\bm{F}}(\bm{k})\tilde{S}_{0}(k_{z};d_{1})\frac{dk_{z}}{2\pi}=\frac{\mu_{0}M_{{\rm
s1}}M_{{\rm
s2}}d_{1}d_{2}}{2k_{\parallel}}N(k_{\parallel};d_{1})N(k_{\parallel};d_{2})\begin{bmatrix}-k_{x}^{2}&-k_{x}k_{y}&ik_{x}k_{\parallel}\\\
-k_{x}k_{y}&-k_{y}^{2}&ik_{y}k_{\parallel}\\\
ik_{x}k_{\parallel}&ik_{y}k_{\parallel}&k^{2}_{\parallel}\end{bmatrix},$ (9)

where $N$ has a form of $N(k;d)=(1-e^{-kd})/(kd)$. We note that in Eq. (7),
the moments $\tilde{\bm{m}}_{1}$ and $\tilde{\bm{m}}_{2}$ are expressed in the
global coordinate system. After switching to local coordinate systems, we
would get a new matrix $\bm{g}(\bm{k}_{\parallel})$ instead of
$\bm{G}(\bm{k}_{\parallel})$ to express $E_{\rm d}$:

$\begin{aligned} E_{{\rm
d}}=\int\tilde{\bm{m}}_{2}^{\dagger}(\bm{k}_{\parallel})\bm{g}(\bm{k}_{\parallel})\tilde{\bm{m}}_{1}(\bm{k}_{\parallel})\frac{d^{2}\bm{k}_{\parallel}}{(2\pi)^{2}}\end{aligned}.$
(10)

Defining the rotation matrix as

$R_{j}=\begin{bmatrix}{\rm sin}\theta_{j}&{\rm cos}\theta_{j}&0\\\ -{\rm
cos}\theta_{j}&{\rm sin}\theta_{j}&0\\\ 0&0&1\end{bmatrix},$ (11)

and only considering magnons confined to transport in the $x$ axis, which
means $k_{\parallel}=|k_{x}|$ and $k_{y}=0$ (i.e., the simplified case
mentioned in the main text, where each layer is a infinitely long strip along
the $x$ axis and the magnetic moment distribution along the $y$ axis is
uniform), we get

$\displaystyle\bm{g}(k_{x})=R^{-1}_{2}\bm{G}(k_{x})R_{1}=\frac{1}{2}\mu_{0}M_{{\rm
s1}}M_{{\rm
s2}}d_{1}d_{2}N(|k_{x}|;d_{1})N(|k_{x}|;d_{2})\begin{bmatrix}-|k_{x}|{\rm
sin}\theta_{1}{\rm sin}\theta_{2}&-|k_{x}|{\rm cos}\theta_{1}{\rm
sin}\theta_{2}&ik_{x}{\rm sin}\theta_{2}\\\ -|k_{x}|{\rm sin}\theta_{1}{\rm
cos}\theta_{2}&-|k_{x}|{\rm cos}\theta_{1}{\rm cos}\theta_{2}&ik_{x}{\rm
cos}\theta_{2}\\\ ik_{x}{\rm sin}\theta_{1}&ik_{x}{\rm
cos}\theta_{1}&|k_{x}|\end{bmatrix}.$ (12)

Right now, the interlayer dipolar energy is written as

$\displaystyle E_{{\rm
d}}=\int\tilde{\bm{m}}_{2}^{\dagger}(k_{x})\bm{g}(k_{x})\tilde{\bm{m}}_{1}(k_{x})\frac{dk_{x}}{2\pi}=-\mu_{0}M_{{\rm
s}1}d_{1}\int\tilde{\bm{h}}_{1}^{{\rm
d}\dagger}(k_{x})\tilde{\bm{m}}_{1}(k_{x})\frac{dk_{x}}{2\pi}=-\mu_{0}M_{{\rm
s}2}d_{2}\int\tilde{\bm{m}}_{2}^{\dagger}(k_{x})\tilde{\bm{h}}_{2}^{{\rm
d}}(k_{x})\frac{dk_{x}}{2\pi},$ (13)

and the interlayer dipolar field $\tilde{\bm{h}}_{j}^{{\rm d}}(k_{x})$ in each
layer is given by

$\begin{gathered}\tilde{\bm{h}}_{1}^{{\rm d}}(k_{x})=-\frac{1}{\mu_{0}M_{{\rm
s}1}d_{1}}\bm{g}^{\dagger}(k_{x})\tilde{\bm{m}}_{2}(k_{x}),\\\
\tilde{\bm{h}}_{2}^{{\rm d}}(k_{x})=-\frac{1}{\mu_{0}M_{{\rm
s}2}d_{2}}\bm{g}(k_{x})\tilde{\bm{m}}_{1}(k_{x}),\end{gathered}$ (14)

which can be written in the form of

$\begin{gathered}\begin{bmatrix}\tilde{h}^{{\rm d}}_{1x}(k_{x})\\\
\tilde{h}^{{\rm d}}_{1z}(k_{x})\end{bmatrix}=f_{1}\begin{bmatrix}|k_{x}|{\rm
sin}\theta_{1}{\rm sin}\theta_{2}&ik_{x}{\rm sin}\theta_{1}\\\ ik_{x}{\rm
sin}\theta_{2}&-|k_{x}|\end{bmatrix}\begin{bmatrix}\tilde{m}_{2x}(k_{x})\\\
\tilde{m}_{2z}(k_{x})\end{bmatrix},\\\ \begin{bmatrix}\tilde{h}^{{\rm
d}}_{2x}(k_{x})\\\ \tilde{h}^{{\rm
d}}_{2z}(k_{x})\end{bmatrix}=f_{2}\begin{bmatrix}|k_{x}|{\rm
sin}\theta_{1}{\rm sin}\theta_{2}&-ik_{x}{\rm sin}\theta_{2}\\\ -ik_{x}{\rm
sin}\theta_{1}&-|k_{x}|\end{bmatrix}\begin{bmatrix}\tilde{m}_{1x}(k_{x})\\\
\tilde{m}_{1z}(k_{x})\end{bmatrix},\end{gathered}$ (15)

where $f_{1(2)}=M_{{\rm s}2(1)}d_{2(1)}N(|k_{x}|;d_{1})N(|k_{x}|;d_{2})/2$.

There are two additional points we would like to mention:

1) If the moment distribution along the $z$ direction is a higher-order
standing mode within each individual layer, which means
$\bm{m}_{1}(\bm{r})=\bm{m}_{1}(\bm{r}_{\parallel})S_{m}(z;d_{1})$,
$\bm{m}_{2}(\bm{r})=\bm{m}_{2}(\bm{r}_{\parallel})S_{n}(z;d_{2})$ with
$m,n=1,2,...$, and

$\begin{gathered}S_{m}(z;d_{1})=\left\\{\begin{aligned} {\rm cos}(m\pi
z/d_{1})\quad&(-d_{1}<z<0)\\\ 0\quad\quad\;\;\quad&({\rm
otherwise}),\end{aligned}\right.\\\ S_{n}(z;d_{2})=\left\\{\begin{aligned}
{\rm cos}(n\pi z/d_{2})\quad\;&(0<z<d_{2})\\\ 0\quad\quad\;\;\quad\;&({\rm
otherwise}),\end{aligned}\right.\end{gathered}$ (16)

we can then use the similar analysis as we do for the uniform mode. The only
difference is that the factor $N(|k_{x}|;d_{1})$ $N(|k_{x}|;d_{2})$ included
in $f_{j}$ is substituted by $N^{m}(|k_{x}|;d_{1})$ $N^{n}(|k_{x}|;d_{2})$,
where $N^{m}(k;d)=kd[1-(-1)^{m}e^{-kd}]/[(kd)^{2}+(m\pi)^{2}]$. Therefore, a
higher-order standing mode will only decrease the interlayer coupling
strength, but will not qualitatively affect the conclusions on the tunable
Chern bands and chiral surface spin currents.

2) If two layers in Fig. S1(c) are separated by a gap of $\delta d$, a factor
of ${\rm exp}(-k_{\parallel}\delta d)$ will be additionally included in
$f_{j}$. Due to this exponential decay character, we only include the dipolar
interaction between nearest neighboring layers.

##  Appendix B  2D Magnonic Hamiltonian with only the interlayer dipolar
coupling

We continue focusing on the magnetic bilayers. Since each magnetic moment
precesses around its equilibrium orientation, in the small-angle precession
limit, we can express the moment as
$\bm{m}_{j}(\bm{r},t)=m_{jx}(\bm{r})e^{i\omega
t}\hat{x}_{j}+1\cdot\hat{y}_{j}+m_{jz}(\bm{r})e^{i\omega t}\hat{z}$, with the
angular frequency $\omega$. The dynamics of precession is characterized by the
Landau-Lifshitz (LL) equation

$\begin{aligned} \frac{\partial\bm{m}_{j}}{\partial
t}=-\gamma\mu_{0}\bm{m}_{j}\times\bm{h}_{j}\end{aligned},$ (17)

where $\gamma$ is electron’s gyromagnetic ratio, $\mu_{0}$ is vacuum
permeability, $\bm{h}_{j}(\bm{r},t)=h_{jx}(\bm{r})e^{i\omega
t}\hat{x}_{j}+H_{j}\hat{y}_{j}+h_{jz}(\bm{r})e^{i\omega t}\hat{z}$ is the
effective field experienced by each layer, with the static component $H_{j}$
including the external field $H$ and the static interfacial exchange field
$H^{\rm ex}_{j}$. The interfacial exchange energy can be written as [2]

$\displaystyle E^{\rm ex}_{\rm int}$
$\displaystyle=-J\int\bm{m}_{1}(\bm{r}_{\parallel})\cdot\bm{m}_{2}(\bm{r}_{\parallel})d^{2}\bm{r}_{\parallel},$
(18)

where $J$ is the interfacial exchange constant, and the integral is on the
whole interface. The static field contributed by the interfacial exchange is
therefore derived as $H^{\rm ex}_{j}=J{\rm
cos}(\theta_{1}-\theta_{2})/(\mu_{0}M_{{\rm s}j}d_{j})$. In the antiparallel
aligned multilayers, i.e., $J<0$ and $\theta_{2}=\theta_{1}\pm\pi$, the static
interfacial exchange field in each layer is $H^{\rm
ex}_{j}=2|J|/(\mu_{0}M_{{\rm s}j}d_{j})$, when including both the interfaces.

Linearizing Eq. (17) in the small-angle precession limit and transforming it
into the reciprocal space, we have

$\begin{gathered}\frac{i\omega}{\gamma\mu_{0}}\tilde{m}_{jx}=-\tilde{h}_{jz}+H_{j}\tilde{m}_{jz},\\\
\frac{i\omega}{\gamma\mu_{0}}\tilde{m}_{jz}=\tilde{h}_{jx}-H_{j}\tilde{m}_{jx}.\end{gathered}$
(19)

In the absence of interlayer interactions, each individual layer has its
unperturbed eigenfrequencies and eigenstates. The intralayer exchange field is
$\bm{h}^{{\rm ex0}}_{j}=[2A^{{\rm ex}}_{j}/(\mu_{0}M_{{\rm
s}j})]\nabla^{2}\bm{m}_{j}$ [3], with the exchange stiffness constant $A^{{\rm
ex}}_{j}$. In Eq. (19), this corresponds to:

$\tilde{h}^{{\rm ex0}}_{jx(z)}=-\frac{2A^{{\rm
ex}}_{j}k_{x}^{2}}{\mu_{0}M_{{\rm s}j}}\tilde{m}_{jx(z)},$ (20)

The unperturbed eigenfrequencies and eigenstates of each layer are therefore
given by

$\begin{gathered}\omega_{j\pm}=\pm(A_{j}k_{x}^{2}+\Omega_{j}),\\\
\tilde{m}_{j\pm}=\frac{1}{\sqrt{2}}(\tilde{m}_{jx}\pm
i\tilde{m}_{jz}),\end{gathered}$ (21)

where $A_{j}=2A^{{\rm ex}}_{j}\gamma/M_{{\rm s}j}$ characterizes the
intralayer exchange strength and $\Omega_{j}=\gamma\mu_{0}H_{j}$ is the Larmor
precession frequency. From Eq. (21), we see that the unperturbed eigenstate
$\tilde{m}_{j\pm}$ is a right- (left-) handed circularly polarized magnon
mode, which corresponds to the positive (negative) eigenfrequency,
representing a particle (hole).

Applying the intralayer exchange field [Eq. (20)] and interlayer dipolar field
[Eq. (15)] into the linearized LL equation [Eq. (19)], we can get the
eigenvalue equation

$\displaystyle\omega\begin{bmatrix}\tilde{m}_{1+}\\\ \tilde{m}_{2+}\\\
\tilde{m}_{1-}\\\
\tilde{m}_{2-}\end{bmatrix}=\begin{bmatrix}\omega_{1+}&\Delta_{1+}&0&\delta_{1+}\\\
\Delta_{2+}&\omega_{2+}&\delta_{2+}&0\\\
0&\delta_{1-}&\omega_{1-}&\Delta_{1-}\\\
\delta_{2-}&0&\Delta_{2-}&\omega_{1-}\end{bmatrix}\begin{bmatrix}\tilde{m}_{1+}\\\
\tilde{m}_{2+}\\\ \tilde{m}_{1-}\\\ \tilde{m}_{2-}\end{bmatrix},$ (22)

where the diagonal terms are exactly unperturbed eigenfrequencies [Eq. (21)]
and the off-diagonal terms are originated from the interlayer dipolar
coupling, which can be written as

$\begin{gathered}\Delta_{j\pm}=\frac{1}{2}\gamma\mu_{0}f_{j}[k_{x}({\rm
sin}\theta_{1}-{\rm sin}\theta_{2})\pm|k_{x}|(1-{\rm sin}\theta_{1}{\rm
sin}\theta_{2})],\\\ \delta_{j\pm}=\frac{1}{2}\gamma\mu_{0}f_{j}[k_{x}({\rm
sin}\theta_{1}+{\rm sin}\theta_{2})\mp|k_{x}|(1+{\rm sin}\theta_{1}{\rm
sin}\theta_{2})].\end{gathered}$ (23)

Since we are interested in bilayers with antiferromagnetic spin textures with
$\theta_{1}=\pi/2$, $\theta_{2}=-\pi/2$, or $\theta_{1}=-\pi/2$,
$\theta_{2}=\pi/2$, $\delta_{j\pm}$ vanishes and Eq. (22) becomes block
diagonalized.

We can further express the magnonic excitations using the Holstein-Primakoff
(HP) transformation [4]

$\displaystyle\tilde{m}_{j+}\approx\sqrt{\frac{2\gamma\hbar}{M_{{\rm
s}j}S_{xy}d_{j}}}b_{j,k_{x}},\quad\tilde{m}_{j-}\approx\sqrt{\frac{2\gamma\hbar}{M_{{\rm
s}j}S_{xy}d_{j}}}b_{j,k_{x}}^{\dagger}$ (24)

where $S_{xy}$ is the area of the interface, and $b_{j,k_{x}}$
($b_{j,k_{x}}^{\dagger}$) annihilates (creates) a circularly polarized magnon
with wave vector $k_{x}$ in the $j$th layer. If we define
$\bm{b}_{k_{x}}=[b_{1,k_{x}},b_{2,k_{x}}]^{T}$, Eq. (22) can be equivalently
written as

$\displaystyle i\frac{d}{dt}\begin{bmatrix}\bm{b}_{k_{x}}\\\
\bm{b}_{-k_{x}}^{\dagger}\end{bmatrix}=\left[\begin{bmatrix}\bm{b}_{k_{x}}\\\
\bm{b}_{-k_{x}}^{\dagger}\end{bmatrix},\mathcal{\hat{H}}_{\rm
2L}\right]=\bm{\sigma}_{3}\bm{H}_{\rm
2L}(k_{x})\begin{bmatrix}\bm{b}_{k_{x}}\\\
\bm{b}_{-k_{x}}^{\dagger}\end{bmatrix},$ (25)

with

$\begin{gathered}\mathcal{\hat{H}}_{\rm
2L}=\frac{1}{2}\sum_{k_{x}}\begin{bmatrix}\bm{b}_{k_{x}}^{\dagger}&\bm{b}_{-k_{x}}\end{bmatrix}\bm{H}_{\rm
2L}(k_{x})\begin{bmatrix}\bm{b}_{k_{x}}\\\
\bm{b}_{-k_{x}}^{\dagger}\end{bmatrix},\\\
\bm{\sigma}_{3}\equiv\begin{bmatrix}\bm{1}&\bm{0}\\\
\bm{0}&-\bm{1}\end{bmatrix},\quad\bm{H}_{\rm
2L}(k_{x})=\begin{bmatrix}\bm{h}_{\rm 2L}(k_{x})&\bm{0}\\\ \bm{0}&\bm{h}_{\rm
2L}^{*}(-k_{x})\end{bmatrix},\end{gathered}$ (26)

where

$\begin{gathered}\bm{h}_{\rm 2L}(k_{x})=\begin{bmatrix}\omega_{1}&\Delta\\\
\Delta&\omega_{2}\end{bmatrix},\\\
\omega_{j}=\omega_{j+}=A_{j}k_{x}^{2}+\Omega_{j},\\\
\Delta=\sqrt{\Delta_{1+}\Delta_{2+}}=\frac{1}{2}\gamma\mu_{0}\sqrt{f_{1}f_{2}}[k_{x}({\rm
sin}\theta_{1}-{\rm sin}\theta_{2})+|k_{x}|(1-{\rm sin}\theta_{1}{\rm
sin}\theta_{2})].\end{gathered}$ (27)

The bosonic Bogoliubov-de Gennes (BdG) Hamiltonian $\bm{H}_{\rm 2L}(k_{x})$
only possesses the diagonal blocks, so the bilayer Hamiltonian operator can be
simplified as

$\mathcal{\hat{H}}_{\rm 2L}=\sum_{k_{x}}\bm{b}_{k_{x}}^{\dagger}\bm{h}_{\rm
2L}(k_{x})\bm{b}_{k_{x}}.$ (28)

As for the antiparallelly aligned multilayers shown in Fig. 1 of the main
text, we can regard it as the stack of bilayers along the $z$ direction.
Particularly, when treating the dipolar coupling between two layers within a
unit cell, we can set $\theta_{1}=\pi/2$, $\theta_{2}=-\pi/2$, while we set
$\theta_{1}=-\pi/2$, $\theta_{2}=\pi/2$ for two layers belonging to
neighboring cells. Therefore, we can get the intracell (intercell) coupling
term $\Delta_{\rm S(D)}$ by assigning the corresponding values of $\theta_{1}$
and $\theta_{2}$ into the expression of $\Delta$ in Eq. (27):

$\begin{gathered}\Delta_{\rm
S}=\Delta|_{\theta_{1}=\pi/2,\theta_{2}=-\pi/2}=B\cdot(|k_{x}|+k_{x}),\\\
\Delta_{\rm
D}=\Delta|_{\theta_{1}=-\pi/2,\theta_{2}=\pi/2}=B\cdot(|k_{x}|-k_{x}),\end{gathered}$
(29)

with $B=\mu_{0}\gamma\sqrt{f_{1}f_{2}}$ characterizing the coupling strength.
The multilayer magnonic Hamiltonian is therefore written as

$\displaystyle\hat{\mathcal{H}}=\hat{\mathcal{H}}_{0}+\hat{\mathcal{H}}_{\rm
dip}=\sum_{k_{x},j,n}\omega_{j}b_{jn,k_{x}}^{\dagger}b_{jn,k_{x}}+\sum_{k_{x},j,n}\left(\Delta_{\rm
S}b_{1n,k_{x}}^{\dagger}b_{2n,k_{x}}+\Delta_{\rm
D}b_{1n,k_{x}}^{\dagger}b_{2,n-1,k_{x}}+\mathrm{H.c.}\right).$ (30)

We can get the bulk Hamiltonian [Eq. (1) of the main text], with the periodic
boundary condition along the $z$ direction, and we can solve for the explicit
surface states [Figs. 2(a) to 2(c) of the main text], with the open boundary
condition along the $z$ direction.

##  Appendix C  3D Magnonic Hamiltonian with all interactions included

In the generic case, magnons can propagate within the whole $xy$ plane in each
individual layer, i.e., with the in-plane momenta
$\bm{k_{\parallel}}=k_{x}\hat{x}+k_{y}\hat{y}$. Besides, both the intralayer
dipolar interaction and the dynamic interfacial exchange interaction should be
included.

First, let’s derive the interlayer dipolar field with finite $k_{y}$. For two
layers within a unit cell, we have $\theta_{1}=\pi/2$ and $\theta_{2}=-\pi/2$.
Substituting them into Eq. (11) gives

$\displaystyle R_{1}^{\prime}=\begin{bmatrix}1&0&0\\\ 0&1&0\\\
0&0&1\end{bmatrix},\quad R_{2}^{\prime}=\begin{bmatrix}-1&0&0\\\ 0&-1&0\\\
0&0&1\end{bmatrix}.$ (31)

Applying these rotation matrices into Eq. (9), we get

$\displaystyle\bm{g}^{\prime}(\bm{k}_{\parallel})=$
$\displaystyle(R_{2}^{\prime})^{-1}\bm{G}(\bm{k}_{\parallel})R_{1}^{\prime}=\frac{\mu_{0}M_{{\rm
s1}}M_{{\rm
s2}}d_{1}d_{2}}{2k_{\parallel}}N(k_{\parallel};d_{1})N(k_{\parallel};d_{2})\begin{bmatrix}k_{x}^{2}&k_{x}k_{y}&-ik_{x}k_{\parallel}\\\
k_{x}k_{y}&k_{y}^{2}&-ik_{y}k_{\parallel}\\\
ik_{x}k_{\parallel}&ik_{y}k_{\parallel}&k^{2}_{\parallel}\end{bmatrix}.$ (32)

Right now, the interlayer dipolar energy is written as

$\displaystyle E_{{\rm d}}^{\prime}=$
$\displaystyle\int\tilde{\bm{m}}_{2}^{\dagger}(\bm{k}_{\parallel})\bm{g}^{\prime}(\bm{k}_{\parallel})\tilde{\bm{m}}_{1}(\bm{k}_{\parallel})\frac{d^{2}\bm{k}_{\parallel}}{(2\pi)^{2}}=-\mu_{0}M_{{\rm
s}1}d_{1}\int\tilde{\bm{h}}_{1}^{{\rm
d}\dagger}(\bm{k}_{\parallel})\tilde{\bm{m}}_{1}(\bm{k}_{\parallel})\frac{d^{2}\bm{k}_{\parallel}}{(2\pi)^{2}}=-\mu_{0}M_{{\rm
s}2}d_{2}\int\tilde{\bm{m}}_{2}^{\dagger}(\bm{k}_{\parallel})\tilde{\bm{h}}_{2}^{{\rm
d}}(\bm{k}_{\parallel})\frac{d^{2}\bm{k}_{\parallel}}{(2\pi)^{2}},$ (33)

and the interlayer dipolar field in each layer is give by

$\begin{gathered}\tilde{\bm{h}}_{1}^{{\rm
d}}(\bm{k}_{\parallel})=-\frac{1}{\mu_{0}M_{{\rm
s}1}d_{1}}\bm{g}^{\prime\dagger}(\bm{k}_{\parallel})\tilde{\bm{m}}_{2}(\bm{k}_{\parallel}),\\\
\tilde{\bm{h}}_{2}^{{\rm d}}(\bm{k}_{\parallel})=-\frac{1}{\mu_{0}M_{{\rm
s}2}d_{2}}\bm{g}^{\prime}(\bm{k}_{\parallel})\tilde{\bm{m}}_{1}(\bm{k}_{\parallel}),\end{gathered}$
(34)

which can be written in the form of components

$\displaystyle\begin{bmatrix}\tilde{h}^{{\rm d}}_{jx}(\bm{k}_{\parallel})\\\
\tilde{h}^{{\rm d}}_{jz}(\bm{k}_{\parallel})\end{bmatrix}$
$\displaystyle=f_{j}^{\prime}\begin{bmatrix}-k_{x}^{2}/k_{\parallel}&ik_{x}\\\
-ik_{x}&-k_{\parallel}\end{bmatrix}\begin{bmatrix}\tilde{m}_{vx}(\bm{k}_{\parallel})\\\
\tilde{m}_{vz}(\bm{k}_{\parallel})\end{bmatrix}.$ (35)

where $v=2$ for $j=1$ and $v=1$ for $j=2$. $f_{j}^{\prime}$ can be obtained by
replacing $|k_{x}|$ with $|\bm{k}_{\parallel}|$ in the expression of $f_{j}$.

Next, let’s derive the intralayer dipolar field, which is similar to what
we’ve done for the interlayer dipolar field. The difference is that there is
only one continuous body now, so $\bm{r}$ and $\bm{r^{\prime}}$ are in the
same layer, as shown in Fig. S1(d). The intralayer dipolar energy of the $j$th
layer is

$\displaystyle E_{{\rm
d}0,j}=\frac{1}{2}\int\tilde{\bm{m}}_{j}^{\dagger}(\bm{k}_{\parallel})\bm{G}_{j}(\bm{k}_{\parallel})\tilde{\bm{m}}_{j}(\bm{k}_{\parallel})\frac{d^{2}\bm{k}_{\parallel}}{(2\pi)^{2}}.$
(36)

$\bm{G}_{j}({\bm{k}_{\parallel}})$ is calculated to be

$\displaystyle\bm{G}_{j}({\bm{k}_{\parallel}})=\frac{1}{3}\mu_{0}M_{{\rm
s}j}^{2}d_{j}\bm{1}-\mu_{0}M_{{\rm
s}j}^{2}d_{j}\frac{1}{k_{\parallel}^{2}}\begin{bmatrix}k_{x}^{2}(1-N_{j}^{\prime})&k_{x}k_{y}(1-N_{j}^{\prime})&0\\\
k_{x}k_{y}(1-N_{j}^{\prime})&k_{y}^{2}(1-N_{j}^{\prime})&0\\\
0&0&k^{2}_{\parallel}N_{j}^{\prime}\end{bmatrix},$ (37)

where $N_{j}^{\prime}=N(|k_{\parallel}|;d_{j})$. After switching to the local
coordinate system, we find
$\bm{g}_{j}(\bm{k}_{\parallel})=(R_{j}^{\prime})^{-1}\bm{G}_{j}(\bm{k}_{\parallel})R_{j}^{\prime}=\bm{G}_{j}(\bm{k}_{\parallel})$,
where $R_{j}^{\prime}$ is shown in Eq. (31). Right now, the intralayer dipolar
energy is written as

$\displaystyle E_{{\rm d}0,j}=$
$\displaystyle\frac{1}{2}\int\tilde{\bm{m}}_{j}^{\dagger}(\bm{k}_{\parallel})\bm{g}_{j}(\bm{k}_{\parallel})\tilde{\bm{m}}_{j}(\bm{k}_{\parallel})\frac{d^{2}\bm{k}_{\parallel}}{(2\pi)^{2}}=-\frac{1}{2}\mu_{0}M_{{\rm
s}j}d_{j}\int\tilde{\bm{m}}_{j}^{\dagger}(\bm{k}_{\parallel})\tilde{\bm{h}}_{j}^{{\rm
d0}}(\bm{k}_{\parallel})\frac{d^{2}\bm{k}_{\parallel}}{(2\pi)^{2}},$ (38)

and the intralayer dipolar field $\tilde{\bm{h}}^{{\rm
d0}}_{j}(\bm{k}_{\parallel})$ in each layer is given by

$\displaystyle\tilde{\bm{h}}_{j}^{{\rm
d0}}(\bm{k}_{\parallel})=-\frac{1}{\mu_{0}M_{{\rm
s}j}d_{j}}\bm{g}_{j}(\bm{k}_{\parallel})\tilde{\bm{m}}_{j}(\bm{k}_{\parallel}),$
(39)

which can be written in the form of

$\displaystyle\begin{bmatrix}\tilde{h}^{{\rm d0}}_{jx}(\bm{k}_{\parallel})\\\
\tilde{h}^{{\rm d0}}_{jz}(\bm{k}_{\parallel})\end{bmatrix}$
$\displaystyle=M_{{\rm
s}j}\begin{bmatrix}(1-N_{j}^{\prime})k_{x}^{2}/k_{\parallel}^{2}+\frac{1}{3}&0\\\
0&N_{j}^{\prime}+\frac{1}{3}\end{bmatrix}\begin{bmatrix}\tilde{m}_{jx}(\bm{k}_{\parallel})\\\
\tilde{m}_{jz}(\bm{k}_{\parallel})\end{bmatrix}.$ (40)

On the other hand, from Eq. (18), we can get the dynamic interfacial exchange
field $\tilde{\bm{h}}_{j}^{{\rm ex}}$ in each layer:

$\displaystyle\begin{bmatrix}\tilde{h}^{{\rm ex}}_{jx}(\bm{k}_{\parallel})\\\
\tilde{h}^{{\rm ex}}_{jz}(\bm{k}_{\parallel})\end{bmatrix}=H^{\rm
ex}_{j}\begin{bmatrix}1/2&0\\\
0&-1/2\end{bmatrix}\begin{bmatrix}\tilde{m}_{vx}(\bm{k}_{\parallel})\\\
\tilde{m}_{vz}(\bm{k}_{\parallel})\end{bmatrix}.$ (41)

Applying the intralayer exchange field [replacing $k_{x}$ by
$\bm{k}_{\parallel}$ in Eq. (20)], interlayer dipolar field [Eq. (35)],
intralayer dipolar field [Eq. (40)], and dynamic interfacial exchange field
[Eq. (41)], into the linearized LL equation [Eq. (19)], we can get the
eigenvalue equation

$\displaystyle\omega\begin{bmatrix}\tilde{m}_{1+}\\\ \tilde{m}_{2+}\\\
\tilde{m}_{1-}\\\
\tilde{m}_{2-}\end{bmatrix}=\begin{bmatrix}\omega_{1}^{\prime}&\Delta_{1+}^{\prime}&\delta_{1}&\delta_{1}^{\prime}\\\
\Delta_{2+}^{\prime}&\omega_{2}^{\prime}&\delta_{2}^{\prime}&\delta_{2}\\\
-\delta_{1}&-\delta_{1}^{\prime}&-\omega_{1}^{\prime}&\Delta_{1-}^{\prime}\\\
-\delta_{2}^{\prime}&-\delta_{2}&\Delta_{2-}^{\prime}&-\omega_{2}^{\prime}\end{bmatrix}\begin{bmatrix}\tilde{m}_{1+}\\\
\tilde{m}_{2+}\\\ \tilde{m}_{1-}\\\ \tilde{m}_{2-}\end{bmatrix},$ (42)

with

$\begin{gathered}\omega_{j}^{\prime}=A_{j}\bm{k}_{\parallel}^{2}+\Omega_{j}+\frac{\gamma\mu_{0}M_{{\rm
s}j}}{2}\left[(1-N_{j}^{\prime})\frac{k_{x}^{2}}{k_{\parallel}^{2}}+N_{j}^{\prime}+\frac{2}{3}\right],\\\
\Delta_{j\pm}^{\prime}=\frac{1}{2}\gamma\mu_{0}f_{j}^{\prime}\left[\pm\left(k_{\parallel}+\frac{k_{x}^{2}}{k_{\parallel}}\right)+2k_{x}\right],\\\
\delta_{j}=\frac{\gamma\mu_{0}M_{{\rm
s}j}}{2}\left[(1-N_{j}^{\prime})\frac{k_{x}^{2}}{k_{\parallel}^{2}}-N_{j}^{\prime}\right],\\\
\delta_{j}^{\prime}=\frac{\gamma\mu_{0}}{2}\left[f_{j}^{\prime}\cdot\left(k_{\parallel}-\frac{k_{x}^{2}}{k_{\parallel}}\right)+H^{\rm
ex}_{j}\right].\end{gathered}$ (43)

Eqs. (42) and (43) give the eigenvalue equation for two layers within the same
unit cell ($\theta_{1}=\pi/2$ and $\theta_{2}=-\pi/2$). As for two layers
belonging to neighboring cells ($\theta_{1}=-\pi/2$ and $\theta_{2}=\pi/2$),
the only difference in the eigenvalue equation is

$\Delta_{j\pm}^{\prime}=\frac{1}{2}\gamma\mu_{0}f_{j}^{\prime}\left[\pm\left(k_{\parallel}+\frac{k_{x}^{2}}{k_{\parallel}}\right)-2k_{x}\right].$
(44)

Then we can do the HP transformation on $\tilde{m}_{j+}$ and $\tilde{m}_{j-}$
[Eq. (24)] and apply the periodic boundary condition along the $z$ direction.
We finally get the bulk magnonic Hamiltonian for the 3D multilayers [Eqs. (3)
to (6) of the main text].

## References

  * Cohen and Keffer [1955] M. H. Cohen and F. Keffer, Phys. Rev. 99, 1128 (1955).
  * Berkowitz and Takano [1999] A. Berkowitz and K. Takano, J. Magn. Magn. Mater. 200, 552 (1999).
  * Klingler  _et al._ [2014] S. Klingler, A. V. Chumak, T. Mewes, B. Khodadadi, C. Mewes, C. Dubs, O. Surzhenko, B. Hillebrands, and A. Conca, J. Phys. D: Appl. Phys. 48, 015001 (2014).
  * Holstein and Primakoff [1940] T. Holstein and H. Primakoff, Phys. Rev. 58, 1098 (1940).


# Driven bright solitons on a mid-infrared laser chip

Dmitry Kazakov<EMAIL_ADDRESS>Harvard John A. Paulson School of
Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA
Theodore P. Letsou  Harvard John A. Paulson School of Engineering and Applied
Sciences, Harvard University, Cambridge, MA 02138, USA  Department of
Electrical Engineering and Computer Science, Massachusetts Institute of
Technology, Cambridge, MA 02142, USA  Marco Piccardo  Harvard John A. Paulson
School of Engineering and Applied Sciences, Harvard University, Cambridge, MA
02138, USA  Department of Physics, Instituto Superior Técnico, Universidade de
Lisboa, 1049-001 Lisbon, Portugal  Instituto de Engenharia de Sistemas e
Computadores – Microsistemas e Nanotecnologias (INESC MN), 1000-029 Lisbon,
Portugal  Lorenzo L. Columbo  Dipartimento di Elettronica e Telecomunicazioni,
Politecnico di Torino, 10129 Torino, Italy  Massimo Brambilla  Dipartimento di
Fisica Interateneo, Università e Politecnico di Bari, 70126 Bari, Italy  CNR-
Istituto di Fotonica e Nanotecnologie, 70126 Bari, Italy  Franco Prati
Dipartimento di Scienza e Alta Tecnologia, Università dell’Insubria, 22100
Como, Italy  Sandro Dal Cin  Institute of Solid State Electronics, TU Wien,
1040 Vienna, Austria  Maximilian Beiser  Institute of Solid State Electronics,
TU Wien, 1040 Vienna, Austria  Nikola Opačak  Institute of Solid State
Electronics, TU Wien, 1040 Vienna, Austria  Pawan Ratra  Harvard John A.
Paulson School of Engineering and Applied Sciences, Harvard University,
Cambridge, MA 02138, USA  Department of Electrical and Electronic Engineering,
Imperial College London, London SW7 2BX, United Kingdom  Michael Pushkarsky
DRS Daylight Solutions, San Diego, CA 92127, USA  David Caffey  DRS Daylight
Solutions, San Diego, CA 92127, USA  Timothy Day  DRS Daylight Solutions, San
Diego, CA 92127, USA  Luigi A. Lugiato  Dipartimento di Scienza e Alta
Tecnologia, Università dell’Insubria, 22100 Como, Italy  Benedikt Schwarz
Harvard John A. Paulson School of Engineering and Applied Sciences, Harvard
University, Cambridge, MA 02138, USA  Institute of Solid State Electronics, TU
Wien, 1040 Vienna, Austria  Federico Capasso<EMAIL_ADDRESS>Harvard
John A. Paulson School of Engineering and Applied Sciences, Harvard
University, Cambridge, MA 02138, USA

###### Abstract

Despite the ongoing progress in integrated optical frequency comb technology
[1], compact sources of short bright pulses in the mid-infrared wavelength
range from 3 $\upmu$m to 12 $\upmu$m so far remained beyond reach. The state-
of-the-art ultrafast pulse emitters in the mid-infrared are complex, bulky,
and inefficient systems based on the downconversion of near-infrared or
visible pulsed laser sources. Here we show a purely DC-driven semiconductor
laser chip that generates one picosecond solitons at the center wavelength of
8.3 $\upmu$m at GHz repetition rates. The soliton generation scheme is akin to
that of passive nonlinear Kerr resonators [2]. It relies on a fast bistability
in active nonlinear laser resonators, unlike traditional passive mode-locking
which relies on saturable absorbers [3] or active mode-locking by gain
modulation in semiconductor lasers [4]. Monolithic integration of all
components — drive laser, active ring resonator, coupler, and pump filter —
enables turnkey generation of bright solitons that remain robust for hours of
continuous operation without active stabilization. Such devices can be readily
produced at industrial laser foundries using standard fabrication protocols.
Our work unifies the physics of active and passive microresonator frequency
combs, while simultaneously establishing a technology for nonlinear integrated
photonics in the mid-infrared [5].

Short optical pulses have revolutionized a range of applications, from high-
resolution imaging and ultrafast spectroscopy to optical communications [6],
laser-based medical procedures [7], and light ranging [8]. They play a pivotal
role in nonlinear optics, enabling a variety of phenomena such as
supercontinuum generation [9] and optical frequency conversion that are at the
heart of optical atomic clocks [10]. Over the past two decades, there has been
a remarkable progress in miniaturizing pulsed optical sources and
transitioning them from tabletop experimental setups to compact photonic
integrated chips. The two primary technologies driving this effort are
semiconductor mode-locked lasers (SMLLs) and nonlinear microresonator
frequency combs [1]. These technologies may reduce the size and complexity of
optical atomic clocks [11, 12], pave the way to terabit-per-second
telecommunication links [13, 3], and provide new tools for linear and
nonlinear absorption spectroscopy [14, 15]. The near-infrared wavelength range
($0.8-2.5$ $\upmu$m) has an abundance of compact pulsed sources based both on
SMLLs and nonlinear resonators. In contrast, there has been a scarce number of
demonstrations of integrated photonic pulse generators in the mid-infrared
wavelength range ($3-12$ $\upmu$m), strategic to applications in gas sensing
and spectroscopy. The only devices in this range that are compatible with
photonic integration and capable of directly generating coherent mid-infrared
radiation are interband cascade lasers (ICLs) and quantum cascade lasers
(QCLs).

While recently there have been significant advancements in pulse generation
using these platforms, major limitations still exist that hinder their
deployment in applications. ICLs, that operate at $3-6$ $\upmu$m, have shown
modest milliwatt-level average optical power and the pulse durations limited
to several picoseconds [4]. QCLs, which are most efficient in the $4.5-12$
$\upmu$m range, now routinely reach Watt-level output power [16, 17] and have
recently entered the femtosecond regime [18]. However, achieving subpicosecond
pulse duration required spectral broadening and external pulse shaping
involving strong radiofrequency modulation of the laser bias, a tabletop delay
line, and an optical isolator to prevent laser destabilization due to optical
feedback. These aspects pose a challenge for miniaturization of the system.

We address this challenge by introducing a new feedback-insensitive and
isolator-free semiconductor laser architecture, capable of generating bright
picosecond, mid-infrared pulses directly on a laser chip without any external
compression. Pulse generation is enabled by a new driving scheme that does not
rely on active or passive mode-locking as in SMLLs; instead, it is inspired by
the soliton excitation techniques developed for passive Kerr microresonators
and is enabled by an experimental unification of the physics of dissipative
solitons in passive and active cavities based on a generalized Lugiato-Lefever
equation (GLLE) [19]. This equation extends the LLE — which was originally
formulated for resonators without population inversion [20] — to ring lasers
driven above threshold [21, 22]. Although we perform our demonstration using a
QCL emitting at $8.3$ $\upmu$m, this approach applies to any device of such
class across the mid-infrared. The chip-scale pulse generators shown here thus
extend the spectral gamut covered by the SMLLs and Kerr microresonators to the
entire mid-infrared range (Fig. 1a).

Figure 1: Pulse generation in integrated resonator systems. a. Spectral
coverage of the three photonic integrated technology stacks for on-chip short
pulse generation. Semiconductor mode-locked laser (SMLL) and Kerr
microresonator spectral coverage data are taken from ref. [1]. b. SMLL
comprises a forward biased gain section and a reverse biased saturable
absorber (SA) section. Reverse biasing the SA section
($\mathrm{P}_{\mathrm{SA}}=-0.8\mathrm{P}_{\mathrm{th}}$) and increasing the
pumping of the forward-biased gain section allow a transition from a CW
waveform to a mode-locked pulsed operation. The panels show the simulated
waveforms and optical spectra of the intracavity field in a SMLL at different
values of the gain section pumping. c. Passive Kerr resonator optically pumped
by an external laser experiences a bistable behavior in the output intensity
as a function of the detuning $\uptheta$ of the laser wavelength from the Kerr
cavity resonance. Here the middle panel shows the mean intracavity intensity
as function of detuning. The intensity values are artificially flipped by
multiplying with $-1$ to facilitate the visual comparison with the
experimental plots later in the manuscript, where not the intracavity
intensity, but the measured output intensity is shown. On a forward scan
(blue-detuned to red-detuned, orange curve in the middle panel) the resonator
transmission shows a series of steps signifying generation of one or several
cavity solitons, which cannot be generated on a backward scan (gray curve in
the middle panel). The simulations are based on the GLLE [19]. State-of-the-
art fabrication technology now allows monolithic integration of the pump laser
with a passive resonator [23]. d. Active quantum cascade (QC) resonator with
an on-chip integrated drive laser. The soliton excitation scheme is similar to
that of passive Kerr resonators and entails scanning the frequency of the
drive laser through the lasing resonance of the active racetrack resonator.
The associated output intensity (orange curve — forward scan, gray curve —
backward scan) exhibits optical bistability and a characteristic step in
transmission when generating one or several cavity solitons on a forward
wavelength scan. Figure 2: Driven bright solitons in an active resonator. This
demonstration relies on an off-chip external drive laser. a. Experimental
setup for external soliton driving. SWEEP, function generator used to apply
triangular sawtooth modulation to the integrated heater (see Methods for
further details on resonance tuning). ECL, external cavity laser. POL,
polarizer. RT, racetrack resonator. WG, waveguide coupler. IFG, interferogram
signal. OUT, output intensity detector. DAQ, data acquisition board. The high-
resolution optical spectrum shows RT and ECL in proximity to each other (ECL
on the blue side of the RT) at the onset of the detuning scan. b. Output
intensity (orange) and interferogram (gray) acquired upon a forward wavelength
scan (blue-detuned to red-detuned). c. Short-time FFT of the output intensity
trace on the forward scan showing the heterodyne notes at the onset of
injection locking of the RT lasing frequency to that of the ECL. d. Output
intensity and interferogram on a backward wavelength scan (red-detuned to
blue-detuned). e. Short-time FFT of the output intensity trace on the backward
scan showing the characteristic step signifying soliton generation. f.
Overlaid output intensity curves over one thousand backward scans, showing
multistability of the soliton states. The traces are aligned so that the first
steps in the output intensity coincide. g. Reconstructed output temporal
waveform of the driven soliton generator over two consecutive roundtrips. h.
Optical spectrum of the waveform in g. The inset shows the corresponding RF
spectrum of the intermode beatnote. i. Optical spectra of multiple soliton
states, obtained on the backward scan, corresponding to the different steps in
the output intensity seen in f. Figure 3: On-chip pump filtering. a. Optical
microscope image of the active resonator photonic integrated chip. The chip
contains two identical racetrack resonators connected via a bus waveguide —
one for soliton generation, another one — for pump filtering. Integrated
heaters are implemented as concentric racetracks on the inner side of the
soliton generator and of the filter. b. Spectrogram of the soliton state from
the left racetrack resonator in a, as the filter response is tuned with the
pump current. Critically coupled ring resonator progressively filters out comb
teeth one by one, shown by the dashed line. c. The filter response at the pump
mode wavelength (1227 cm-1), showing the extinction ratio (ER) of 42 dB. d.
Optical spectra of the soliton generator when the filter is switched off and
switched on. The suppression of the pump frequency by 45 dB is achieved by
fine-tuning the current of the filter and of the integrated heater. e.
Reconstructed temporal waveforms of the soliton generator over two consecutive
cavity roundtrips with filter off and filter on. Figure 4: Integrated turnkey
soliton generator. a. Optical microscope image of the QC photonic integrated
chip. Shown chip contains two identical devices each comprised by four
components: a Fabry-Perot drive laser (FP), a waveguide coupler (WG), a
resistive heater (HT) and a racetrack resonator (RT). Insets show blown up
micrographs of the coupling regions between the FP and the WG (butt coupler)
and between the WG and the RT (directional coupler). b. Top: output intensity
(orange), interferogram (gray), and the optical spectra of the RT resonator on
a backward and on a forward scan of the FP through the RT lasing resonance.
Bottom: spectrogram of the radio frequency beat note on the backward and on
the forward scan as function of the detuning. c. Power spectral density (PSD)
and the reconstructed temporal waveform of the soliton state obtained on the
forward wavelength scan. d. Interferogram (gray) and output intensity (orange)
when the soliton generator is turned on and off. The same state is recovered
at each power cycle following the on-off algorithm that involves adjusting the
biases of RT, FP, and HT. e. PSD in the soliton state, acquired at equally
distributed time intervals, showing an uninterrupted and intact soliton state
over a span of four hours.

Pulsed operation of a semiconductor laser requires both a mechanism for pulse
formation and a mechanism for pulse stabilization. In SMMLs the pulse is
formed and stabilized by an intracavity saturable absorber, that opens a short
temporal window of net gain only for the optical pulse, and attenuates the
background, preventing its growth and pulse destabilization (Fig. 1b). A
competing destabilizing effect of the optical gain, that acts such as to
amplify the noise spikes on the low-intensity background, is weak if the
active medium is slow as the gain does not have the time to recover after
being depleted by the pulse. In contrast, in lasers with a fast gain recovery
a saturable absorber is not able to stabilize the pulse. It is for this
reason, that mid-infrared QCLs, where the gain recovery time is typically less
than a picosecond, up to now could not be passively mode-locked to emit
pulses. In free-running standing-wave cavity Fabry-Perot QCLs, the locked
frequency comb state corresponds to a frequency-modulated (FM) wave with a
nearly flat intensity in the time domain [24, 25, 26, 27]. In traveling-wave
ring QCLs frequency comb operation yields gray or dark — instead of bright —
pulses, on a high intensity background [28, 29].

An alternative pulse stabilizing mechanism is optical bistability that may
occur in resonators where the optical susceptibility of the host material
depends on the intracavity field intensity. In passive resonators, the
susceptibility can be taken as purely real in the limit of large detuning of
the frequency of optical field from the material absorption frequency. The
intensity-dependent real part of susceptibility, the refractive index, is at
the origin of stable soliton formation in passive systems with no optical
gain, such as fiber loop resonators [30] and waveguide-based microresonators
[2]. In these systems, a coherent CW laser field injected into the resonator
may break up into one or multiple pulses that are formed by the competing
action of the nonlinearity and the cavity dispersion. The typical soliton
excitation scheme involves scanning the wavelength of the external pump laser
through resonance from the blue-detuned side to the red-detuned side (forward
scan) [31]. During the scan, the CW pump undergoes an instability and forms
one or several cavity solitons (Fig. 1c). Once excited, these solitons remain
stable against noise when the driving field is effectively red-detuned from
the cavity resonance. Solitons cannot be excited by scanning the laser
wavelength from the red-detuned to the blue-detuned side of the resonance
(backward scan) — a manifestation of the bistable behavior of this optical
system (Fig. 1c).

Here we experimentally demonstrate hybrid nonlinear photonic devices by
applying the framework of soliton generation in passive Kerr microresonators
to a semiconductor laser above its threshold. We consider a traveling-wave
ring cavity filled with the optical gain medium where the bistability stems
from the resonant third-order nonlinearity. In a laser cavity the medium
susceptibility is as well intensity-dependent, by virtue of gain saturation.
Temporal changes in the field intensity lead to changes in the imaginary part
of the susceptibility. The real part of the susceptibility, inherently
dependent on its imaginary part, will equally follow changes in the
intracavity field intensity. This behaviour leads to the refractive
nonlinearity that is akin to the Kerr nonlinearity of passive microresonators.
In a semiconductor with fast gain, the refractive nonlinearity can be
quantified by the laser linewidth enhancement factor (LEF) [32], which in QCLs
has been shown to originate from the Bloch gain [33]. Unlike the bulk crystal
Kerr nonlinearity, the nonlinearity induced by the LEF is strongly dispersive,
and the nonlinear coefficient may even attain both positive and negative
values at different frequencies within the gain bandwidth [34].

The nonlinear active resonators have two important differences with respect to
passive Kerr resonators. First, whereas the gain medium is fast, with the
typical gain recovery time on the order of a few hundred femtoseconds, the
nonlinearity itself that arises from the saturation of the fast gain is still
slower than the almost instantaneous Kerr nonlinearity of passive resonators.
The slow nonlinearity limits the bandwidth over which it acts, albeit stronger
than a fast nonlinearity [35]. Second, unlike in an optically pumped
microresonator, the optical pump field is generated directly inside the ring
laser cavity above its threshold. In a free-running laser, the pump frequency
is tightly linked to the position of the laser cold cavity resonance, and its
detuning from the resonance cannot be freely controlled; the laser operates at
an optimal wavelength located between the peak of the cavity resonance and the
peak of the saturated gain [36]. As a result, the bistability cannot be
reached without a control over the field detuning, and thus, bright solitons
cannot be generated in free-running lasers [29].

Injecting an external control optical field from a wavelength-tunable laser
enables the bistability and allows one to overcome this limitation. In such an
externally driven scheme, the intracavity field will injection lock to the
drive field when the frequency detuning between the two is sufficiently small
(Fig. 1d). The intracavity field frequency will follow the frequency of the
drive field — as long as the drive field remains within the finite locking
range of frequencies — thus effectively decoupling it from the cavity
resonance. Previously we theoretically proposed that such a scheme would allow
the generation of optical solitons in a laser once the intracavity field is
injection locked to the drive field and the intensity of the driving field is
first increased, and then decreased [19]. Alternatively, akin to the soliton
generation scheme of Kerr microresonators, we find that solitons can be
generated by performing wavelength scans of the drive field through the lasing
resonance (Fig. 1d). The intracavity field injection locks to the drive and
generates one or several cavity solitons — in direct analogy with passive Kerr
microresonators.

We generate solitons in active racetrack (RT) QC resonators, operated above
their lasing threshold [37]. The RT operates as a single-mode laser when a
drive signal injected through an on-chip directional coupler from an external
cavity laser (ECL) is detuned far off the RT lasing resonance (Fig. 2a). The
RT QCL operates in a unidirectional regime — an essential prerequisite to the
soliton generation in active laser resonators [19]. Furthermore,
unidirectional operation makes the ring laser essentially feedback-insensitive
— lifting the requirement of an optical isolator. The back-reflected wave
quickly dies out inside of the resonator since the gain is strongly saturated
by the forward propagating intracavity field [29]. By sweeping the current of
the integrated heater (HT) we can move the RT lasing resonance with respect to
the drive to perform the wavelength scans [38]. In the backward wavelength
scan (RT lasing resonance scanning from the red-detuned side to the blue-
detuned side of the ECL frequency), the RT remains in single-mode operation,
when being both in resonance and out of resonance (Fig. 2b, c). On the forward
wavelength scan (RT tuned from blue-detuned to red-detuned), the output
intensity features a characteristic step-like drop, at which a soliton is
generated (Fig. 2d, g). This step is a defining feature of soliton generation
in both passive Kerr resonators and fiber resonators [31, 39]. Importantly,
the RT remains injection-locked to the ECL drive field throughout the backward
scan, (Fig. 2e). The soliton spectrum is comprised of a strong pump and drive
field line and a family of equidistant modes on one side of the pump line
(Fig. 2h). Performing identical backward scans multiple times reveals the
existence of several possible steps in the output intensity (Fig. 2f)
signifying the multistability of the soliton states — a feature inherent to
solitons. Single- and multisoliton states can be excited in a stochastic
manner by performing the same scan multiple times (Fig. 2h, i).

We found the salient spectral asymmetry to be a general feature of our active
resonator solitons. It stands in a strong contrast with the highly symmetric
spectra of passive resonator solitons. While the spectral asymmetry is not
captured by the GLLE model, we hint towards three factors that may be causing
it. First, taking into account higher order dispersion in the GLLE yields
asymmetric states in simulations. Second, the highly dispersive nature of the
nonlinearity induced by the LEF may as well a contributing factor to the
observed spectral shape. Third, the GLLE assumes laser operation with no
detuning of the pump field from the gain peak. In reality, it is hardly ever
the case, and the intracavity field is likely detuned from the gain peak by a
finite value. We have found the detuning of the drive laser frequency from the
gain peak to be a contributing factor to the spectral asymmetry.

One impeding characteristic of the driven cavity solitons is the presence of a
strong background driving field that manifests itself in a mode $20-30$ dB
greater than the rest of the soliton spectrum. This feature necessitates the
use of detectors with high dynamic range and may cause gain depletion of
optical amplifiers by the unnecessarily strong laser line. In passive Kerr
resonators the pump is typically filtered out by an external fiber Bragg
grating [40] or by outcoupling the intracavity field through an add-drop
waveguide [41]. The first approach is not compatible with photonic integration
as the mid-infrared optics is based on free-space components. The second
approach would not work with the active resonators presented here, since the
pump field that drives soliton formation, despite being controlled by an
externally injected field, arises directly inside the cavity and would be
equally outcoupled through the drop port.

Instead, we achieve pump suppression with a photonic integrated notch filter.
An active ring resonator driven below the lasing threshold can act as such a
filter, where the coupling condition, the quality factor, and the resonance
frequency can be tuned independently via electrical pumping [37]. Such a notch
filter targeting only one laser line can be naturally integrated on the same
chip with the soliton generator. In our implementation, both the soliton
generator and the notch filter are identical racetrack resonators connected
via a bus waveguide. The resonators have nominally equal free spectral range
(FSR) which can be finely tuned by adjusting the pump current in both
resonators and the integrated heaters (Fig. 3a). The soliton is generated by
an off-chip ECL in the first RT QCL driven above threshold, as in Fig. 2, and
the pump line is subsequently filtered by the second RT QCL, operated below
threshold (Fig. 3a). Fig. 3b shows the spectrogram of the output of this chip
while tuning the bias of the notch filter. The filter is close to the critical
coupling, and its FSR is slightly detuned from the intermode spacing of the
soliton frequency comb, so that, relying on a Vernier effect, it filters out
one comb line at each current setting. The filter suppresses the pump line by
45 dB when tuned in resonance with the driving field (Fig. 3c). Fig. 3d shows
the spectra of the output state when the filter is turned off and on,
signifying the highly selective suppression of the driving field, while mildly
affecting the rest of the comb spectrum. A soliton on a high-intensity
background turns into a background-free bright pulse when the filter is
switched on. (Fig. 3e).

The requirement of an additional external laser to drive the formation of
solitons and the reliance on the precise optical alignment of the system would
hinder the practical deployment of these devices in real applications. To
address this issue we further demonstrate a mid-infrared photonic chip where
the pump laser is monolithically integrated with the active racetrack
resonator. Fig. 4a shows such a chip where the drive signal is derived from a
Fabry-Perot (FP) cavity laser that is butt-coupled to the waveguide
directional coupler. Bringing the FP and the RT above their lasing thresholds,
both in single-mode operation, and tuning the current of the heater (HT)
integrated next to RT, we perform the wavelength sweep of the FP through the
RT resonance in the same way as was done with an external drive laser. On a
red-detuned to blue-detuned backward scan over a finite range of detuning the
output intensity exhibits low-frequency fluctuations (Fig. 4b), strongly
reminiscent of the modulation instability comb regime in passive Kerr
resonators [31]. As in Kerr resonators, the associated RF spectrum of the
intermode beat note features multiple frequency tones. On the forward scan
(blue-detuned to red-detuned), soliton formation is associated with the
disappearance of the intensity fluctuations in the laser output, while the
driver laser is still in resonance, and the emergence of a stable microwave
tone at the roundtrip frequency, signifying the high degree of the phase
coherence of the multimode state (Fig. 4b). The spectrum of the soliton state
features tens of lines, all mutually locked, giving rise to a stable bright
pulse emission on top of the CW background (Fig. 4c).

Another superior attribute of the integrated soliton generator, besides its
compactness, is the long-term stability of the soliton state and its ability
to recover upon power cycling. We demonstrate turnkey soliton generation and
perfect state recovery in a sequence of power cycles by switching on and off
the drive currents of the integrated components in a prescribed manner (Fig.
4d). The soliton state initiated in this way remains stable for hours of
continuous operation (Fig. 4e). These features of the fully integrated system
are as well prevalent in Si3N4 Kerr resonators heterogeneously integrated with
III-V semiconductor pump lasers [42].

The demonstrated devices are prototypical components for more elaborate mid-
infrared photonic integrated architectures for short pulse generation. We
showed the first DC-driven, chip-scale pulse generator in this wavelength
range. A straightforward, at the same time, key improvement will be to raise
the peak pulse intensity to the level suitable for supercontinuum generation
in nonlinear crystals, chalcogenide glass fibers, and passive III-V/Si/Ge
waveguides [43]. Power gains can be achieved by improved chip thermal
management (buried heterostructure regrowth, epidown mounting) and by
combining the designs shown in Fig. 3 and Fig. 4 with an additional
integration of a semiconductor optical amplifier at the end of the waveguide
coupler. Although the current demonstration is performed at 8.3 $\upmu$m, our
approach generalizes throughout the mid-infrared ($4-12$ $\upmu$m) by
engineering the QCL active region, and potentially, to shorter wavelengths
using other semiconductor laser gain media — interband cascade, quantum well,
quantum dot, and quantum dash. Soliton mode-locking reliant on a bistability
may be a compelling alternative to passive mode-locking with a saturable
absorber in these laser systems as the absorber induces additional cavity loss
that limits the average optical power and the pulse energy.

The maturity of both semiconductor lasers, active resonators, and nonlinear
Kerr cavities, passive resonators, invites one to consider the shift in the
near future towards hybrid passive-active nonlinear integrated systems. On one
hand, passive Si3N4 nonlinear resonators can be endowed with active functions
— amplification and lasing — by ion implantation [44, 45] or by heterogeneous
and monolithic integration of the III-V active regions and passive waveguides
[23]. On the other hand, as we show here, concepts initially coined in the
context of passive cavities can be extended to resonators with gain, allowing
observation of new phenomena, such as direct on-chip bright pulse generation,
up to now regarded as improbable in QCLs due to the absence of a suitable
mode-locking mechanism.

In this work, all components, including the routing waveguides, are made of
active media, necessitating electrical biasing to reduce optical loss. Further
monolithic integration of III-V passive semiconductor waveguides with QC
active media [46, 47] should enable larger-scale integrated architectures for
linear and nonlinear absorption spectroscopy, absolute time and frequency
metrology in the $4-12$ $\upmu$m window, and free-space optical communications
[48]. For instance, on-chip dual-comb QCL spectrometers is a long-standing
vision proponed multiple times the past [49]. The fact that it is now possible
to have multiple active components on a chip — including lasers above
threshold — biased independently and working simultaneously, each fulfilling
their own function, is another non-trivial and important result of this work.
The intrinsic immunity of the ring resonator based QCL frequency combs to
delayed optical feedback [29] further justifies the on-chip integration of
such components to create functional spectroscopic and metrological systems,
eliminating the need for an optical isolator.

## Methods

Generalized model. Numerical simulations are carried out by integrating the
GLLE model on a regular CPU. The simulation parameters are identical to those
in Ref. [50]. The simulations of both passive and active devices based on the
GLLE, were performed using two different integration schemes, verifying the
robustness of our theoretical predictions. One code is based on the
Exponential Time Stepping scheme [51] with large spatial grids (500-10000
points) and time steps optimized for stability, convergence and numerical
efficiency. The second scheme is based on the modal decomposition of the
electric field $E$ in terms of the longitudinal modes of the cold cavity. The
field is expanded into $2N+1$ modes (the central one and $N$ modes for each
side of the spectrum). Typically, $N=100$.

We further provide an application (MAC OS and Windows) for a custom GUI-
enhanced GLLE solver that we used for rapid prototyping of the longer
parameter sweeps.

Device fabrication and operation. The lasers emit at around 8.3 $\upmu$m and
have a structure consisting of GaInAs/AlInAs layers on an InP substrate. The
active region band structure design is based on a single-phonon continuum
depopulation scheme. The waveguides are dry etched (6 $\upmu$m depth) using
the standard ridge process with optical lithography. The dry etch results in
nearly 90∘ vertical waveguide sidewalls. Waveguide width is 10 $\upmu$m, the
curved section of the racetrack is a semicircle with a radius of 500 $\upmu$m,
and the length of the straight section is $785$ $\upmu$m for the devices in
Fig. 2 and 3 (FSR 18.6 GHz), and $1.5$ mm for the device in Fig. 4 (FSR 13.8
GHz). The air gap in the directional coupler section is 1 $\upmu$m wide. The
facets of the Fabry-Perot pump lasers and of the waveguide couplers are left
as cleaved or are dry etched in places where the optical coupling occurs on
the chip. Laser dies are indium soldered episide up on copper carriers and
individual contact sections are wirebonded to PCBs. All chips are placed on
temperature-stabilized heat sinks at 16∘C. All components (laser resonators,
waveguide couplers, heaters) are individually biased with low-noise current
sources (Wavelength Electronics QCL1500 or QCL2000).

Optical injection setup. The coherent optical drive signal is derived from a
custom-built external cavity tunable QCL (ECL) by DRS Daylight Solutions. The
beam from the ECL is focused onto the WG facet with an aspheric antireflection
coated lens (NA$=0.56$) and the output radiation is collected at the opposite
WG facet with an identical lens. Part of the output is focused directly onto
an MCT photodetector (VIGO photonics, 300 MHz cutoff) for the output
power/transmission measurement. Another part is sent through a Michelson
interferometer onto an MCT photodetector (VIGO photonics, 1 GHz cutoff) for
acquisition of an interferogram and of an optical spectrum. The detector
outputs are sampled at 1.8 GSa/s (output intensity) and 7 MSa/s
(interferogram) using data acquisition (DAQ) and oscilloscope (OSC) modules of
the Zurich Instruments MF and UHF lock-in amplifiers. We control the resonance
detuning and perform the resonance sweep scans by fixig the wavelength of the
ECL and applying a triangular wave modulation with a function generator (SWEEP
box) to the electrical drive current of heater integrated next to the
racetrack resonator. An RF spectrum of the intermode beatnote is acquired by
bringing an RF probe in contact with a laser die. The temporal waveforms are
reconstructed using SWIFTS [52].

## Acknowledgements

This material is based on work supported by the National Science Foundation
under Grant No. ECCS-2221715. T. P. Letsou thanks the Department of Defense
(DoD) through the National Defense Science and Engineering Graduate (NDSEG)
Fellowship Program. N. Opa$\mathrm{\check{c}}$ak and B. Schwarz are supported
by the European Research Council (853014). D. Kazakov and T. Letsou thank
Prof. Kiyoul Yang and Tony Song (Harvard University) for enlightening
discussions.

## Data availability

The data that support the findings of this study are available upon request
from the authors.

## Author contributions

M. Piccardo and D.K. conceived the project. D.K. designed the devices, M.B.,
S.D.C., B.S., fabricated them. D.K. designed the experiments. D.K., T.P.L.,
N.O, and P.R. carried out measurements and analyzed data. L.L.C, M.B., F.P.,
and L.A.L. developed the theory of driven active resonators and carried out
numerical simulations. D.K. and M. Piccardo developed the interactive GLLE
solver. M. Pushkarsky, D.C., and T.D. manufactured an external cavity tunable
laser. All authors participated in the production of the manuscript. B.S., M.
Piccardo, and F.C. supervised the project.

## Competing interests

The authors declare the existence of a financial competing interest. The
Office of Technology Development of Harvard University has begun the process
of filing a patent application based on the materials of this work and is
exploring commercialization opportunities for the presented technology.

## References

  * Chang  _et al._ [2022] L. Chang, S. Liu, and J. E. Bowers, Integrated optical frequency comb technologies, Nature Photonics 16, 95 (2022).
  * Kippenberg  _et al._ [2018] T. J. Kippenberg, A. L. Gaeta, M. Lipson, and M. L. Gorodetsky, Dissipative Kerr solitons in optical microresonators, Science 361, 10.1126/science.aan8083 (2018).
  * Liu  _et al._ [2019] S. Liu, X. Wu, D. Jung, J. C. Norman, M. J. Kennedy, H. K. Tsang, A. C. Gossard, and J. E. Bowers, High-channel-count 20 GHz passively mode-locked quantum dot laser directly grown on Si with 4.1 Tbit/s transmission capacity, Optica 6, 128 (2019).
  * Hillbrand  _et al._ [2019] J. Hillbrand, M. Beiser, A. M. Andrews, H. Detz, A. Schade, R. Weih, S. Höfling, B. Schwarz, and G. Strasser, Picosecond pulses from a mid-infrared interband cascade laser, Optica 6, 1334 (2019).
  * Ren  _et al._ [2023] D. Ren, C. Dong, and D. Burghoff, Integrated nonlinear photonics in the longwave-infrared: A roadmap, MRS Communications , 1 (2023).
  * Marin-Palomo  _et al._ [2019] P. Marin-Palomo, J. N. Kemal, P. Trocha, S. Wolf, K. Merghem, F. Lelarge, A. Ramdane, W. Freude, S. Randel, and C. Koos, Comb-based WDM transmission at 10 Tbit/s using a DC-driven quantum-dash mode-locked laser diode, Optics Express 27, 31110 (2019).
  * Dausinger  _et al._ [2004] F. Dausinger, H. Lubatschowski, and F. Lichtner,  _Femtosecond Technology for Technical and Medical Applications_, edited by F. Dausinger, H. Lubatschowski, and F. Lichtner, Topics in Applied Physics (Springer Berlin Heidelberg, Berlin, Heidelberg, 2004).
  * Riemensberger  _et al._ [2020] J. Riemensberger, A. Lukashchuk, M. Karpov, W. Weng, E. Lucas, J. Liu, and T. J. Kippenberg, Massively parallel coherent laser ranging using a soliton microcomb, Nature 581, 164 (2020).
  * Dudley  _et al._ [2006] J. M. Dudley, G. Genty, and S. Coen, Supercontinuum generation in photonic crystal fiber, Reviews of Modern Physics 78, 1135 (2006).
  * Takamoto  _et al._ [2020] M. Takamoto, I. Ushijima, N. Ohmae, T. Yahagi, K. Kokado, H. Shinkai, and H. Katori, Test of general relativity by a pair of transportable optical lattice clocks, Nature Photonics 14, 411 (2020).
  * Drake  _et al._ [2019] T. E. Drake, T. C. Briles, J. R. Stone, D. T. Spencer, D. R. Carlson, D. D. Hickstein, Q. Li, D. Westly, K. Srinivasan, S. A. Diddams, and S. B. Papp, Terahertz-Rate Kerr-Microresonator Optical Clockwork, Physical Review X 9, 031023 (2019).
  * Spencer  _et al._ [2018] D. T. Spencer, T. Drake, T. C. Briles, J. Stone, L. C. Sinclair, C. Fredrick, Q. Li, D. Westly, B. R. Ilic, A. Bluestone, N. Volet, T. Komljenovic, L. Chang, S. H. Lee, D. Y. Oh, M. G. Suh, K. Y. Yang, M. H. Pfeiffer, T. J. Kippenberg, E. Norberg, L. Theogarajan, K. Vahala, N. R. Newbury, K. Srinivasan, J. E. Bowers, S. A. Diddams, and S. B. Papp, An optical-frequency synthesizer using integrated photonics, Nature 557, 81 (2018).
  * Marin-Palomo  _et al._ [2017] P. Marin-Palomo, J. N. Kemal, M. Karpov, A. Kordts, J. Pfeifle, M. H. Pfeiffer, P. Trocha, S. Wolf, V. Brasch, M. H. Anderson, R. Rosenberger, K. Vijayan, W. Freude, T. J. Kippenberg, and C. Koos, Microresonator-based solitons for massively parallel coherent optical communications, Nature 546, 274 (2017).
  * Cundiff and Mukamel [2013] S. T. Cundiff and S. Mukamel, Optical multidimensional coherent spectroscopy, Physics Today 66, 44 (2013).
  * Picqué and Hänsch [2019] N. Picqué and T. W. Hänsch, Frequency comb spectroscopy, Nature Photonics 13, 146 (2019).
  * Wang  _et al._ [2020] F. Wang, S. Slivken, D. H. Wu, Q. Y. Lu, and M. Razeghi, Continuous wave quantum cascade lasers with 5.6 W output power at room temperature and 41% wall-plug efficiency in cryogenic operation, AIP Advances 10, 10.1063/5.0003318/1021666 (2020).
  * Schwarz  _et al._ [2017] B. Schwarz, C. A. Wang, L. Missaggia, T. S. Mansuripur, P. Chevalier, M. K. Connors, D. McNulty, J. Cederberg, G. Strasser, and F. Capasso, Watt-Level Continuous-Wave Emission from a Bifunctional Quantum Cascade Laser/Detector, ACS Photonics 4, 1225 (2017).
  * Täschler  _et al._ [2021] P. Täschler, M. Bertrand, B. Schneider, M. Singleton, P. Jouy, F. Kapsalidis, M. Beck, and J. Faist, Femtosecond pulses from a mid-infrared quantum cascade laser, Nature Photonics 15, 919 (2021).
  * Columbo  _et al._ [2021] L. Columbo, M. Piccardo, F. Prati, L. A. Lugiato, M. Brambilla, A. Gatti, C. Silvestri, M. Gioannini, N. Opačak, B. Schwarz, and F. Capasso, Unifying Frequency Combs in Active and Passive Cavities: Temporal Solitons in Externally Driven Ring Lasers, Physical Review Letters 126, 173903 (2021).
  * Lugiato and Lefever [1987] L. A. Lugiato and R. Lefever, Spatial Dissipative Structures in Passive Optical Systems, Physical Review Letters 58, 2209 (1987).
  * Lugiato  _et al._ [1988] L. A. Lugiato, C. Oldano, and L. M. Narducci, Cooperative frequency locking and stationary spatial structures in lasers, JOSA B 5, 879 (1988).
  * Piccardo  _et al._ [2020] M. Piccardo, B. Schwarz, D. Kazakov, M. Beiser, N. Opačak, Y. Wang, S. Jha, J. Hillbrand, M. Tamagnone, W. T. Chen, A. Y. Zhu, L. L. Columbo, A. Belyanin, and F. Capasso, Frequency combs induced by phase turbulence, Nature 582, 360 (2020).
  * Xiang  _et al._ [2023] C. Xiang, W. Jin, O. Terra, B. Dong, H. Wang, L. Wu, J. Guo, T. J. Morin, E. Hughes, J. Peters, Q. X. Ji, A. Feshali, M. Paniccia, K. J. Vahala, and J. E. Bowers, 3D integration enables ultralow-noise isolator-free lasers in silicon photonics, Nature 620, 78 (2023).
  * Singleton  _et al._ [2018] M. Singleton, P. Jouy, M. Beck, and J. Faist, Evidence of linear chirp in mid-infrared quantum cascade lasers, Optica 5, 948 (2018).
  * Opačak and Schwarz [2019] N. Opačak and B. Schwarz, Theory of Frequency-Modulated Combs in Lasers with Spatial Hole Burning, Dispersion, and Kerr Nonlinearity, Physical Review Letters 123, 243902 (2019).
  * Burghoff [2020] D. Burghoff, Unraveling the origin of frequency modulated combs using active cavity mean-field theory, Optica 7, 1781 (2020).
  * Täschler  _et al._ [2023] P. Täschler, A. Forrer, M. Bertrand, F. Kapsalidis, M. Beck, and J. Faist, Asynchronous Upconversion Sampling of Frequency Modulated Combs, Laser & Photonics Reviews , 220063 (2023).
  * Meng  _et al._ [2021] B. Meng, M. Singleton, J. Hillbrand, M. Franckié, M. Beck, and J. Faist, Dissipative Kerr solitons in semiconductor ring lasers, Nature Photonics 16, 142 (2021).
  * Opačak  _et al._ [2023] N. Opačak, D. Kazakov, L. L. Columbo, M. Beiser, T. P. Letsou, F. Pilat, M. Brambilla, F. Prati, M. Piccardo, F. Capasso, and B. Schwarz, Nozaki-Bekki optical solitons, arXiv:2304.10796v1 (2023).
  * Leo  _et al._ [2010] F. Leo, S. Coen, P. Kockaert, S. P. Gorza, P. Emplit, and M. Haelterman, Temporal cavity solitons in one-dimensional Kerr media as bits in an all-optical buffer, Nature Photonics 4, 471 (2010).
  * Herr  _et al._ [2013] T. Herr, V. Brasch, J. D. Jost, C. Y. Wang, N. M. Kondratiev, M. L. Gorodetsky, and T. J. Kippenberg, Temporal solitons in optical microresonators, Nature Photonics 8, 145 (2013).
  * Henry [1982] C. H. Henry, Theory of the Linewidth of Semiconductor Lasers, IEEE Journal of Quantum Electronics 18, 259 (1982).
  * Opačak  _et al._ [2021a] N. Opačak, S. D. Cin, J. Hillbrand, and B. Schwarz, Frequency Comb Generation by Bloch Gain Induced Giant Kerr Nonlinearity, Physical Review Letters 127, 093902 (2021a).
  * Opačak  _et al._ [2021b] N. Opačak, F. Pilat, D. Kazakov, S. D. Cin, G. Ramer, B. Lendl, F. Capasso, and B. Schwarz, Spectrally resolved linewidth enhancement factor of a semiconductor frequency comb, Optica 8, 1227 (2021b).
  * Khurgin [2023] J. B. Khurgin, Nonlinear optics from the viewpoint of interaction time, Nature Photonics 17, 545 (2023).
  * Shimoda [1986] K. Shimoda,  _Introduction to Laser Physics_, Springer Series in Optical Sciences (Springer Berlin Heidelberg, Berlin, Heidelberg, 1986).
  * Kazakov  _et al._ [2022] D. Kazakov, T. P. Letsou, M. Beiser, Y. Zhi, N. Opačak, M. Piccardo, B. Schwarz, and F. Capasso, Semiconductor ring laser frequency combs with active directional couplers, arXiv:2206.03379v2 10.48550/arxiv.2206.03379 (2022).
  * Joshi  _et al._ [2016] C. Joshi, J. K. Jang, K. Luke, X. Ji, S. A. Miller, A. Klenner, Y. Okawachi, M. Lipson, and A. L. Gaeta, Thermally controlled comb generation and soliton modelocking in microresonators, Optics Letters 41, 2565 (2016).
  * Englebert  _et al._ [2021] N. Englebert, C. Mas Arabí, P. Parra-Rivas, S.-P. Gorza, and F. Leo, Temporal solitons in a coherently driven active resonator, Nature Photonics 15, 536 (2021).
  * Trocha  _et al._ [2018] P. Trocha, M. Karpov, D. Ganin, M. H. Pfeiffer, A. Kordts, S. Wolf, J. Krockenberger, P. Marin-Palomo, C. Weimann, S. Randel, W. Freude, T. J. Kippenberg, and C. Koos, Ultrafast optical ranging using microresonator soliton frequency combs, Science 359, 887 (2018).
  * Wang  _et al._ [2016] P.-H. Wang, J. A. Jaramillo-Villegas, Y. Xuan, X. Xue, C. Bao, D. E. Leaird, M. Qi, and A. M. Weiner, Intracavity characterization of micro-comb generation in the single-soliton regime, Optics Express 24, 10890 (2016).
  * Xiang  _et al._ [2021] C. Xiang, J. Liu, J. Guo, L. Chang, R. N. Wang, W. Weng, J. Peters, W. Xie, Z. Zhang, J. Riemensberger, J. Selvidge, T. J. Kippenberg, and J. E. Bowers, Laser soliton microcombs heterogeneously integrated on silicon, Science 373, 99 (2021).
  * Brès  _et al._ [2023] C. S. Brès, A. Della Torre, D. Grassani, V. Brasch, C. Grillet, and C. Monat, Supercontinuum in integrated photonics: Generation, applications, challenges, and perspectives, Nanophotonics 12, 1199 (2023).
  * Liu  _et al._ [2022] Y. Liu, Z. Qiu, X. Ji, A. Lukashchuk, J. He, J. Riemensberger, M. Hafermann, R. Ning Wang, J. Liu, C. Ronning, and T. J. Kippenberg, A photonic integrated circuit-based erbium-doped amplifier, Science 376, 1309 (2022).
  * Gaafar  _et al._ [2023] M. A. Gaafar, M. Ludwig, K. Wang, T. Wildi, T. Voumard, M. Sinobad, J. Lorenzen, H. Francis, S. Zhang, T. Bi, P. DeľHaye, M. Geiselmann, N. Singh, F. X. Kärtner, S. M. Garcia-Blanco, and T. Herr, Femtosecond pulse amplification on a chip, arXiv:2311.04758v1 (2023).
  * Jung  _et al._ [2019] S. Jung, D. Palaferri, K. Zhang, F. Xie, Y. Okuno, C. Pinzone, K. Lascola, and M. A. Belkin, Homogeneous photonic integration of mid-infrared quantum cascade lasers with low-loss passive waveguides on an InP platform, Optica 6, 1023 (2019).
  * Wang  _et al._ [2022] R. Wang, P. Täschler, Z. Wang, E. Gini, M. Beck, and J. Faist, Monolithic Integration of Mid-Infrared Quantum Cascade Lasers and Frequency Combs with Passive Waveguides, ACS Photonics 9, 426 (2022).
  * Dely  _et al._ [2022] H. Dely, T. Bonazzi, O. Spitz, E. Rodriguez, D. Gacemi, Y. Todorov, K. Pantzas, G. Beaudoin, I. Sagnes, L. Li, A. Giles Davies, E. H. Linfield, F. Grillot, A. Vasanelli, and C. Sirtori, 10 Gbit/s Free Space Data Transmission at 9 µm Wavelength With Unipolar Quantum Optoelectronics, Laser & Photonics Reviews 16, 2100414 (2022).
  * Villares  _et al._ [2015] G. Villares, J. Wolf, D. Kazakov, M. J. Süess, A. Hugi, M. Beck, and J. Faist, On-chip dual-comb based on quantum cascade laser frequency combs, Applied Physics Letters 107, 251104 (2015).
  * Prati  _et al._ [2021] F. Prati, M. Brambilla, M. Piccardo, L. L. Columbo, C. Silvestri, M. Gioannini, A. Gatti, L. A. Lugiato, and F. Capasso, Soliton dynamics of ring quantum cascade lasers with injected signal, Nanophotonics 10, 195 (2021).
  * Cox and Matthews [2002] S. M. Cox and P. C. Matthews, Exponential Time Differencing for Stiff Systems, Journal of Computational Physics 176, 430 (2002).
  * Burghoff  _et al._ [2020] D. Burghoff, D. Ren, and Z. Han, Sensitivity of SWIFT spectroscopy, Optics Express 28, 6002 (2020).