# Enhancing Short-Term Wind Speed Forecasting using Graph Attention and
Frequency-Enhanced Mechanisms

Hao Liu, Huimin Ma, Tianyu Hu H. Liu, T. Hu and H. Ma are with the School of
Computer and Communication Engineering, University of Science and Technology
Beijing, Beijing, China.

###### Abstract

The safe and stable operation of power systems is greatly challenged by the
high variability and randomness of wind power in large-scale wind-power-
integrated grids. Wind power forecasting is an effective solution to tackle
this issue, with wind speed forecasting being an essential aspect. In this
paper, a Graph-attentive Frequency-enhanced Spatial-Temporal Wind Speed
Forecasting model based on graph attention and frequency-enhanced mechanisms,
i.e., GFST-WSF, is proposed to improve the accuracy of short-term wind speed
forecasting. The GFST-WSF comprises a Transformer architecture for temporal
feature extraction and a Graph Attention Network (GAT) for spatial feature
extraction. The GAT is specifically designed to capture the complex spatial
dependencies among wind speed stations to effectively aggregate information
from neighboring nodes in the graph, thus enhancing the spatial representation
of the data. To model the time lag in wind speed correlation between adjacent
wind farms caused by geographical factors, a dynamic complex adjacency matrix
is formulated and utilized by the GAT. Benefiting from the effective spatio-
temporal feature extraction and the deep architecture of the Transformer, the
GFST-WSF outperforms other baselines in wind speed forecasting for the 6-24
hours ahead forecast horizon in case studies.

###### Index Terms:

Wind speed forecast, deep learning, spatial-temporal correlations,
Transformer, GAT.

##  I Introduction

In recent years, the rapid development of modern industry and the exponential
growth of population have caused environmental pollution and global energy
crisis to become increasingly serious. As a result, renewable energy has
gained widespread attention for its potential to solve environmental and
energy problems. Among clean renewable energy sources, wind power has seen
overwhelming growth over the past decade, with global wind power capacity
reaching 837 GW in 2021 [1]. However, the variability of wind speed can lead
to intermittent power generation, which poses challenges to the safety and
stability of smart grid energy management systems. Studies have shown that
wind power forecasting is one of the most cost-effective and efficient methods
to minimize this risk [2, 3, 4]. Therefore, accurate wind power forecasting is
urgently needed, as well as wind speed forecasting, which actually serves as
the foundation for wind power forecasting [5, 6].

Numerous wind speed forecasting methods have been proposed, including
physical, statistical, machine learning, and deep learning methods. Physical
methods mostly rely on Numeric Weather Prediction (NWP), which uses weather
data such as temperatures and pressure to solve complex mathematical and
physical models and forecast wind speed [7]. Statistical methods, on the other
hand, construct a linear or nonlinear function from historical data to
forecast future wind speed values [8], e.g., Auto-Regressive Moving Average
(ARMA) [9], Auto-Regressive Integrated Moving Average (ARIMA) [10], and
Eupport Vector Regression (SVR) [11], etc.

Artificial intelligence technologies have led to the emergence and development
of numerous machine learning and deep learning models for wind speed
forecasting. AI-based short-term wind power forecasting methods include
Extreme Learning Machines (ELM) [12], Light Gradient Boosting Machine
(LightGBM) [13], Artificial Neural Network (ANN) [14], Convolutional Neural
Network (CNN) [15], Recurrent Neural Network (RNN) [16], Long Short-term
Memory Network (LSTM) [17], among others. However, these methods have mainly
relied on historical data from individual wind farms, hardly considering
cross-farm spatial correlation [18]. As multiple wind farms close to each
other are typically within the same wind belt, wind speeds are strongly
correlated spatially. As a result, several studies considering both temporal
and spatial correlation between adjacent wind fields have emerged in recent
years, and have further improved forecasting performance [19, 20, 21, 22].

In [19], a spatial vector consisting of wind speed values from multiple nodes
was used as the spatial feature input to the forecasting model, while the
temporal features were extracted by LSTM. [20] proposed a model for wind speed
forecasting with spatio-temporal correlation, integrating CNN and a multi-
layer perceptron (MLP). Here, the spatial features were extracted by CNN, and
the MLP captured the temporal dependencies among these extracted spatial
features. In [21], a deep architecture named the Predictive Spatio-Temporal
Network (PSTN) was used for wind speed forecasting, integrating CNN and LSTM.
Similar to the CNN-MLP model, the temporal feature extraction module of this
model was changed from MLP to LSTM. Finally, in [22], a model was proposed to
forecast the wind farm cluster by combining Graph Convolution Networks (GCN)
with LSTM. The performance improvements observed in these methods, whether
using CNN or Graph Neural Network (GNN), demonstrate the effectiveness of
utilizing spatial correlation for wind speed forecasting. Therefore, we
utilize the GAT[23] to extract spatial correlation features between
neighboring wind farms, aiming to further enhance the forecasting performance
of our model.

Recently, the Transformer architecture has gained attention in deep learning
for its ability to model long-range dependencies, e.g., ChatGPT, which makes
it suitable for time series modeling [24]. In this paper, we leverage the
high-performing Frequency Enhanced Decomposed Transformer (FEDformer) [25] to
model the auto-correlation of wind speed series. However, we observe that
FEDformer primarily operates on frequency domain information obtained through
Fourier transformation and lacks traditional processing of time-domain
information. To address this limitation, we introduce a multi-head self-
attention mechanism that allows us to effectively handle time-domain
information. By incorporating this mechanism, we enhance the model’s ability
to capture temporal features and improve its forecasting performance.

In this study, we analyzed wind speed data from 25 wind farms to determine
their correlation. We found that the spatial distribution of wind farms
varies, resulting in potential time lags in the correlation of wind speeds.
Additionally, wind speeds can vary with different wind directions. To address
these challenges, we propose a dynamic complex adjacency matrix that can
simultaneously consider the spatial correlation and time lags between wind
speeds from different wind farms. Furthermore, it can reflect real-time
changes in the correlation between wind speeds, capturing their evolving
relationship.

To capture the spatio-temporal features among neighboring wind farms, we
integrate GAT and FEDformer models and incorporate a multi-head self-attention
mechanism to effectively process time-domain information. As a result, we
propose the GFST-WSF model.

The main contributions of this paper are summarized as follows:

(a) We propose a complex adjacency matrix to represent the spatio-temporal
correlation among neighboring wind farms. Each element of the matrix is a
complex number, where the real part represents the strength of the
correlation, and the imaginary part represents the time difference of the
correlation’s occurrence. This approach provides a more comprehensive
representation of the spatio-temporal correlation, leading to better
forecasting performance in the proposed GFST-WSF model.

(b) We enhance the GFST-WSF model with the GAT to effectively utilize the
complex adjacency matrix and capture intricate spatial dependencies among
neighboring wind farms’ wind speed. This approach improves the model’s ability
to extract spatial features, which is crucial for wind speed forecasting
accuracy.

(c) We introduce the FEDformer architecture, which leverages a frequency-
enhanced mechanism to improve the temporal information extraction process in
wind speed forecasting. We have further improved the FEDformer by adding the
multi-headed attention mechanism to enhance its ability to capture long and
short-term dependencies in the wind speed series. The deep architecture and
multi-head attention of the GFST-WSF model enables it to effectively extract
the temporal correlation within the wind speed from neighboring wind farms.

The rest of this paper is organized as follows: Section II presents the
problem description; Section III presents the framework of GFST-WSF; Section
IV presents the case studies; conclusions are drawn in Section V.

##  II Problem Description

We focus on the problem of short-term spatial-temporal wind speed forecasting
for multiple wind farms. The meteorological data (including wind speed, wind
direction, air pressure, temperature, etc.) of the $i$-th wind farm at time
period $t$ can be denoted as $x_{t}^{i}$. Thus, the data vector of all the
wind farms at time $t$ is
$\mathbf{x}_{t}=\\{x_{t}^{1},x_{t}^{2},x_{t}^{3},...,x_{t}^{n}\\}$, where $n$
represents the total number of wind farms. We use multiple variables to
forecast the wind velocity vector $p_{t+r}$ of a single wind farm, where $r\in
Z^{+}$ represents the lead time.

Therefore, the forecasting model can be formulated as follows:

$f(\mathbf{x}_{t},\mathbf{x}_{t-1},\mathbf{x}_{t-2},...\mathbf{x}_{t-(m-1)},\theta,\mathcal{G})=\hat{p}_{t+r}$
(1)

where $f(\cdot)$ represents mapping function, $\theta$ represents the
parameters of $f(\cdot)$, $\mathcal{G}$ represents the spatial relationships
among wind farms, $\hat{p}_{t+r}$ is the forecasted value of $p_{t+r}$, $m$ is
the historical period for forecasting.

Figure 1: The GFST-WSF Structure. The Self-Attention Block utilizes the
original multi-headed attention mechanism in Transformer. The GAT Block is
composed of the time series transform module (TST) and the graph attention
network (GAT). The Frequency Enhanced Block (FEB) and Frequency Enhanced
Attention (FEA) are utilized to perform representation learning in frequency
domain. The series decomposition blocks (SeriesDecomp) are employed to
decompose the series into trend-cyclical and seasonal parts.

##  III The Framework of GFST-WSF

Inspired by recent innovations in the field of time series forecasting, GFST-
WSF adopts an encoder-decoder structure, which includes graph neural network
encoding, multi-head attention from the original Transformer, sequence
decomposition using average pooling, frequency-enhanced mechanism from
FEDformer for sequence-level connections, and feedforward networks, as shown
in Fig.1.

Figure 2: GAT Block. Constructing a complex adjacency matrix based on
historical data (illustrated by a simple example in the figure), where the
upper part represents the real components and the lower part represents the
imaginary components.

###  III-A GAT Block

The GAT Block, illustrated in Fig.2, consists of two components: a GAT and a
Time Series Transform module (TST). Initially, the dynamic complex adjacency
matrix $\Lambda$ is decomposed into a real part matrix $\mathcal{A}$ and an
imaginary part matrix $\mathcal{B}$. The elements of $\mathcal{A}$ reflect the
correlation between wind farms, while the elements of $\mathcal{B}$ represent
the time difference between them. In the TST module, the original time series
features are transformed according to the elements of $\mathcal{B}$. Then, the
updated time series features are used as input to the GAT, with the adjacency
matrix of GAT being $\mathcal{A}$. Finally, the learned spatio-temporal
information is utilized as input to the encoder.

####  III-A1 Building a Dynamic Complex Adjacency Matrix

In this study, we propose a novel approach to represent the correlation
between all wind farms in a given graph. Traditionally, the correlation
between wind farms is represented using a simple adjacency matrix, where the
nodes are the wind farms and the edges represent the correlation between them.
However, due to the variability of wind speed, distance between wind farms,
and constantly changing wind direction, the correlation between adjacent wind
farms may exist at different time lags. To address this issue, we introduce a
dynamic complex adjacency matrix denoted as $\Lambda$ to represent the time
difference in the correlation between wind farms. This matrix contains complex
values of the form $a+bi$, as illustrated in Fig.2. The real part of the
complex value, $a$, represents the correlation between wind farms. If $a>0$,
it indicates a positive correlation between wind farms. On the other hand, the
imaginary part of the complex value, $b$, represents the time difference in
the correlation between wind farms. If there is no correlation between wind
farms, the value of $b$ is invalid. The complex adjacency matrix is defined as
follows:

$\Lambda(s_{i},s_{j})=\begin{cases}a+bi&\text{if}\ \ a>0\\\
0&\text{otherwise}\end{cases}$ (2)

And the definition of the parameter $a$ is given as follows:

$\begin{aligned} a=\begin{cases}1&if\ \ \ F(s_{i},s_{j})>\beta\\\
0&otherwise\end{cases}\end{aligned}$

(3)

where $\beta$ is the threshold of the correlation coefficient, $s_{i}$ and
$s_{j}$ represent the wind speed series of the $i$-th and $j$-th wind farm,
and $F(\cdot)$ is a function that calculates the correlation between the two
wind farms. In our study, we use the Pearson correlation coefficient, which is
defined as follows:

$F(s_{i},s_{j})=\frac{\sum{(s_{i}-\mathaccentV{bar}0\symboldoperators
16{s}_{i})(s_{j}-\mathaccentV{bar}0\symboldoperators
16{s}_{j})}}{\sqrt{\sum{(s_{i}-\mathaccentV{bar}0\symboldoperators
16{s}_{i})^{2}}\sum{(s_{j}-\mathaccentV{bar}0\symboldoperators
16{s}_{j})^{2}}}}$ (4)

To calculate the time difference $b$ between any two wind farms, we use the
latitude and longitude information in the dataset to determine the azimuth
information. We also apply the Haversine formula [26] to calculate the
distance between any two points. Based on these values, we factor in the
azimuth, wind direction, and distance to determine the time difference $b$.

####  III-A2 Graph Attention Networks

Wind is a rapidly propagating air flow that affects a large area, resulting in
intrinsic correlations between the wind speeds of adjacent wind farms. To
capture these relationships, we model each wind farm as a graph node and the
meteorological information series for the wind farm as the node’s
characteristic quantity. The correlation between wind farms is captured as the
attribute of the graph’s edges.

GNN [27] refers to a general class of models that apply neural networks to
graphs. These models can be categorized into different types based on their
underlying techniques. However, as wind direction changes over time, the
correlation between adjacent wind farms can also change. To account for this
dynamic relationship, we employ GAT [23], which leverage masked self-
attentional layers to overcome the limitations of prior methods, such as GCN
and Graph Sample and Aggregate (GraphSAGE), that cannot handle dynamic graph
relationships.

The input of GAT is a set of node features,
$\mathbf{h}=\\{\vec{h}_{1},\vec{h}_{2},...,\vec{h}_{n}\\}$, where $n$ is the
number of nodes. The importance of node $j$’s features to node $i$ is defined
as:

$\begin{aligned}
e_{ij}=a(\mathbf{W}\vec{h}_{i},\mathbf{W}\vec{h}_{j})\end{aligned}$

(5)

where $\mathbf{W}$ is a weight matrix, and the attention mechanism is a
single-layer feedforward neural network. Meanwhile, use the softmax function
to normalize them:

$\begin{aligned}
\alpha_{ij}=softmax(e_{ij})=\frac{exp(e_{ij})}{\sum_{k\in\mathcal{N}_{i}}exp(e_{ik})}\end{aligned}$

(6)

Furthermore, GAT introduces the idea of multi-headed attention mechanism into
graph neural networks. Specifically, perform $k$ separate attention
calculations and then concatenate them:

$\begin{aligned}
\vec{h}_{i}^{\prime}=\|_{k=1}^{K}\sigma\left(\sum_{j\in\mathcal{N}_{i}}\alpha_{ij}^{k}\omega^{k}\vec{h}_{j}\right)\end{aligned}$

(7)

where $\|$ is the concatenation operation and $\omega^{k}$ is the
corresponding input linear transformation’s weight matrix.

###  III-B Encoder

The encoder comprises multiple layers, with the output of the $l$-th layer
denoted as $\mathcal{X}_{en}^{l}=Encoder(\mathcal{X}_{en}^{l-1})$, where
$l\in\\{1,2,..,N\\}$. The input to the encoder is $\mathcal{X}_{en}^{0}$,
which is the historical series that has been embedded using a graph attention
network and one-dimensional convolution. The spatio-temporal information
learned from this embedding is fed to the encoder to enable learning of the
sequence through frequency enhanced blocks, which are then decomposed using
average pooling. To capture interdependent features in the sequence and
increase model variability, multi-head attention is used. The overall learning
process is formulated as follows:

$\displaystyle\mathcal{X}_{en}^{0}$
$\displaystyle=Conv1D(GAT(\mathcal{X},\mathcal{G}))$ (8)
$\displaystyle\mathcal{I}_{en}^{l}$
$\displaystyle=Attention(\mathcal{X}_{en}^{l-1})+\mathcal{X}_{en}^{l-1}$ (9)
$\displaystyle\mathcal{S}_{en}^{l,1},\\_$
$\displaystyle=SeriesDecomp(FEB(\mathcal{I}_{en}^{l})+\mathcal{I}_{en}^{l})$
(10) $\displaystyle\mathcal{S}_{en}^{l,2},\\_$
$\displaystyle=SeriesDecomp(FeedForward(\mathcal{S}_{en}^{l,1})+\mathcal{S}_{en}^{l,1})$
(11) $\displaystyle\mathcal{X}_{en}^{l}$
$\displaystyle=\mathcal{S}_{en}^{l,2}$ (12)

where $\mathcal{X}$ represents the meteorological information of all wind
farms, $\mathcal{G}$ represents the graph network code of the sites,
$\mathcal{I}_{en}^{l}$ is the result of multi-head attention, and
$\mathcal{S}_{en}^{l,i},i\in{1,2}$ represents the seasonal component after the
$i$-th series decomposition block in the $l$-th layer respectively. The
frequency enhanced block (FEB) module is implemented using a Discrete Fourier
Transform (DFT) mechanism.

####  III-B1 Frequency Enhanced Block with Fourier Transform (FEB)

The structure presented in this paper leverages the Discrete Fourier Transform
(DFT), which was first introduced in [25]. By employing the Fast Fourier
Transform (FFT) and randomly selecting a fixed number of Fourier components,
the computational complexity of the structure can be reduced to $O(N)$, making
it more efficient.

The FEB block takes an input $\bm{x}\in\mathbb{R}^{N\times D}$, which is first
linearly projected using $\bm{w}\in\mathbb{R}^{D\times D}$, and then converted
to the frequency domain using the Fourier transform. Specifically, the input
is transformed as follows: $\bm{q}=\bm{x}\cdot\bm{w},\ \
\bm{Q}=\mathcal{F}(\bm{q})$, where $\mathcal{F}$ represents the Fourier
transform, and $\bm{Q}\in\mathbb{C}^{N\times D}$. In the frequency domain,
only a randomly selected subset of $M(M<<N)$ modes are retained, resulting in
$\bm{\widetilde{Q}}=Select(\bm{Q})$, where
$\bm{\widetilde{Q}}\in\mathbb{C}^{M\times D}$. Thus, the FEB is defined as

$\begin{aligned}
FEB(\bm{q})=\mathcal{F}^{-1}(Padding(\widetilde{\bm{Q}}\odot\bm{R}))\end{aligned}$

(13)

where $\bm{R}\in\mathbb{C}^{D\times D\times M}$ is a parameterized kernel
initialized randomly, $\mathcal{F}^{-1}$ denotes the inverse Fourier
transform. The operator $\odot$ is the complex matrix multiplication. I.e.,
the result is zero-padded to $\mathbb{C}^{N\times D}$ and converted back to
the time dimension by $\mathcal{F}^{-1}$.

####  III-B2 Series Decomp

In order to model the complex temporal patterns in wind speed sequences, the
decomposition approach is adopted to divide the time series into two
components: the trend-cyclical component and the seasonal component, which
respectively reflect the long-term progress and periodicity of the time
series. In practice, the moving average method is used to smooth out the
cyclical fluctuations and highlight the long-term trend. Specifically, for an
input sequence $\mathcal{X}$ of length $L$, the decomposition process is as
follows:

$\displaystyle\mathcal{X}_{t}$ $\displaystyle=AvgPool(Padding(\mathcal{X}))$
(14) $\displaystyle\mathcal{X}_{s}$
$\displaystyle=\mathcal{X}-\mathcal{X}_{t}$ (15)

where $\mathcal{X}_{t},\mathcal{X}_{s}$ denote the trend-cyclical and the
seasonal part respectively. We use the $AvgPool(\textperiodcentered)$
operation with padding for moving average and summarize the above equations
using $\mathcal{X}_{t},\mathcal{X}_{s}=SeriesDecomp(\mathcal{X})$.

###  III-C Decoder

The decoder also adopts a multilayer structure as :
$\mathcal{X}_{de}^{l},\mathcal{T}_{de}^{l}=Decoder(\mathcal{X}_{de}^{l-1},\mathcal{T}_{de}^{l-1})$,
where $l\in\\{1,2,...,M\\}$ is the output of $l$-th decoder layer. The decoder
takes two inputs: $\mathcal{X}_{de}$, which is the spatio-temporal embedding
obtained from the seasonal initialization term through graph attention and
one-dimensional convolution, and $\mathcal{T}_{de}$, which is the direct input
of the trend initialization term. The first layer of the decoder,
$\mathcal{X}_{de}^{0}$, is obtained by applying a one-dimensional convolution
to the output of a graph attention layer that takes as input the concatenation
of the average-pooled $\mathcal{X}$ and the seasonal initialization term
$\mathcal{S}_{init}$, with attention weights given by the graph $\mathcal{G}$.
In each subsequent layer $l$, the input $\mathcal{X}_{de}^{l-1}$ is refined by
an attention mechanism and added to the output of the encoder for the
following decoding. The trend component $\mathcal{T}_{de}$ is accumulated
across the layers to improve the model’s inference ability. Specifically, we
define the decomposition blocks as follows:

$\displaystyle\mathcal{S}_{de}^{0}$
$\displaystyle=Concat(SeriesDecomp(\mathcal{X}),\mathcal{S}_{init})$ (16)
$\displaystyle\mathcal{X}_{de}^{0}$
$\displaystyle=Conv1D(GAT(\mathcal{S}_{de}^{0},\mathcal{G}))$ (17)
$\displaystyle\mathcal{I}_{de}^{l}$
$\displaystyle=Attention(\mathcal{X}_{de}^{l-1})+\mathcal{X}_{de}^{l-1}$ (18)
$\displaystyle\mathcal{T}_{de}^{0}$
$\displaystyle=Concat(SeriesDecomp(\mathcal{X}),\mathcal{T}_{init})$ (19)
$\displaystyle\mathcal{S}_{de}^{l,1},\mathcal{T}_{de}^{l,1}$
$\displaystyle=SeriesDecomp(FEB(\mathcal{I}_{de}^{l})+\mathcal{I}_{de}^{l})$
(20) $\displaystyle\mathcal{S}_{de}^{l,2},\mathcal{T}_{de}^{l,2}$
$\displaystyle=SeriesDecomp(FEA(\mathcal{S}_{de}^{l,1},\mathcal{X}_{en}^{N})+\mathcal{S}_{de}^{l,1})$
(21) $\displaystyle\mathcal{S}_{de}^{l,3},\mathcal{T}_{de}^{l,3}$
$\displaystyle=SeriesDecomp(FeedForward(\mathcal{S}_{de}^{l,2})+\mathcal{S}_{de}^{l,2})$
(22) $\displaystyle\mathcal{T}_{de}^{l}$
$\displaystyle=\mathcal{T}_{de}^{l-1}+\mathcal{W}_{l,1}\cdot\mathcal{T}_{de}^{l,1}+\mathcal{W}_{l,2}\cdot\mathcal{T}_{de}^{l,2}+\mathcal{W}_{l,3}\cdot\mathcal{T}_{de}^{l,3}$
(23) $\displaystyle\mathcal{X}_{de}^{l}$
$\displaystyle=\mathcal{T}_{de}^{l,3}$ (24)

where $\mathcal{S}_{init}$ and $\mathcal{T}_{init}$ represent the
initialization of seasonal and trend terms respectively. After the $i$-th
series decomposition block in the $l$-th layer, the seasonal and trend
components are represented by $\mathcal{S}_{de}^{l,i}$ and
$\mathcal{T}_{de}^{l,i}$, where $i\in\\{1,2,3\\}$. The linear projector for
the $i$-th extracted trend $\mathcal{T}_{de}^{l,i}$ is denoted by
$\mathcal{W}_{l,i}$. The frequency enhanced attention (FEA) is similar to FEB
in that it uses DFT projection with an attention design.

Frequency Enhanced Attention with Fourier Transform(FEA) FEA takes inputs that
are similar to those of the original Transformer, including queries
$\bm{q}\in\mathbb{R}^{L\times D}$, keys $\bm{k}\in\mathbb{R}^{L\times D}$, and
values $\bm{v}\in\mathbb{R}^{L\times D}$. The queries are sourced from the
decoder, while the keys and values come from the encoder :
$\bm{q}=\bm{x}_{de}\cdot\bm{w}_{q},\ \bm{k}=\bm{x}_{en}\cdot\bm{w}_{k},\
\bm{v}=\bm{x}_{en}\cdot\bm{w}_{v}$ , where
$\bm{w}_{q},\bm{w}_{k},\bm{w}_{v}\in\mathbb{R}^{D\times D}$. Then the
canonical attention can be defined as

$\begin{aligned}
Atten(\bm{q},\bm{k},\bm{v})=Softmax(\frac{\bm{q}\bm{k}^{T}}{\sqrt{d_{q}}})\bm{v}\end{aligned}$

(25)

where $d_{q}$ represents the length of $\bm{q}$. In FEA, the queries, keys,
and values ($\bm{q},\bm{k}$, and $\bm{v}$) are first transformed into the
frequency domain using the Fourier transform, after which a subset of $M$
modes is randomly selected in the frequency domain. Next, a similar cross-
attention mechanism is applied in the frequency domain. The FEA can be defined
as follows:

$\displaystyle FEA(\bm{q},\bm{k},\bm{v})$
$\displaystyle=\mathcal{F}^{-1}(Padding(\sigma(\widetilde{\bm{Q}}\cdot\widetilde{\bm{K}})\cdot\widetilde{\bm{V}}))$
(26)

where $\sigma$ is the activation function,
$\widetilde{\bm{Q}}=Select(\mathcal{F}(\bm{q}))$,
$\widetilde{\bm{K}}=Select(\mathcal{F}(\bm{k})),$ and
$\widetilde{\bm{V}}=Select(\mathcal{F}(\bm{v}))$. The result of
$\sigma(\widetilde{\bm{Q}}\cdot\widetilde{\bm{K}})\cdot\widetilde{\bm{V}}$
also need to be zero-padded to $\mathbb{C}^{L\times D}$ before converting back
to the time dimension.

##  IV Case Studies

###  IV-A Evaluation Metrics

The Mean Squared Error (MSE) and normalized Mean Absolute Error (MAE) are
commonly used as evaluation metrics in the task of forecasting wind speed at a
single site. The metrics are defined as follows:

$\displaystyle MSE=\frac{1}{T}\sum_{t=1}^{T}(p_{t+r}-\hat{p}_{t+r})^{2}$ (27)
$\displaystyle MAE=\frac{1}{T}\sum_{t=1}^{T}|p_{t+r}-\hat{p}_{t+r}|$ (28)

Figure 3: The geographic distribution of the wind farms. These wind farms are
located in southern India and have been clearly marked on the map.

###  IV-B Data Description

The Wind Integration National Dataset (WIND) Toolkit [28] is a software
package developed by the National Renewable Energy Laboratory (NREL) in the
United States. It provides researchers with access to a wide range of wind
power data. For this study, we utilized the India Wind Dataset from WIND,
which was developed in part through the India Renewable Integration Study.
This dataset contains wind speed, wind direction, temperature, and pressure at
heights of 40m, 80m, 100m, and 120m above the ground. The data has a spatial
resolution of 3 km and a temporal resolution of 5 minutes. From this dataset,
we selected 25 sites in Tamil Nadu, India for short-term wind speed
forecasting, as shown in Fig.3. The annual dataset for 2014 was chosen and the
time resolution was adjusted to 15 minutes. The wind farm being forecasted is
located at the center of the selected sites and is the largest wind farm in
southern India, named the Muppandal Wind Farm (the green point). The distance
between sites ranges from 5 kilometers to 40 kilometers, as shown in the
Fig.3.

###  IV-C Comparison on Test Set

To verify the superiority of the GFST-WSF, we employed state-of-the-art time
series forecasting methods and deep learning methods as baselines. There are
following:

1) The Persistence method (Pers.): The Persistence method forecasts the next
time step by taking the current time step’s observed value as the forecasted
value, i.e., all forecasted values are equal to the last observed value at the
known time step.

2) LightGBM: LightGBM is a mainstream implementation of the Gradient Boosting
Decision Tree (GBDT) algorithm, which utilizes iteratively trained weak
classifiers (decision trees) to obtain an optimized model. It is commonly used
for forecasting tasks.

3) Long Short-Term Memory (LSTM): LSTM, as a deep learning algorithm used for
processing sequence data, features memory cells and gate mechanisms,
exhibiting strong abilities in modeling long-term dependencies and effectively
addressing the issues of vanishing and exploding gradients.

4) Transformer: Transformer is a neural network architecture based on self-
attention mechanism used for sequence-to-sequence tasks and is one of the
state-of-the-art models in this field.

5) FEDformer: FEDformer is a variant of the Transformer model that introduces
frequency domain enhancement mechanisms, and currently demonstrates excellent
performance in time series forecasting tasks.

(a) Results on test set (6-hour ahead)

(b) Results on test set (12-hour ahead)

(c) Results on test set (24-hour ahead)

Figure 4: Comparison between GFST-WSF and baselines.

(a) Results on test set (6-hour ahead)

(b) Results on test set (12-hour ahead)

(c) Results on test set (24-hour ahead)

Figure 5: Comparison of results from ablation experiments. TABLE I: MSE of
wind speed in the next 6h, 12h and 24h by different models.

Models | Forecasting Horizon
---|---
24(6-hour) | 48(12-hour) | 96(24-hour)
Pers. | $0.215$ | $0.378$ | $0.415$
LightGBM | $0.164$ | $0.243$ | $0.284$
LSTM | $0.175$ | $0.272$ | $0.366$
Transformer | $0.165$ | $0.236$ | $0.339$
FEDformer | $0.165$ | $0.225$ | $0.244$
GFST-WSF | $\mathbf{0.141}$ | $\mathbf{0.197}$ | $\mathbf{0.226}$

* All of the data has been normalized.

* The experimental results are the mean of multiple experiments.

TABLE II: MAE of wind speed in the next 6h, 12h and 24h by different models.

Models | Forecasting Horizon
---|---
24(6-hour) | 48(12-hour) | 96(24-hour)
Pers. | $0.337$ | $0.466$ | $0.499$
LightGBM | $0.299$ | $0.379$ | $0.414$
LSTM | $0.314$ | $0.405$ | $0.489$
Transformer | $0.317$ | $0.393$ | $0.473$
FEDformer | $0.315$ | $0.373$ | $0.390$
GFST-WSF | $\mathbf{0.282}$ | $\mathbf{0.344}$ | $\mathbf{0.374}$

* All of the data has been normalized.

* The experimental results are the mean of multiple experiments.

In TABLE I and TABLE II, we present a comparison of the performance of our
GFST-WSF model with other time series forecasting models. The results clearly
demonstrate that GFST-WSF outperforms all considered methods in terms of the
lowest MAE and MSE values. As a commonly used machine learning method,
LightGBM performs well in short-term forecasting but exhibits poor performance
in long-term forecasting. Among the baselines, LSTM performed the worst in
terms of forecasting performance, with its error increasing sharply as the
length of the forecasting horizon increased. However, the Transformer model
showed improved performance, owing to its effective capture of the attention
mechanism on the long-term and short-term dependence of time series. Building
upon the Transformer architecture, FEDformer achieved significant performance
improvements in forecasting long-term time series. Specifically, in the
12-hour forecasting task, FEDformer lowered MSE and MAE by approximately 5%,
while in the 24-hour forecasting task, it demonstrated a 28% reduction in MSE
and a 17% reduction in MAE. This is because FEDformer replaced the multi-head
self-attention mechanism in Transformer with a frequency-enhanced mechanism,
and applied cross-attention to the frequency-domain information. Additionally,
FEDformer utilized a temporal decomposition module, which improved its
performance in long sequence forecasting.

GFST-WSF is a model based on FEDformer, which combines GAT module to extract
spatial features and uses multi-head attention mechanism to extract temporal
information, thus further improving the forecasting performance. At the same
time, GFST-WSF uses dynamic complex adjacency matrix to model the time lag
relationship between adjacent wind farms, in order to better capture the wind
speed correlation. Through these improvements, GFST-WSF can effectively
extract various wind speed features and capture complex patterns in wind speed
data, thereby performing well in short-term wind speed forecasting tasks.
Compared with other baseline models, GFST-WSF has stronger modeling and
representation capabilities, and therefore has better application prospects in
ensuring the safe and stable operation of wind power integrated grids.

The experimental results on the test set are shown in Fig.4, where Fig.4a,
Fig.4b, and Fig.4c represent the experimental results of all comparative
models in forecasting future 6, 12, and 24 hours, respectively. It can be
clearly observed that GFST-WSF is the closest to the ground truth in the
results of different forecasting lengths. Especially as the forecasting length
increases, the performance advantage of GFST-WSF becomes more pronounced.

###  IV-D Ablation experiments

To verify the effectiveness and necessity of each module in GFST-WSF, we also
conducted an additional ablation experiment, and the results are listed in
TABLE III. After removing the multi-head self-attention, the new model is
called GAT-FEDformer; then remove the GAT Block, and the model becomes the
original FEDformer. TABLE III shows the following.

1) FEDformer has shown remarkable performance in forecasting longer time
series, with its frequency-enhanced mechanism playing a crucial role.

2) GAT-FEDformer outperforms FEDformer by 5%-12% in forecasting wind speed at
different lengths, indicating the effectiveness of the GAT and the designed
dynamic complex adjacency matrix in capturing spatial features.

3) GFST-WSF surpasses GAT-FEDformer in terms of forecasting accuracy,
demonstrating that the multi-head self-attention on temporal information can
enhance the model’s representation ability and improve its overall
performance.

TABLE III: Performance Comparison of Ablation Experiments

Models | Metric | Lead step
---|---|---
24(6h) | 48(12h) | 96(24h)
Pers. | MSE | $0.215$ | $0.378$ | $0.415$
MAE | $0.337$ | $0.466$ | $0.499$
FEDformer | MSE | $0.165$ | $0.225$ | $0.244$
MAE | $0.315$ | $0.373$ | $0.390$
GAT-FEDformer | MSE | $0.157$ | $0.198$ | $0.232$
MAE | $0.299$ | $0.345$ | $0.375$
GFST-WSF | MSE | $\mathbf{0.141}$ | $\mathbf{0.197}$ | $\mathbf{0.226}$
MAE | $\mathbf{0.282}$ | $\mathbf{0.344}$ | $\mathbf{0.374}$

Based on the results of the ablation experiment, which are shown in Fig.5 and
TABLE III, it is evident that the GFST-WSF model exhibits the highest forecast
performance. Specifically, the results indicate that each module plays a
significant role in the model’s overall performance, and that the removal of
any one module significantly impairs the forecast capabilities of the model.
Thus, the GFST-WSF model stands out as a superior choice for this task, owing
to its comprehensive architecture and the high degree of synergy among its
constituent parts.

(a) Distribution of MSE for 100 instances.

(b) Distribution of MAE for 100 instances.

Figure 6: Convergence and stability analysis for GFST-WSF on the test sets.

###  IV-E Convergence and Stability of GFST-WSF

In order to verify the convergence and stability of GFST-WSF, we re-
implemented the training of 100 instances of GFST-WSF. The final distribution
of MSE and MAE on the test set are shown in Fig.6a and Fig.6b, respectively,
demonstrating that: the performance of GFST-WSF in terms of convergence and
stability is satisfactory, with a small range of fluctuation in both MSE and
MAE, as shown in TABLE IV. Case studies and convergence analysis reveal that
GFST-WSF has good performance and high training stability.

TABLE IV: Convergence and stability analysis

Lead step | Metric | Avg | Max | Min | Std.
---|---|---|---|---|---
6-hour | MSE | $0.141$ | $0.154$ | $0.136$ | $0.297\%$
MAE | $0.282$ | $0.296$ | $0.273$ | $0.434\%$
12-hour | MSE | $0.197$ | $0.209$ | $0.188$ | $0.438\%$
MAE | $0.344$ | $0.363$ | $0.334$ | $0.522\%$
24-hour | MSE | $0.226$ | $0.243$ | $0.217$ | $0.474\%$
MAE | $0.374$ | $0.392$ | $0.365$ | $0.475\%$

###  IV-F Summary of Case Studies

In summary of the case studies conducted on the test set, the following
phenomena can be observed:

1) LSTM performs poorly compared to deep Transformer-based forecasting models.
As the forecasting length increases, the performance of LSTM decreases, while
Transformer-based models perform better.

2) FEDformer, which is the latest variant of Transformer, exhibits
increasingly significant performance improvements as the forecasting length
increases.

3) the introduction of graph neural networks (GNNs) enhances the model’s
performance compared to methods that only utilize single wind farm data. This
is because GNNs can extract spatial features and improve the model’s ability
to handle uncertainty in wind speed data.

4) GFST-WSF outperforms all benchmark experiments because it not only extracts
spatio-temporal features of wind speed data but also performs representation
learning in both time and frequency domains.

##  V Conclusion

In this paper, we proposes a novel GFST-WSF model for short-term wind speed
forecasting based on spatio-temporal information. It is the first work to
apply Transformer to wind speed forecasting and design a dynamic complex
adjacency matrix for GAT. The set of wind farms is modeled as a graph, where
nodes with high correlation in wind speed and direction are connected by an
edge. The connectivity of nodes in the graph is represented by a complex
adjacency matrix, which clearly characterizes the wind speed correlation and
time lag between different nodes. GAT can better capture the spatial features
of wind speeds based on this matrix. In addition, the variant of Transformer,
FEDformer, with its representation learning ability for time series, enables
GFST-WSF to better forecast wind speed sequences. The spatio-temporal features
obtained by the model can handle noise and uncertainty in wind speed data. The
proposed model is verified by comparison with various benchmark testing
methods.

Case studies on multi-wind farm data demonstrated the superiority of GFST-WSF
in wind speed forecasting. On the testing dataset, compared with state-of-the-
art deep learning models, the GFST-WSF achieved a 14%, 12%, and 8% decrease in
MSE when forecasting wind speeds for 6, 12, and 24 hours, respectively. These
results not only demonstrate the successful application of GFST-WSF in
capturing spatial features, but also suggest that deep learning models with
deep architectures can better capture complex wind speed variations.

## References

  * [1] G. W. E. Council, “Global wind report 2022. 2022,” _URL https://gwec. net/global-wind-report-2022/,(Accessed: 2022-08-17)_.
  * [2] X. Wang, P. Guo, and X. Huang, “A review of wind power forecasting models,” _Energy procedia_ , vol. 12, pp. 770–778, 2011.
  * [3] H.-z. Wang, G.-q. Li, G.-b. Wang, J.-c. Peng, H. Jiang, and Y.-t. Liu, “Deep learning based ensemble approach for probabilistic wind power forecasting,” _Applied energy_ , vol. 188, pp. 56–70, 2017.
  * [4] Y. Wang, R. Zou, F. Liu, L. Zhang, and Q. Liu, “A review of wind speed and wind power forecasting with deep neural networks,” _Applied Energy_ , vol. 304, p. 117766, 2021.
  * [5] Y. Zhang, J. Wang, and X. Wang, “Review on probabilistic forecasting of wind power generation,” _Renewable and Sustainable Energy Reviews_ , vol. 32, pp. 255–270, 2014.
  * [6] D. Nohara, M. Ohba, T. Watanabe, and S. Kadokura, “Probabilistic wind power prediction based on ensemble weather forecasting,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 12 151–12 156, 2020.
  * [7] J. J. Traiteur, D. J. Callicutt, M. Smith, and S. B. Roy, “A short-term ensemble wind speed forecasting system for wind power applications,” _Journal of Applied Meteorology and Climatology_ , vol. 51, no. 10, pp. 1763–1774, 2012.
  * [8] J. Wang, Y. Song, F. Liu, and R. Hou, “Analysis and application of forecasting models in wind power integration: A review of multi-step-ahead wind speed forecasting models,” _Renewable and Sustainable Energy Reviews_ , vol. 60, pp. 960–981, 2016.
  * [9] S. Rajagopalan and S. Santoso, “Wind power forecasting and error analysis using the autoregressive moving average modeling,” in _2009 IEEE power & energy society general meeting_. IEEE, 2009, pp. 1–6.
  * [10] R. G. Kavasseri and K. Seetharaman, “Day-ahead wind speed forecasting using f-arima models,” _Renewable Energy_ , vol. 34, no. 5, pp. 1388–1393, 2009\.
  * [11] K. Chen and J. Yu, “Short-term wind speed prediction using an unscented kalman filter based state-space support vector regression approach,” _Applied energy_ , vol. 113, pp. 690–705, 2014.
  * [12] X. Luo, J. Sun, L. Wang, W. Wang, W. Zhao, J. Wu, J.-H. Wang, and Z. Zhang, “Short-term wind speed forecasting via stacked extreme learning machine with generalized correntropy,” _IEEE Transactions on Industrial Informatics_ , vol. 14, no. 11, pp. 4963–4971, 2018.
  * [13] G. Ke, Q. Meng, T. Finley, T. Wang, W. Chen, W. Ma, Q. Ye, and T.-Y. Liu, “Lightgbm: A highly efficient gradient boosting decision tree,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [14] E. Cadenas and W. Rivera, “Wind speed forecasting in three different regions of mexico, using a hybrid arima–ann model,” _Renewable Energy_ , vol. 35, no. 12, pp. 2732–2738, 2010.
  * [15] Y.-Y. Hong and T. R. A. Satriani, “Day-ahead spatiotemporal wind speed forecasting using robust design-based deep learning neural network,” _Energy_ , vol. 209, p. 118441, 2020.
  * [16] Q. Cao, B. T. Ewing, and M. A. Thompson, “Forecasting wind speed with recurrent neural networks,” _European Journal of Operational Research_ , vol. 221, no. 1, pp. 148–154, 2012.
  * [17] Y.-L. Hu and L. Chen, “A nonlinear hybrid wind speed forecasting model using lstm network, hysteretic elm and differential evolution algorithm,” _Energy conversion and management_ , vol. 173, pp. 123–142, 2018.
  * [18] T. Hu, W. Wu, Q. Guo, H. Sun, L. Shi, and X. Shen, “Very short-term spatial and temporal wind power forecasting: A deep learning approach,” _CSEE Journal of Power and Energy Systems_ , vol. 6, no. 2, pp. 434–443, 2019.
  * [19] M. He, L. Yang, J. Zhang, and V. Vittal, “A spatio-temporal analysis approach for short-term forecast of wind farm generation,” _IEEE Transactions on Power Systems_ , vol. 29, no. 4, pp. 1611–1622, 2014.
  * [20] Q. Zhu, J. Chen, L. Zhu, X. Duan, and Y. Liu, “Wind speed prediction with spatio–temporal correlation: A deep learning approach,” _Energies_ , vol. 11, no. 4, p. 705, 2018.
  * [21] Q. Zhu, J. Chen, D. Shi, L. Zhu, X. Bai, X. Duan, and Y. Liu, “Learning temporal and spatial correlations jointly: A unified framework for wind speed prediction,” _IEEE Transactions on Sustainable Energy_ , vol. 11, no. 1, pp. 509–523, 2019.
  * [22] R. Chen, J. Liu, F. Wang, H. Ren, and Z. Zhen, “Graph neural network-based wind farm cluster speed prediction,” in _2020 IEEE 3rd Student Conference on Electrical Machines and Systems (SCEMS)_. IEEE, 2020, pp. 982–987.
  * [23] P. Velickovic, G. Cucurull, A. Casanova, A. Romero, P. Lio, Y. Bengio _et al._ , “Graph attention networks,” _stat_ , vol. 1050, no. 20, pp. 10–48 550, 2017.
  * [24] Q. Wen, T. Zhou, C. Zhang, W. Chen, Z. Ma, J. Yan, and L. Sun, “Transformers in time series: A survey,” _arXiv preprint arXiv:2202.07125_ , 2022.
  * [25] T. Zhou, Z. Ma, Q. Wen, X. Wang, L. Sun, and R. Jin, “Fedformer: Frequency enhanced decomposed transformer for long-term series forecasting,” in _International Conference on Machine Learning_. PMLR, 2022, pp. 27 268–27 286.
  * [26] C. C. Robusto, “The cosine-haversine formula,” _The American Mathematical Monthly_ , vol. 64, no. 1, pp. 38–40, 1957.
  * [27] J. Zhou, G. Cui, S. Hu, Z. Zhang, C. Yang, Z. Liu, L. Wang, C. Li, and M. Sun, “Graph neural networks: A review of methods and applications,” _AI open_ , vol. 1, pp. 57–81, 2020.
  * [28] C. Draxl, A. Clifton, B.-M. Hodge, and J. McCaa, “The wind integration national dataset (wind) toolkit,” _Applied Energy_ , vol. 151, pp. 355–366, 2015.


# Topology with memory in nonlinear driven-dissipative photonic lattices

Subhaskar Mandal<EMAIL_ADDRESS>Gui-Geng Liu  Baile Zhang
<EMAIL_ADDRESS>Division of Physics and Applied Physics, School of Physical
and Mathematical Sciences, Nanyang Technological University, Singapore 637371,
Singapore

###### Abstract

We consider a photonic lattice of nonlinear lossy resonators subjected to a
coherent drive, where the system remembers its topological phase. Initially,
the system is topologically trivial. After the application of an additional
coherent pulse, the intensity is increased, which modifies the couplings in
the system and then induce a topological phase transition. However, when the
effect of the pulse dies out, the system does not go back to the trivial
phase. Instead, it remembers the topological phase and maintains its topology
acquired during the pulse application. The pulse can be used as a switch to
trigger amplification of the topological modes. We further show that the
amplification takes place at a different frequency as well as at a different
position from those of the pulse, indicating frequency conversion and
intensity transfer. Our work can be useful in triggering the different
functionalities of active topological photonic devices.

## Introduction

The intriguing properties of topological photonics have enabled widespread
applications in modern optical devices, such as robust signal transport Wang
_et al._ (2009); Poo  _et al._ (2011); Shalaev  _et al._ (2019); He  _et al._
(2019); Yang  _et al._ (2020), optical delay line Hafezi  _et al._ (2011),
quantum interface Barik  _et al._ (2018), quantum light source Mittal  _et
al._ (2018), robust splitters Cheng  _et al._ (2016), topological lasers St-
Jean  _et al._ (2017); Bahari  _et al._ (2017); Bandres  _et al._ (2018); Zeng
_et al._ (2020); Smirnova  _et al._ (2020a), etc. Topological photonics is
also promising for optical information processing technologies. For example,
valley photonic crystals are identified as an excellent candidate for robust
information transfer in next-generation devices Shalaev  _et al._ (2019); He
_et al._ (2019); Yang  _et al._ (2020). Similar to transferring the
information, the ability to store them in memory, is an equally important task
in information processing. However, optical memories along with topological
protection have not been explored till now.

Nonlinearity is at the core of memory devices. Interplay between the
nonlinearity and the topology has made way for many novel effects such as
topological solitons Leykam and Chong (2016); Kartashov and Skryabin (2016);
Mukherjee and Rechtsman (2020, 2021); Pernet  _et al._ (2022), high harmonic
generation Kruk  _et al._ (2019); Wang  _et al._ (2019), topological phase
transitions Bardyn  _et al._ (2016); Hadad  _et al._ (2018); Zhou  _et al._
(2017); Mandal  _et al._ (2019); Zangeneh-Nejad and Fleury (2019); Maczewsky
_et al._ (2020), and others Lumer  _et al._ (2013); Kartashov and Skryabin
(2017); Banerjee  _et al._ (2020); Xia  _et al._ (2021); Kirsch  _et al._
(2021); Ezawa (2021) (see Ref. Smirnova  _et al._ (2020b) for a comprehensive
review). However, none of the previous works can show the memory feature: once
the key ingredient, that induces the functionalities, is removed from the
scheme, the systems can no longer continue to exhibit such effects.

In this work, we introduce for the first time a topological phase with memory
in a lattice of lossy resonators having local onsite Kerr-nonlinearity, where
the system remembers its topological phase. The lossy nature of the system
leads to a steady-state in the presence of a coherent drive $F$. However, due
to the nonlinearity, our system subjected to a properly designed $F$ shows not
only one but two steady states: low and high-intensity states. Our system is
topologically trivial and after $F$ is introduced, it attains the low-
intensity steady-state. The introduction of an additional coherent pulse
increases the intensity of the system. At higher intensities, the nonlinear
interaction modifies the couplings and the otherwise-trivial system becomes
topological. However, at longer times when the effect of the pulse dies, the
system does not go back to its previous trivial phase. Instead, it remembers
the topological phase and maintains its topology acquired during the pulse
application. As an application of this effect, we show a unique amplification
phenomenon, where the amplification is triggered by a pulse.

Figure 1: (a) Schematic of a coherently driven nonlinear resonator. (b)
Analytically and numerically calculated bistable curves in green and red,
respectively. Parameters: $\Delta=-3$, $\Gamma=1$.

We start by considering a nonlinear optical resonator subjected to a coherent
drive $F$ (see Figure 1a), which is represented by the following nonlinear
schrödinger equation (NLSE)

$\displaystyle i\frac{\partial\psi}{\partial
t}=\left(\omega_{0}-i\Gamma\right)\psi+|\psi|^{2}\psi+F\exp(-i\omega_{p}t).$
(1)

Here $\omega_{0}$ is the onsite potential and $\Gamma$ is the linear decay.
The next term represents the defocusing Kerr nonlinearity, where the nonlinear
coefficient is set to 1. $F$ is a coherent drive having frequency
$\omega_{p}$. For the steady-state $\psi_{s}$, where
$\partial\psi_{s}/\partial t=0$, one can obtain

$\displaystyle\left|F\right|^{2}=\left[\left(\Delta+\left|\psi_{s}\right|^{2}\right)^{2}+\Gamma^{2}\right]\left|\psi_{s}\right|^{2},$
(2)

where $\Delta=\left(\omega_{0}-\omega_{p}\right)$. From eq. (2), it is easy to
find that within the gray region for each $|F|^{2}$ three possible
$|\psi_{s}|^{2}$ exist (see Figure 1b). However, in practice, the middle
branch is not stable. This can be confirmed by numerically solving eq. (1),
but letting $F$ vary very slowly in time such that at each time step steady
state can be reached (more details on the numerical calculation of bistability
can be found in the supporting information). The red dashed curve obtained
numerically follows the analytical green curve, however, the middle branch
does not appear. Consequently, the system shows bistability by allowing both
the low and high-intensity stable states for a fixed value of $F$ within the
gray region. The absence of the middle branch of the Bistability curve can be
explained using the stability analysis based on the first Lyapunov Criterion
Yu  _et al._ (2019) (see supporting information).

An important characteristic of bistability is their ability to mimic the
memory: the state of the system is not only determined by the current
parameters (such as $F$) but also by its previous state. For example, let us
consider a system that is initially in the low-intensity state “$A$” as shown
in Figure 1b. An additional coherent drive $\Delta F$ is added such that the
system moves to a high-intensity state “$B$”. Now if $\Delta F$ is removed,
the system does not return back to its original state “$A$”, instead, it
chooses the high-intensity state “$C$”. While determining the final state, the
system memorizes the information (high-intensity) about the intermediate state
“$B$” and in the case where “$B$” is a low-intensity state, the system would
return to the original state “$A$” upon removing $\Delta F$.

## Model

We arrange the nonlinear resonators in a 2D lattice. In order to be close with
experiments, we model the dynamics of the system using the NLSE in the
continuum limit (where the space is taken as continuous). Without the loss of
generality, we work with the dimensionless NLSE, which is expressed as

$\displaystyle i\frac{\partial\psi(x,y)}{\partial t}=$
$\displaystyle\left[-\nabla^{2}+V(x,y)-i\Gamma\right]\psi(x,y)$
$\displaystyle+\left|\psi(x,y)\right|^{2}\psi(x,y)+F(x,y)e^{-i\omega_{p}t}$
$\displaystyle+F_{p}(x,y)\exp\left[-\frac{(t-t_{0})^{2}}{2\tau^{2}}\right]e^{-i\omega_{p}t}.$
(3)

Here $\nabla^{2}\equiv\left({\partial^{2}}/{\partial
x^{2}}+{\partial^{2}}/{\partial y^{2}}\right)$ is the transverse Laplacian
operator, $V$ is the external potential profile corresponding to the
resonators, $\Gamma$ is the linear decay, $F$ is a position-dependent coherent
pump having frequency $\omega_{p}$, and $F_{p}$ is a coherent pulse having
duration $\tau$ centered at time $t_{0}$. Next, we consider circular
resonators having diameter $d_{m}$, which we call main resonators. Two main
resonators separated by $L$ are coupled via an auxiliary larger resonator
having diameter $d_{a}$, where $d_{a}>d_{m}$ (see Figure 2a). The potential is
taken as $V=0$ inside and $V=V_{0}>0$ outside the resonators.

Figure 2: (a) Schematic of two main resonators connected by an auxiliary
resonator. (b,c,d) The lowest three modes of the system shown in (a). (e) The
ratio of the coupling in the nonlinear regime $J_{NL}$ to the linear coupling
$J$ between two main resonators as a function of the intensity of the
auxiliary resonator. Parameters: $V_{0}=236$, $L=2.13$, $d_{m}=1$,
$d_{a}=1.13$. Figure 3: (a) Spatial profile of the coherent pump. (b,c)
Steady-states of the system before and after the coherent pulse $F_{p}$,
respectively. The green circle in (b) represents the width of $F_{p}$. (d) The
total intensity of the system as a function of time showing the bistable
behavior. The two arrows indicate the times at which the states in (b,c) are
plotted. Parameters: $\Gamma=0.13$, $P_{0}=\sqrt{0.5}$, $\omega_{p}=14.74$,
$\sigma=0.3$, $\tau=23.9$, $t_{0}=77$. $F_{p}$ is a gaussian pulse having
strength $20P_{0}$ and width $7.5\sigma$. All other paramters are kept the
same as those in Figure 2.

To capture the role of nonlinearity, we first consider two main resonators
connected by an auxiliary resonator as shown in Figure 2a. The ground state
wave function $\psi_{0}$ is mainly localized at the auxiliary resonator (see
Figure 2b), whereas the first ($\psi_{1}$) and second ($\psi_{2}$) excited
states are localized at the main resonators (see Figure 2c,d). The coupling
strength $J$ between the main resonators can be estimated from the difference
in the eigenvalues of the symmetric ($E_{2}$) and asymmetric ($E_{1}$)
eigenstates, where $J=\left(E_{2}-E_{1}\right)/2$. The important feature that
plays a key role and signifies the nonlinear effect in this work is the
ability to control the coupling between the main resonators by changing the
intensity of the auxiliary resonator. This is captured by choosing an
effective potential $V_{\text{eff}}=V+g\left|\psi_{0}\right|^{2}$, where $g$
corresponds to the peak value of the ground state intensity, and obtaining the
coupling in the nonlinear regime $J_{NL}$ in a self-consistent way. Figure 2e
shows the enhancement of $J_{NL}$ compared to $J$ as a function of $g$.
However, such an enhancement of $J_{NL}$ is limited. Once $g$ becomes larger
than the difference between the fundamental frequencies of the main and the
auxiliary resonators, $J_{NL}$ would start to decrease and for
$g\rightarrow\infty$, $J_{NL}\rightarrow 0$.

Figure 4: (a,b,c) Real eigenfrequencies for different cases. In (a) the
auxiliary resonator band is shown in gray and the main resonator band is shown
in black. Red and green dots in (c) correspond to the corner and edge modes,
respectively. The blue line in (a) represents $\omega_{p}$, with respect to
which (b,c) are rescaled. (d,e) Spatial profiles of an edge and a corner mode,
respectively. All the parameters are kept the same as those in Figure 3.

## Topological memory

Now that we have all the ingredients, we proceed to study the topological
phase in a 2D square lattice formed by the above-mentioned resonators, where
between any two main resonators there is an auxiliary resonator. Recalling
that the intensity of the auxiliary resonators enhances the coupling between
the main resonators, we choose the coherent pump profile in such a way that
the inter-cells are coupled strongly similar to the 2D Su-Schrieffer–Heeger
(SSH) model Xie  _et al._ (2018). The spatial profile of the coherent pump is
expressed as

$\displaystyle
F(x,y)=P_{0}\sum_{X_{n},Y_{n}}\exp\left[-\frac{(x-X_{n})^{2}+(y-Y_{n})^{2}}{2\sigma^{2}}\right],$
(4)

where $P_{0}$ is the strength of the pump and $\\{X_{n},Y_{n}\\}$ are the
coordinates of the center of the pumped auxiliary resonators as shown in
Figure 3a. Such pump profiles are readily achieved in practice using the
spatial light modulators Moreno  _et al._ (2012); Töpfer  _et al._ (2021). We
choose the value of $P_{0}$ such that the auxiliary resonators subjected to
the pump are bistable. To show the bistable behaviour of the whole system, we
solve eq. (Model) without the pulse ($F_{p}=0$) and take zeros as the initial
condition. While plotting the spatial profiles, we keep the background
potential to distinguish the intensities between the main and auxiliary
resonators. Figure 3b shows the steady-state of the system before the
application of the pulse, where intensity is mainly localized at the pumped
auxiliary resonators. The coupling between the resonators results in slightly
non-identical bistability curves of the resonators placed at different
positions (see supporting information). Due to this, the intensity among the
pumped resonators varies a little, but they remain in the low-intensity state,
where $|\psi|^{2}$ is negligible.

Next, we apply a gaussian-shaped coherent pulse $F_{p}$ centered at the bulk
as shown by the green circle in Figure 3b. The addition and removal of the
additional pump $\Delta F$ in Figure 1b is played by the pulse $F_{p}$ here.
In Figure 3c the steady-state of the system after the application of the pulse
is shown. The system indeed remembers the high intensity created by the pulse
and once the effect of the pulse dies out the system chooses to stay at the
high-intensity state. Compared to the low-intensity state in Figure 3b, a much
larger intensity outside the pumped auxiliary resonators exists, which
signifies the enhancement of the coupling due to significant $|\psi|^{2}$. In
Figure 3d the total intensity of the system, $I(t)=\int\psi(x,y,t)~{}dxdy$,
where the integration is over the whole system, is shown as a function of
time, which shows the bistable behaviour of the system. The full dynamics of
the system is shown in movie1.

We have performed all the calculations corresponding to the 2D lattice on a
$2^{9}\times 2^{9}$ grid. The Laplacian is taken into account through the FFT
(Fast Fourier Transform) spectral method. It should be noted that the finite
difference (FD) method can also be implemented to express the Laplacian.
However, FD requires larger computational memory and is more time-consuming
compared to the FFT method. The time dynamics is performed using Matlab’s ODE
solvers, which relies on well-established numerical techniques, such as the
explicit Runge-Kutta (4,5) formula, the Dormand-Prince pair Dormand and Prince
(1980).

Having established the memory effect in the 2D lattice, here we show the
topology associated with it. In Figure 4a the eigenfrequencies of the linear
system are shown, which can be found by putting the nonlinear and pumping
terms to zero in eq. (Model) and diagonalizing its corresponding Hamiltonian.
The lower band (shown in gray) has the main contributions from to the
auxiliary resonators, whereas the upper band has the main contributions from
the main resonators. For the rest of the work, we shall focus on the main
resonator band, which is topologically trivial and gapless in the linear
regime. To include the nonlinear effect, we study the Bogoliubov fluctuations
on top of the steady-state Pitaevskii and Stringari. (2003):

$\displaystyle\psi(x,y)=\psi_{s}(x,y)+u_{n}(x,y)e^{-i\omega_{n}t}+v_{n}^{*}(x,y)e^{i\omega_{n}^{*}t}.$
(5)

Here $\psi_{s}$ represent the low and high-intensity steady states shown in
Figure 3b,c, respectively. $u_{n}$ and $v_{n}$ represent the fluctuations
having frequency $\omega_{n}$. Substituting eq. (5) into eq. (Model) and by
ignoring the higher-order terms in $u_{n}$ and $v_{n}$, we obtain the
following eigenvalue equation

$\begin{bmatrix}H_{0}+2\left|\psi_{s}\right|^{2}&\psi_{s}^{2}\\\
-{\psi_{s}^{2}}^{*}&-H_{0}^{*}-2\left|\psi_{s}\right|^{2}\end{bmatrix}\begin{bmatrix}u_{n}\\\
v_{n}\end{bmatrix}=\omega_{n}\begin{bmatrix}u_{n}\\\ v_{n}\end{bmatrix},$ (6)

where $H_{0}=\left[-\nabla^{2}+V(x,y)-i\Gamma-\omega_{p}\right]$, which is
rescaled with respect to the pump frequency $\omega_{p}$.

The fluctuation Hamiltonian has particle-hole symmetry, which makes the
eigenfrequencies appear in pairs ($\omega_{n},~{}-\omega_{n}^{*}$). The
Hilbert space of the fluctuation Hamiltonian is double in size compared to the
linear one. Consequently, for better visualization we show the
eigenfrequencies near the main resonator band and for Re $\omega_{n}>0$.
Figure 4b shows the eigenfrequencies of the fluctuations before the pulse is
applied. Similar to the linear case, the band is gapless. This is
understandable as for the low-intensity state, $|\psi|^{2}$ is not significant
enough to induce the topological transition. After the application of the
pulse positioned at the bulk, the system switches to the high-intensity state.
In this case, the nonlinear effect becomes significant and the coupling
between the main resonators connected by a pumped auxiliary resonator
increases. The system goes through a topological phase transition from trivial
to second-order topological phase, where a bulk band gap opens and four
topological corner modes appear. An effective tight-binding model based on the
strong coupling induced by the pump can reproduce the topological corner modes
(see supporting information).

In Figure 4c the eigenfrequencies of the fluctuation after the application of
the pulse are shown, where the topological corner modes are marked in red and
the edge modes are shown in green. Figure 4d,e show the spatial profiles of
one of the topological edge modes ($n=146$) and topological corner modes
($n=142$), respectively. In experiment, the topological modes will be hidden
in the high-intensity steady state. However, they can be probed using a weak
additional coherent pump followed by the frequency filtration to subtract the
steady-state.

Higher order topological phases have been an intense area of research
Benalcazar  _et al._ (2017a); Liu  _et al._ (2019); Xie  _et al._ (2018). We
note that the effect of the on-site Kerr-nonlinearity on these system have
been studied recently Banerjee  _et al._ (2020); Ezawa (2021); Hu  _et al._
(2021); Kirsch  _et al._ (2021). However, the previous works consider the
effect of nonlinearity on the already existing linear topological band
structure. The presented result is the first example, where onsite Kerr-
nonlinearity alone induces higher-order topological phase transition. Although
here we have focused on the topological corner modes, it would be interesting
to investigate further whether different truncations can lead to different
types of edge states, such as the ones shown in Ref. Cheng  _et al._ (2022).

## Calculation of the topological invariant

Here we calculate the bulk polarization, which characterizes the topological
corner modes Benalcazar  _et al._ (2017b) using a bi-orthogonal Wilson loop
Luo and Zhang (2019). We define a Wilson line along the $y$-direction as

$\displaystyle\Xi^{m,n}_{y}(k_{x},k_{y})$
$\displaystyle=\langle\Psi^{L,m}_{k_{x},k_{y}+\Delta_{y}}|\Psi^{R,n}_{k_{x},k_{y}}\rangle$
$\displaystyle=\int_{\text{unit
cell}}\Psi^{L,m*}_{k_{x},k_{y}+\Delta_{y}}(x,y)\Psi^{R,n}_{k_{x},k_{y}}(x,y)~{}dxdy,$
(7)

where $\Psi^{n}_{k_{x},k_{y}}$ represents the Bloch eigenvectors of the
fluctuations (see supporting information for details on the Bloch eigenvector
calculation), $L$ and $R$ correspond to the left and right Bloch eigenstates,
respectively. $\Delta_{y}=2\pi/N_{y}$, where $N_{y}$ is the total number of
points used along the $y$ direction of the Brillouin zone for calculation.
Taking the periodicity as one, the Wilson loop along the $y$ direction is
defined as

$\displaystyle W^{y}(k_{x},k_{y})=$
$\displaystyle\Xi^{m,n}_{y}(k_{x},k_{y}+2\pi).....\Xi^{m,n}_{y}(k_{x},k_{y}+2\Delta_{y})$
$\displaystyle\Xi^{m,n}_{y}(k_{x},k_{y}+\Delta_{y})~{}\Xi^{m,n}_{y}(k_{x},k_{y}).$
(8)

The Wannier Hamiltonian is given by

$\displaystyle
H_{W}^{y}(k_{x},k_{y})=-\frac{i}{2\pi}\log\left[W^{y}(k_{x},k_{y})\right].$
(9)

The eigenvalues of $H_{W}^{y}$ form the Wannier bands, which are shown in
Figure 5. The bulk polarization $P=(P_{x},P_{y})$ is the same as the Wannier
center. For a topologically trivial system $P=(0,0)$. The $x$ component of the
bulk polarization $P_{x}$ is given by summing all the eigenvalues $\nu_{y}$
corresponding to all the momentum $k_{x}$ Benalcazar  _et al._ (2017b). It can
be seen from Figure 5 that $P_{x}\approx 0.5$. Since the system has $C_{4}$
rotational symmetry, one can obtain a similar result by choosing the Wilson
loop along the $x$-direction, such that the $y$ component of the bulk
polarization $P_{y}\approx 0.5$ and the total polarization becomes
$P\approx(0.5,0.5)$ making the system topologically nontrivial.

Figure 5: Plot of the Wannier bands obtained by diagonalizing eq. (9). Figure
6: (a,b,c) Complex eigenfrequencies for different cases. Steady states from
Figure 3b,c are used for obtaining (b,c). (d,e) Steady states of the system
with gain at the corners before and after the pulse is applied, respectively.
The green circle in (d) shows the position and width of the pulse. Parameters:
Peak value of the gain $G_{0}=3\Gamma$. All other parameters are kept the same
as those in the Figure 3.

## Amplification of the corner modes

As an application, we use our proposed scheme to control the functionality of
an active topological photonic device. We introduce gain at the main four
corner resonators of the 2D lattice. The gain is modeled by adding a term
$+iG(x,y)\psi(x,y)$ at the right-hand side of eq. (Model), where $G(x,y)$ is
composed of four Gaussians centered at the four corners having width $\sigma$
and peak value $G_{0}$. In this case, $H_{0}$ in eq. (6) is updated to
$H_{0}\rightarrow H_{0}+iG(x,y)$. To signify the role of nonlinearity, it is
important that the gain alone can not induce lasing in the linear regime.
Consequently, we obtain the complex eigenfrequencies for the linear system.
The system stays below the lasing threshold Im $\omega_{n}<0$ (see Figure 6a).
Next, we take the steady states corresponding to Figure 3b,c and obtain the
same plot for the fluctuations before and after the pulse is applied. Similar
to the linear case, the modes have Im $\omega_{n}<0$ before the application of
the pulse (see Figure 6b). However, after the pulse is applied topological
corner modes appear and due to the significant overlap with the gain, only
they have Im $\omega_{n}>0$, while all other modes remains at Im
$\omega_{n}<0$ (see Figure 6c). This has a significant effect on the steady-
states. As predicted from the complex eigenfrequencies, the gain at the four
corners does not alter the steady-state before the pulse, which is the same as
the one obtained without the gain in Figure 3b. However, after the pulse is
applied, a large intensity at the corners is observed along with the high-
intensity steady state at the bulk (see movie2).

To confirm that the intensity at the corners correspond to the topological
corner modes, we further obtain their spatial profile from the time-dynamics,
which does not rely on the linear Bogoliubov theory (which does not include
the higher-order terms in $u$ and $v$). We store the solutions at the
intermediate time steps corresponding to the steady states in Figure 6d,e and
then Fourier transform them along the time axis to move to the frequency
dimension. At this stage, we can plot the intensity as a function of
frequency. In order to obtain the spatial profile corresponding to a
particular frequency, we first multiply with a gaussian to filter the desired
frequency and plot the intensity. In Figure 7a, the above-mentioned steps are
given. In Figure 7b,c,d,e the frequency-dependent intensities are shown for
different cases. Without the gain, we obtain a peak at the pumped frequency
$\omega_{p}$ for both before and after the application of the pulse. It is
important to note that no other peak is observed. Now we move on to the case
where we introduce gain at the four corners. In this case also, before the
application of the pulse a single peak at $\omega_{p}$ is obtained. However,
after the pulse is applied a new peak along with the pumped peak is obtained
(see Figure 7e). It proves that the amplification at a different frequency has
taken place after the application of the pulse. To prove that the amplified
modes are the corner modes, we plot the spatial profile of the wave function
around the amplified frequency by filtering out the desired frequency using
the recipe mentioned in Figure 7a. From the spatial profile shown in Figure
7f, it is clear that the amplified modes are indeed the corner modes. Due to
the presence of the $\alpha|\psi|^{2}$ term the corner modes are shifted in
frequency from their linear Bogoliubov spectrum, which does not include the
higher-order terms in $u$ and $v$.

Figure 7: (a) Steps to obtain spatial profiles of $\psi$ at a particular
frequency by solving the NLSE. (b,c,d,e) Frequency-dependent intensity. The
second peak in (e) corresponds to the amplification. (f) The spatial profile
corresponds to the amplification peak in (e), which has the profile of corner
modes. In each case we normalize the intensity with respect to the intensity
at $\omega_{p}$.

There are a few unique feature associated with the amplification process
presented above. Firstly, the pulse is positioned away from the corners at the
bulk. However, the amplification takes place at the corners of the system,
which do not overlaps with the pulse. This can be interpreted as the intensity
transfer from the bulk to the corner. However, we stress that such a process
is completely different from the signal transfer typically associated with the
1D edge modes of first-order topological insulators. Secondly, the amplified
corner states have frequency different from that of the coherent drive, which
makes this system suitable for frequency conversion. This is similar to the
recent works on the high harmonic generation using topological systems Kruk
_et al._ (2019); Wang  _et al._ (2019); Zograf  _et al._ (2022); Smirnova  _et
al._ (2019), although here our underlying linear system is topologically
trivial. Lastly, the amplification, which is triggered by a pulse, is
different from all the previous cases. For example, in lasers if the pumping
term (or any other ingredient), which induces amplification (lasing) is
removed, understandably the amplification would stop. Due to the memory
effect, the functionality in our scheme remains, although the pulse
disappears.

We note that a different type of topological memory is proposed recently in
dynamic PT-symmetric optical resonators having saturable nonlinearity Choi
_et al._ (2022); Yu  _et al._ (2021). Unlike the present case, a lattice is
not involved there, instead, the system consists of a single diatomic or
triatomic resonator. The topology introduced in those works is different from
the topology associated with band structures in usual topological lattices.
There oscillating quenching states are protected against the defects of the
parameter space that preserve dynamical state trajectories.

## Proposal for experimental realization

The NLSE in eq. (Model) is generally used to describe the topological physics
in photonic waveguide arrays Mukherjee and Rechtsman (2020, 2021). Bistability
is realized experimentally in waveguide arrays Khomeriki and Leon (2005),
whereas the Bogoliubov fluctuations in eq. (6) can be arranged using the
parametric down-conversion Solntsev  _et al._ (2014); Kruse  _et al._ (2013);
Peano  _et al._ (2016); Yang  _et al._ (2021); Doyle  _et al._ (2022).
Alternatively, the system of exciton-polaritons, where cavity photons exhibit
Kerr-nonlinearity by coupling strongly with the quantum well excitons, is also
a promising platform for realizing our scheme. They are well known for
studying topological photonics Bardyn  _et al._ (2015); Nalitov  _et al._
(2015); Klembt  _et al._ (2018); Banerjee  _et al._ (2018); Kartashov and
Skryabin (2019); Mandal  _et al._ (2020, 2022); Harder  _et al._ (2021)..
Bistability is well established for exciton-polaritons Baas  _et al._ (2004);
Paraïso  _et al._ (2010); Banerjee and Liew (2020). Bogoliubov fluctuations
naturally arise in polariton system Utsunomiya  _et al._ (2008); Pernet  _et
al._ (2022). By choosing the proper physical units, our present parameters can
be related to the exciton-polariton lattices.

Here we provide the physical parameters based on the exciton-polariton
system.Typically, the system of exciton-polaritons is expressed using the NLSE
(also known as the nonlinear Gross-Pitaevskii equation) Carusotto and Ciuti
(2013)

$\displaystyle i\hbar\frac{\partial\phi(X,Y)}{\partial T}=$
$\displaystyle\left[-\frac{\hbar^{2}}{2m}\left(\frac{\partial^{2}}{\partial
X^{2}}+\frac{\partial^{2}}{\partial Y^{2}}\right)\right.$
$\displaystyle\left.+V_{e}(X,Y)-i\frac{\hbar}{2\lambda}+\alpha|\phi(X,Y)|^{2}\right]\phi(X,Y)$
$\displaystyle+f(X,Y)\exp(-i\Omega_{p}t).$ (10)

We restrict ourselves at the bottom of the lower polariton dispersion, where
the dispersion can be approximated to be a parabola having effective mass $m$.
$V_{e}$ is the external potential for the polaritons, which consists of
microcavity pillars Klembt  _et al._ (2018); Harder  _et al._ (2020); Dusel
_et al._ (2020). $\lambda$ is the effective lifetime of the polaritons, which
can be controlled by adjusting the quality factor of the cavity. $\alpha>0$ is
the effective polariton-polariton interaction coefficient, which is typically
repulsive in nature. $f$ is the coherent drive. Next, we perform the following
transformation in order to transform to the eq. (Model)

$\displaystyle X\rightarrow xa,~{}Y\rightarrow ya,~{}V_{e}\rightarrow
V\varepsilon_{u},~{}T\rightarrow
tt_{u},~{}\Omega_{p}\rightarrow\omega_{p}/t_{u},$
$\displaystyle\lambda\rightarrow
2t_{u}/\Gamma,~{}\phi\rightarrow\psi\sqrt{\frac{\varepsilon_{u}}{\alpha}},~{}\text{and}~{}f\rightarrow
F\frac{\varepsilon_{u}^{3/2}}{\sqrt{\alpha}},$ (11)

where $x,~{}y,~{}V,~{}\Gamma,~{}\psi,~{}\text{and}~{}F$ are the dimensionless
quantities used in eq. (Model). Here $a$ is the length unit,
$\varepsilon_{u}=\hbar^{2}/2ma^{2}$ is the energy unit, $t_{u}=2ma^{2}/\hbar$
is the time unit, $\sqrt{\varepsilon_{u}/\alpha}$ is the wave function unit,
and $\varepsilon^{3/2}_{u}/\sqrt{\alpha}$ is the pump unit. After the above
transformation in eq. (Proposal for experimental realization), eq. (Proposal
for experimental realization) can be matched exactly with eq. (Model).

We set $a=3~{}\mu$m, which makes the diameter of the main pillars 3 $\mu$m and
that for the auxiliary pillars 3.4 $\mu$m. The spacing between the main
pillars becomes 6.4 $\mu$m. Such a system of micropillars is readily achieved
in experiments Klembt  _et al._ (2018); Harder  _et al._ (2020); Dusel  _et
al._ (2020). By choosing $m=5\times 10^{-5}m_{e}$, where $m_{e}$ is the free
electron mass, the energy unit becomes $\varepsilon_{u}\approx 0.085$ meV.
Fixing the energy unit also fixes all the energy scales. The topological
bandgap becomes around 0.12 meV, the effective potential depth for the
polaritons becomes 20 meV, and the lifetime becomes 30 ps. The energy scale
can be increased by choosing smaller sized micropillars and reducing $a$,
which increases $\varepsilon_{u}$. Alternatively, by adjusting the detuning
between the exciton and photon branches, the effective mass of the polaritons
can be reduced, thereby increasing $\varepsilon_{u}$.

Another important parameter is the nonlinear interaction constant $\alpha$.
The value of $\alpha$ can be controlled by adjusting the exciton fraction of
the polaritons. However, measuring $\alpha$ exactly in experiments is
difficult and still an ongoing research. The measurable quantity in
experiments is the blueshift $\alpha|\phi|^{2}$. For the above-mentioned
parameters the blueshift for the low-intensity state is around 0.03 meV and
that for the high-intensity state is around 0.35 meV. Such values of blueshift
are routine observation in experiments Sun  _et al._ (2017); Estrecho  _et
al._ (2019).

For the introduction of the gain at the four corners, an optical pump
positioned at the corner having much higher energy from the polariton
resonance can be used. The pump creates free electron-hole pairs, which relax
down and form the excitonic reservoir. The density of the excitons in the
excitonic reservoir act as gain to the polaritons Carusotto and Ciuti (2013).
Alternatively, it is also possible to arrange gain using electrical pumping
Schneider  _et al._ (2013); Klaas  _et al._ (2017); Gagel  _et al._ (2022).

## Conclusion

To conclude, we have presented a new concept, where the topological phase can
be induced as a memory. In particular, we show that a nonlinear system of
photonic lossy resonators goes through a topological phase transition under
the application of a coherent pulse. The system continues to maintain its
topology, although the effect of the pulse disappears. The topological modes
can show fair robustness against realistic disorder (see supporting
information). This scheme is independent of the dimension and similar effects
can be found in 1D lattices (see supporting information). Our scheme is the
first example, where onsite Kerr-nonlinearity induces higher-order topology
and can be useful in triggering different functionalities of active
topological photonic devices.

### acknowledgement

S. M. thanks Dr. Rimi Banerjee for helpful discussions.

### Funding Sources

This research is supported by the Singapore National Research Foundation
Competitive Research Program under Grant No. NRF-CRP23-2019-0007, and the
Singapore Ministry of Education Academic Research Fund Tier 3 Grant
MOE2016-T3-1-006.

## Supporting Information Available

This material is available free of charge via the internet at
http://pubs.acs.org. The supporting information contains some additional
calculations and two movies. The additional calculations include the
discussion on bistability, effective tight binding model, calculation of the
topological invariant, robustness of the scheme, and 1D topological memory.
Movie 1 shows the memory effect without the gain at the corner. Movie 2 shows
the triggering of the amplification of the corner modes with a pulse.

## References

  * Wang  _et al._ (2009) Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljačić, “Observation of unidirectional backscattering-immune topological electromagnetic states,” Nature 461, 772–775 (2009).
  * Poo  _et al._ (2011) Y. Poo, R. Wu, Z. Lin, Y. Yang, and C. T. Chan, “Experimental realization of self-guiding unidirectional electromagnetic edge states,” Phys. Rev. Lett. 106, 093903 (2011).
  * Shalaev  _et al._ (2019) M. I. Shalaev, W. Walasik, A. Tsukernik, Y. Xu, and N. M. Litchinitser, “Robust topologically protected transport in photonic crystals at telecommunication wavelengths,” Nat. Nanotechnol. 14, 31–34 (2019).
  * He  _et al._ (2019) X. T. He, E. T. Liang, J. J. Yuan, H. Y. Qiu, X. D. Chen, F. L. Zhao, and J. W. Dong, “A silicon-on-insulator slab for topological valley transport,” Nat. Commun. 10, 872 (2019).
  * Yang  _et al._ (2020) Y. Yang, Y. Yamagami, X. Yu, P. Pitchappa, J. Webber, B. Zhang, M. Fujita, T. Nagatsuma, and R. Singh, “Terahertz topological photonics for on-chip communication,” Nat. Photonics 14, 446–451 (2020).
  * Hafezi  _et al._ (2011) M. Hafezi, E. A. Demler, M. D. Lukin, and J. M. Taylor, “Robust optical delay lines with topological protection,” Nat. Phys. 7, 907–912 (2011).
  * Barik  _et al._ (2018) S. Barik, A. Karasahin, C. Flower, T. Cai, H. Miyake, W. DeGottardi, M. Hafezi, and E. Waks, “A topological quantum optics interface,” Science 359, 666–668 (2018).
  * Mittal  _et al._ (2018) S. Mittal, E. A. Goldschmidt, and M. Hafezi, “A topological source of quantum light,” Nature 561, 502–506 (2018).
  * Cheng  _et al._ (2016) X. Cheng, C. Jouvaud, X. Ni, S. H. Mousavi, A. Z. Genack, and A. B. Khanikaev, “Robust reconfigurable electromagnetic pathways within a photonic topological insulator,” Nat. Mater. 15, 542–548 (2016).
  * St-Jean  _et al._ (2017) P. St-Jean, V. Goblot, E. Galopin, A. Lemaître, T. Ozawa, L. L. Gratiet, I. Sagnes, J. Bloch, and A. Amo, “Lasing in topological edge states of a one-dimensional lattice,” Nat. Photonics 11, 651–656 (2017).
  * Bahari  _et al._ (2017) B. Bahari, A. Ndao, F. Vallini, A. El Amili, Y. Fainman, and B. Kanté, “Nonreciprocal lasing in topological cavities of arbitrary geometries,” Science 358, 636–640 (2017).
  * Bandres  _et al._ (2018) M. A. Bandres, S. Wittek, G. Harari, M. Parto, J. Ren, M. Segev, D. N. Christodoulides, and M. Khajavikhan, “Topological insulator laser: Experiments,” Science 359, eaar4005 (2018).
  * Zeng  _et al._ (2020) Y. Zeng, U. Chattopadhyay, B. Zhu, B. Qiang, J. Li, Y. Jin, L. Li, A. G. Davies, E. H. Linfield, B. Zhang, Y. Chong, and Q. J. Wang, “Electrically pumped topological laser with valley edge modes,” Nature 578, 246–250 (2020).
  * Smirnova  _et al._ (2020a) D. Smirnova, A. Tripathi, S. Kruk, M. Hwang, H. Kim, H. Park, and Y. Kivshar, “Room-temperature lasing from nanophotonic topological cavities,” Light Sci. Appl. 9, 127 (2020a).
  * Leykam and Chong (2016) D. Leykam and Y. D. Chong, “Edge solitons in nonlinear-photonic topological insulators,” Phys. Rev. Lett. 117, 143901 (2016).
  * Kartashov and Skryabin (2016) Y. V. Kartashov and D. V. Skryabin, “Modulational instability and solitary waves in polariton topological insulators,” Optica 3, 1228–1236 (2016).
  * Mukherjee and Rechtsman (2020) S. Mukherjee and M. C. Rechtsman, “Observation of floquet solitons in a topological bandgap,” Science 368, 856–859 (2020).
  * Mukherjee and Rechtsman (2021) S. Mukherjee and M. C. Rechtsman, “Observation of unidirectional solitonlike edge states in nonlinear floquet topological insulators,” Phys. Rev. X 11, 041057 (2021).
  * Pernet  _et al._ (2022) Nicolas Pernet, Philippe St-Jean, Dmitry D Solnyshkov, Guillaume Malpuech, Nicola Carlon Zambon, Quentin Fontaine, Bastian Real, Omar Jamadi, Aristide Lemaître, Martina Morassi, Luc Le Gratiet, Téo Baptiste, Abdelmounaim Harouri, Isabelle Sagnes, Alberto Amo, Sylvain Ravets, and Jacqueline Bloch, “Gap solitons in a one-dimensional driven-dissipative topological lattice,” Nature Physics 18, 678–684 (2022).
  * Kruk  _et al._ (2019) S. Kruk, A. Poddubny, D. Smirnova, L. Wang, A. Slobozhanyuk, A. Shorokhov, I. Kravchenko, B. Luther-Davies, and Y. Kivshar, “Nonlinear light generation in topological nanostructures,” Nat. Nanotechnol. 14, 126–130 (2019).
  * Wang  _et al._ (2019) Y. Wang, L. J. Lang, C. H. Lee, B. Zhang, and Y. D. Chong, “Topologically enhanced harmonic generation in a nonlinear transmission line metamaterial,” Nat. Commun. 10, 1102 (2019).
  * Bardyn  _et al._ (2016) C. E. Bardyn, T. Karzig, G. Refael, and T. C. H. Liew, “Chiral bogoliubov excitations in nonlinear bosonic systems,” Phys. Rev. B 93, 020502 (2016).
  * Hadad  _et al._ (2018) Y. Hadad, J. C. Soric, A. B. Khanikaev, and A. Alù, “Self-induced topological protection in nonlinear circuit arrays,” Nat. Electron. 1, 178 (2018).
  * Zhou  _et al._ (2017) Xin Zhou, You Wang, Daniel Leykam, and Y D Chong, “Optical isolation with nonlinear topological photonics,” New J. Phys. 19, 095002 (2017).
  * Mandal  _et al._ (2019) S. Mandal, R. Ge, and T. C. H. Liew, “Antichiral edge states in an exciton polariton strip,” Phys. Rev. B 99, 115423 (2019).
  * Zangeneh-Nejad and Fleury (2019) F. Zangeneh-Nejad and R. Fleury, “Nonlinear second-order topological insulators,” Phys. Rev. Lett. 123, 053902 (2019).
  * Maczewsky  _et al._ (2020) L. J. Maczewsky, M. Heinrich, M. Kremer, S. K. Ivanov, M. Ehrhardt, F. Martinez, Y. V. Kartashov, V. V. Konotop, L. Torner, D. Bauer, and A. Szameit, “Nonlinearity-induced photonic topological insulator,” Science 370, 701–704 (2020).
  * Lumer  _et al._ (2013) Y. Lumer, Y. Plotnik, M. C. Rechtsman, and M. Segev, “Self-localized states in photonic topological insulators,” Phys. Rev. Lett. 111, 243905 (2013).
  * Kartashov and Skryabin (2017) Y. V. Kartashov and D. V. Skryabin, “Bistable topological insulator with exciton-polaritons,” Phys. Rev. Lett. 119, 253904 (2017).
  * Banerjee  _et al._ (2020) R. Banerjee, S. Mandal, and T. C. H. Liew, “Coupling between exciton-polariton corner modes through edge states,” Phys. Rev. Lett. 124, 063901 (2020).
  * Xia  _et al._ (2021) S. Xia, D. Kaltsas, D. Song, I. Komis, J. Xu, A. Szameit, H. Buljan, K. G. Makris, and Z. Chen, “Nonlinear tuning of pt symmetry and non-hermitian topological states,” Science 372, 72–76 (2021).
  * Kirsch  _et al._ (2021) M. S. Kirsch, Y. Zhang, M. Kremer, L. J. Maczewsky, S. K. Ivanov, Y. V. Kartashov, L. Torner, D. Bauer, A. Szameit, and M. Heinrich, “Nonlinear second-order photonic topological insulators,” Nat. Phys. 17, 995–1000 (2021).
  * Ezawa (2021) M. Ezawa, “Nonlinearity-induced transition in the nonlinear su-schrieffer-heeger model and a nonlinear higher-order topological system,” Phys. Rev. B 104, 235420 (2021).
  * Smirnova  _et al._ (2020b) D. Smirnova, D. Leykam, Y. Chong, and Y. Kivshar, “Nonlinear topological photonics,” Appl. Phys. Rev. 7, 021306 (2020b).
  * Yu  _et al._ (2019) Sunkyu Yu, Xianji Piao, and Namkyoo Park, “Neuromorphic functions of light in parity-time-symmetric systems,” Advanced Science 6, 1900771 (2019).
  * Xie  _et al._ (2018) B. Y. Xie, H. F. Wang, H. X. Wang, X. Y. Zhu, J. H. Jiang, M. H. Lu, and Y. F. Chen, “Second-order photonic topological insulator with corner states,” Phys. Rev. B 98, 205147 (2018).
  * Moreno  _et al._ (2012) I. Moreno, J. A. Davis, T. M Hernandez, D. M. Cottrell, and D. Sand, “Complete polarization control of light from a liquid crystal spatial light modulator,” Opt. Express 20, 364–376 (2012).
  * Töpfer  _et al._ (2021) J. D. Töpfer, I. Chatzopoulos, H. Sigurdsson, T. Cookson, Y. G. Rubo, and P. G. Lagoudakis, “Engineering spatial coherence in lattices of polariton condensates,” Optica 8, 106–113 (2021).
  * Dormand and Prince (1980) J. R. Dormand and P. J. Prince, “A family of embedded runge-kutta formulae,” J. Comput. Appl. Math. 6, 19–26 (1980).
  * Pitaevskii and Stringari. (2003) L. P. Pitaevskii and S. Stringari., “Bose Einstein Condensation,” (Oxford University Press, Oxford, England) (2003).
  * Benalcazar  _et al._ (2017a) Wladimir A. Benalcazar, B. Andrei Bernevig, and Taylor L. Hughes, “Quantized electric multipole insulators,” Science 357, 61–66 (2017a).
  * Liu  _et al._ (2019) Tao Liu, Yu-Ran Zhang, Qing Ai, Zongping Gong, Kohei Kawabata, Masahito Ueda, and Franco Nori, “Second-order topological phases in non-hermitian systems,” Phys. Rev. Lett. 122, 076801 (2019).
  * Hu  _et al._ (2021) Z. Hu, D. Bongiovanni, D. Jukić, E. Jajtić, S. Xia, D. Song, J. Xu, R. Morandotti, H. Buljan, and Z. Chen, “Nonlinear control of photonic higher-order topological bound states in the continuum,” Light Sci. Appl. 10, 164 (2021).
  * Cheng  _et al._ (2022) Dali Cheng, Bo Peng, Meng Xiao, Xianfeng Chen, Luqi Yuan, and Shanhui Fan, “Truncation-dependent $\mathcal{PT}$ phase transition for the edge states of a two-dimensional non-hermitian system,” Phys. Rev. B 105, L201105 (2022).
  * Benalcazar  _et al._ (2017b) Wladimir A. Benalcazar, B. Andrei Bernevig, and Taylor L. Hughes, “Electric multipole moments, topological multipole moment pumping, and chiral hinge states in crystalline insulators,” Phys. Rev. B 96, 245115 (2017b).
  * Luo and Zhang (2019) Xi-Wang Luo and Chuanwei Zhang, “Higher-order topological corner states induced by gain and loss,” Phys. Rev. Lett. 123, 073601 (2019).
  * Zograf  _et al._ (2022) George Zograf, Kirill Koshelev, Anastasiia Zalogina, Viacheslav Korolev, Richard Hollinger, Duk-Yong Choi, Michael Zuerch, Christian Spielmann, Barry Luther-Davies, Daniil Kartashov, Sergey V. Makarov, Sergey S. Kruk, and Yuri Kivshar, “High-harmonic generation from resonant dielectric metasurfaces empowered by bound states in the continuum,” ACS Photonics 9, 567–574 (2022).
  * Smirnova  _et al._ (2019) Daria Smirnova, Sergey Kruk, Daniel Leykam, Elizaveta Melik-Gaykazyan, Duk-Yong Choi, and Yuri Kivshar, “Third-harmonic generation in photonic topological metasurfaces,” Phys. Rev. Lett. 123, 103901 (2019).
  * Choi  _et al._ (2022) Seou Choi, Jungmin Kim, Jeonghun Kwak, Namkyoo Park, and Sunkyu Yu, “Topologically protected all-optical memory,” Adv. Electron. Mater. 8, 2200579 (2022).
  * Yu  _et al._ (2021) Sunkyu Yu, Xianji Piao, and Namkyoo Park, “Topologically protected optical signal processing using parity–time-symmetric oscillation quenching,” Nanophotonics 10, 2883–2891 (2021).
  * Khomeriki and Leon (2005) R. Khomeriki and J. Leon, “Light detectors bistable nonlinear waveguide arrays,” Phys. Rev. Lett. 94, 243902 (2005).
  * Solntsev  _et al._ (2014) Alexander S. Solntsev, Frank Setzpfandt, Alex S. Clark, Che Wen Wu, Matthew J. Collins, Chunle Xiong, Andreas Schreiber, Fabian Katzschmann, Falk Eilenberger, Roland Schiek, Wolfgang Sohler, Arnan Mitchell, Christine Silberhorn, Benjamin J. Eggleton, Thomas Pertsch, Andrey A. Sukhorukov, Dragomir N. Neshev, and Yuri S. Kivshar, “Generation of nonclassical biphoton states through cascaded quantum walks on a nonlinear chip,” Phys. Rev. X 4, 031007 (2014).
  * Kruse  _et al._ (2013) R. Kruse, F. Katzschmann, A. Christ, A. Schreiber, S. Wilhelm, K. Laiho, A. Gábris, C. S. Hamilton, I. Jex, and C. Silberhorn, “Spatio-spectral characteristics of parametric down-conversion in waveguide arrays,” New J. Phys. 15, 083046 (2013).
  * Peano  _et al._ (2016) V. Peano, M. Houde, C. Brendel, F. Marquardt, and A. A. Clerk, “Topological phase transitions and chiral inelastic transport induced by the squeezing of light,” Nat. Commun. 7, 10779 (2016).
  * Yang  _et al._ (2021) Y. Yang, Z. Zuo, and D. Cao, “Biphoton topology in a quadratic nonlinear waveguide array under the su-schrieffer-heeger model,” Phys. Rev. A 104, 043710 (2021).
  * Doyle  _et al._ (2022) Cooper Doyle, Wei-Wei Zhang, Michelle Wang, Bryn A. Bell, Stephen D. Bartlett, and Andrea Blanco-Redondo, “Biphoton entanglement of topologically distinct modes,” Phys. Rev. A 105, 023513 (2022).
  * Bardyn  _et al._ (2015) Charles-Edouard Bardyn, Torsten Karzig, Gil Refael, and Timothy C. H. Liew, “Topological polaritons and excitons in garden-variety systems,” Phys. Rev. B 91, 161413 (2015).
  * Nalitov  _et al._ (2015) A. V. Nalitov, D. D. Solnyshkov, and G. Malpuech, “Polariton $\mathbb{Z}$ topological insulator,” Phys. Rev. Lett. 114, 116401 (2015).
  * Klembt  _et al._ (2018) S. Klembt, T. H. Harder, O. A. Egorov, K. Winkler, R. Ge, M. A. Bandres, M. Emmerling, L. Worschech, T. C. H. Liew, M. Segev, C. Schneider, and S. Höfling, “Exciton-polariton topological insulator,” Nature 562, 552–556 (2018).
  * Banerjee  _et al._ (2018) R. Banerjee, T. C. H. Liew, and O. Kyriienko, “Realization of hofstadter’s butterfly and a one-way edge mode in a polaritonic system,” Phys. Rev. B 98, 075412 (2018).
  * Kartashov and Skryabin (2019) Yaroslav V. Kartashov and Dmitry V. Skryabin, “Two-dimensional topological polariton laser,” Phys. Rev. Lett. 122, 083902 (2019).
  * Mandal  _et al._ (2020) S. Mandal, R. Banerjee, Elena A. Ostrovskaya, and T. C. H. Liew, “Nonreciprocal transport of exciton polaritons in a non-hermitian chain,” Phys. Rev. Lett. 125, 123902 (2020).
  * Mandal  _et al._ (2022) Subhaskar Mandal, Rimi Banerjee, and Timothy C. H. Liew, “From the topological spin-hall effect to the non-hermitian skin effect in an elliptical micropillar chain,” ACS Photonics 9, 527–539 (2022), https://doi.org/10.1021/acsphotonics.1c01425 .
  * Harder  _et al._ (2021) Tristan H. Harder, Meng Sun, Oleg A. Egorov, Ihor Vakulchyk, Johannes Beierlein, Philipp Gagel, Monika Emmerling, Christian Schneider, Ulf Peschel, Ivan G. Savenko, Sebastian Klembt, and Sven Höfling, “Coherent topological polariton laser,” ACS Photonics 8, 1377–1384 (2021), https://doi.org/10.1021/acsphotonics.0c01958 .
  * Baas  _et al._ (2004) A. Baas, J. Ph. Karr, H. Eleuch, and E. Giacobino, “Optical bistability in semiconductor microcavities,” Phys. Rev. A 69, 023809 (2004).
  * Paraïso  _et al._ (2010) T. K. Paraïso, M. Wouters, Y. Léger, F. Morier-Genoud, and B. Deveaud-Plédran, “Multistability of a coherent spin ensemble in a semiconductor microcavity,” Nat. Mater. 9, 655–660 (2010).
  * Banerjee and Liew (2020) R. Banerjee and T. C. H. Liew, “Artificial life in an exciton-polariton lattice,” New J. Phys. 22, 103062 (2020).
  * Utsunomiya  _et al._ (2008) S. Utsunomiya, L. Tian, G. Roumpos, C. W. Lai, N. Kumada, T. Fujisawa, M. Kuwata-Gonokami, A. Löffler, S. Höfling, A. Forchel, and Y. Yamamoto, “Observation of Bogoliubov excitations in exciton-polariton condensates,” Nat. Phys. 4, 700–705 (2008).
  * Carusotto and Ciuti (2013) I. Carusotto and C. Ciuti, “Quantum fluids of light,” Rev. Mod. Phys. 85, 299 (2013).
  * Harder  _et al._ (2020) T. H. Harder, O. A. Egorov, J. Beierlein, P. Gagel, J. Michl, M. Emmerling, C. Schneider, U. Peschel, S. Höfling, and S. Klembt, “Exciton-polaritons in flatland: Controlling flatband properties in a lieb lattice,” Phys. Rev. B 102, 121302 (2020).
  * Dusel  _et al._ (2020) M. Dusel, S. Betzold, O. A. Egorov, S. Klembt, J. Ohmer, U. Fischer, S. Höfling, and C. Schneider, “Room temperature organic exciton–polariton condensate in a lattice,” Nat. Commun. 11, 2863 (2020).
  * Sun  _et al._ (2017) Y. Sun, Y. Yoon, M. Steger, G. Liu, L. N. Pfeiffer, K. West, D. W. Snoke, and K. A. Nelson, “Direct measurement of polariton–polariton interaction strength,” Nat. Phys. 13, 870–875 (2017).
  * Estrecho  _et al._ (2019) E. Estrecho, T. Gao, N. Bobrovska, D. Comber-Todd, M. D. Fraser, M. Steger, K. West, L. N. Pfeiffer, J. Levinsen, M. M. Parish, T. C. H. Liew, M. Matuszewski, D. W. Snoke, A. G. Truscott, and E. A. Ostrovskaya, “Direct measurement of polariton-polariton interaction strength in the thomas-fermi regime of exciton-polariton condensation,” Phys. Rev. B 100, 035306 (2019).
  * Schneider  _et al._ (2013) C. Schneider, A. Rahimi-Iman, N. Y. Kim, J. Fischer, I. G Savenko, M. Amthor, M. Lermer, A. Wolf, L. Worschech, V. D Kulakovskii, I. A Shelykh, M. Kamp, S. Reitzenstein, A. Forchel, Y. Yamamoto, and S. Höfling, “An electrically pumped polariton laser,” Nature 497, 348–352 (2013).
  * Klaas  _et al._ (2017) M. Klaas, S. Mandal, T. C. H. Liew, M. Amthor, S. Klembt, L. Worschech, C. Schneider, and S. Höfling, “Optical probing of the coulomb interactions of an electrically pumped polariton condensate,” Applied Phys. Lett. 110, 151103 (2017).
  * Gagel  _et al._ (2022) Philipp Gagel, Tristan H. Harder, Simon Betzold, Oleg A. Egorov, Johannes Beierlein, Holger Suchomel, Monika Emmerling, Adriana Wolf, Ulf Peschel, Sven Höfling, Christian Schneider, and Sebastian Klembt, “Electro-optical switching of a topological polariton laser,” ACS Photonics 9, 405–412 (2022), https://doi.org/10.1021/acsphotonics.1c01605 .