ELGC-Net: Efficient Local-Global Context Aggregation for Remote Sensing Change Detection

Mubashir Noman^1
Mustansar Fiaz^2
Hisham Cholakkal^1
Salman Khan^1,3
and Fahad Shahbaz Khan^1,4

^1Mohamed bin Zayed University of AI   ^2IBM Research   ^3Australian National University
^4Linköping University

Deep learning has shown remarkable success in remote sensing change detection (CD), aiming to identify semantic change regions between co-registered satellite image pairs acquired at distinct time stamps. However, existing convolutional neural network (CNN) and transformer-based frameworks often struggle to  accurately segment semantic change regions.  Moreover, transformers-based methods with standard self-attention suffer from quadratic computational complexity with respect to the image resolution, making them less practical for CD tasks with limited training data. To address these issues, we propose an efficient change detection framework, ELGC-Net, which leverages rich contextual information to precisely estimate change regions while reducing the model size. Our ELGC-Net comprises a Siamese encoder, fusion modules, and a decoder. The focus of our design is the introduction of an Efficient Local-Global Context Aggregator (ELGCA) module within the encoder, capturing enhanced global context and local spatial  information through a novel pooled-transpose (PT) attention and depthwise convolution, respectively. The PT attention employs  pooling operations for robust feature extraction and minimizes computational cost with transposed attention.
Extensive experiments on three challenging CD datasets demonstrate that ELGC-Net outperforms existing methods.

Compared to the recent transformer-based CD approach (ChangeFormer), ELGC-Net achieves a 1.4% gain in intersection over union (IoU) metric on the LEVIR-CD dataset, while significantly reducing trainable parameters.
Our proposed ELGC-Net sets a new state-of-the-art performance in remote sensing change detection benchmarks.  Finally, we also introduce ELGC-Net-LW, a  lighter variant with significantly reduced computational complexity,   suitable for resource-constrained settings, while achieving comparable performance.
Our source code is publicly available at <https://github.com/techmn/elgcnet>.

Remote sensing, change detection, transformers, local and global context aggregation.

§ INTRODUCTION

Remote  sensing change detection is a challenging task where the objective is to detect semantic changes between satellite image pairs taken at different time stamps. It finds diverse applications such as environmental monitoring, urban planning <cit.>, land cover mapping, and disaster assessment <cit.>. With the increasing availability of high-resolution satellite imagery, the demand for efficient and accurate change detection methods to process this huge amount of data has grown significantly in recent years.

Remote sensing change detection faces several challenges, such as environmental variations, complex and irregular geometrical structures, variations in object scales <cit.>, and the high-resolution nature of satellite imagery. These challenges introduce various complexities when accurately identifying changes between image pairs. For instance, Fig. <ref> (columns 1-3) illustrates semantically irrelevant changes that need to be ignored by the change detection (CD) algorithm (indicated by red boxes), such as (i) shadows and illumination variations, (ii) displacement of movable objects like cars and roof color changes, and (iii) seasonal variations in land coverings, while (iv) accurately segmenting both subtle and major semantic changes between image pairs (indicated by green boxes in the last column).
Thus, it is desired to develop an accurate CD method that can effectively learn the semantic information and accurately detect the change regions while suppressing irrelevant ones.
Traditional approaches utilize feature difference technique and obtain binary mask via thresholding  <cit.>. However, such approaches are unable to capture the underlying characteristics of the objects and fail to discriminate between the semantic and noisy changes. Numerous methods have been proposed that are based on hand-crafted features and utilize different classifiers including decision trees <cit.>,   change vector analysis <cit.>, support vector machine <cit.>, and clustering approaches <cit.> for better discrimination.
Nevertheless, these classical approaches exhibit various limitations and are incapable of capturing rich feature representations.

Example image pairs depicting various challenges related to remote sensing change detection and the corresponding ground-truth segmentation masks. It is highly challenging to ignore semantically irrelevant changes (shown in red boxes) such as (i) shadows and illumination variations, (ii) movable cars and building roof changes (iii) seasonal variations, while accurately segmenting   (iv) subtle and large semantic changes (shown in green box).

On the other hand,  convolutional neural network (CNN) based CD methods <cit.> utilize Siamese architecture to address the CD problem. These approaches typically exploit architectural components such as dilated convolutions, multi-scale features, and attention mechanisms in their framework. However, these CNN-based approaches often struggle to capture global contextual representations, which is likely to deteriorate the CD performance.

In contrast, transformer-based approaches <cit.>
employing standard self-attention <cit.> often struggles to effectively incorporate local contextual information for accurate segmentation of subtle changes at irregular boundaries.
Furthermore, while state-of-the-art transformer-based change detection frameworks <cit.> have achieved impressive results, they come with inherent challenges. These approaches generally rely on standard self-attention, which exhibits quadratic complexity with the number of tokens in the encoder. Consequently, the requirement for a large number of parameters, extensive memory usage, and high FLOPs (floating-point operations per second) make these methods less suitable for practical remote sensing change detection applications.

Hence, it is desired to capture both local and global contextual information to effectively detect subtle as well as significant structural changes between image pairs.

In order to handle the above mentioned challenges, we propose a Siamese-based efficient change detection framework named  ELGC-Net to accurately identify semantic changes in satellite image pairs taken at different time stamps.
The focus of our design is to the introduction of a novel efficient local-global context aggregator (ELGCA) module that can capture both local and global contextual information while addressing the limitations of standard self-attention.
Specifically, our ELGCA strives to (i) capture global contextual information through pooled-transpose (PT) attention, (ii) encode local spatial context using depth-wise convolution, and (iii) multi-channel feature aggregation.  The PT attention leverages pooling operations for robust feature extraction and reduces the computational cost with transposed attention.
We extensively evaluate ELGC-Net on three diverse change detection benchmark datasets, validating its superiority over existing methods in literature.

Our method achieves state-of-the-art performance while significantly reducing the model size, as demonstrated by the comprehensive quantitative and qualitative analysis.
In addition to the  ELGC-Net, we introduce its lightweight variant called ELGC-Net-LW that requires only 6.78M parameters while achieving comparable performance to recent ChangeFormer <cit.>, which utilizes 41M parameters (as shown in Fig. 2 in the manuscript). ELGC-Net-LW serves as a resource-efficient alternative, catering to scenarios with limited computational resources without compromising on CD accuracy. The key contributions of our method are summarized as:

* We propose a Siamese-based efficient change detection framework named as ELGC-Net to accurately identify semantic changes between the bi-temporal satellite images.

* Instead of modeling dense dependencies, our ELGCA module proposes to capture both local and global relationships through depthwise convolution and pooled-transpose attention, respectively, while reducing computational complexity.

* Our extensive experiments on three CD datasets demonstrate the robustness and efficiency of our approach while achieving state-of-the-art performance. Moreover, our lightweight variant called ELGC-Net-LW obtains comparable performance compared to recent ChangeFormer <cit.> requiring only 6.78M parameters.

The remainder of this paper is structured as follows. Section II presents the related works and Section III discusses about the baseline framework. Section IV presents an in-depth discussion of the proposed method and our Efficient Local-Global Context Aggregator (ELGCA) module. Section V presents the experimental setup, results, and performance comparison with existing methods, and finally, Section VI concludes the paper.

axis lines = left,
ymin=76, ymax=85,
xmin=0, xmax=200,
xlabel= No. of Parameters (M),
ylabel = (Change Class) IoU %,
grid style=dashed,
(legend) at (axis description cs:1.02,0.07);

[only marks,
mark=otimes*, violet,
mark size=3.5pt

[only marks,
mark=otimes*, green,
mark size=3.5pt

[only marks,
mark=otimes*, pink,
mark size=3.5pt

[only marks,
mark=otimes*, blue,
mark size=3.5pt

[only marks,
mark=otimes*, cyan,
mark size=3.5pt

[only marks,
mark=otimes*, orange,
mark size=3.5pt

[only marks,
mark=otimes*, olive,
mark size=3.5pt

[only marks,
mark=halfdiamond*, magenta,
mark size=5pt
(6.78, 82.36);

[only marks,
mark=pentagon*, magenta,
mark size=5pt
(10.57, 83.83);

[draw=none,fill=white, anchor= south east] at

Method     Params (M)     IoU     F1
<ref> STANet <cit.>     16.93     77.40     87.30

<ref> SNUNet <cit.>     12.03     78.83     88.18

<ref> BiT <cit.>     18.78     80.68     89.31

<ref> DTCDSCN <cit.>     41.07     78.05     87.67

<ref> H-TransCD <cit.>     183.83     81.92     90.06

<ref> ChangeFormer <cit.>     41.03     82.48     90.40

<ref> TransUNetCD <cit.>     95.49     83.67     91.11

<ref> ELGC-Net-LW (Ours)     6.78     82.36     90.33

<ref> ELGC-Net (Ours)     10.57     83.83     91.20

Accuracy (IoU) vs. model size (params) comparison with existing methods on  on LEVIR-CD.  Our ELGCNet achieves state-of-the-art performance

while having $9\times$ less number of model parameters compared to the existing state-of-the-art TransUNetCD <cit.>.

Overall architecture of the proposed CD framework. (a) The complete network architecture is presented, illustrating the  pre- and post-change input images to the shared encoder. The encoder blocks extract features at four stages and these features are merged by a simple fusion module comprising linear projection, feature concatenation, and a $1\times1$ convolution.
These fused feature maps at each stage are then merged in the decoder, which includes several convolutions and transposed convolution layers to obtain upsampled features. The upsampled feature maps are used in the prediction layer to obtain the final change map. (b) The structure of our Encoder block is shown, featuring the Efficient Local-Global Context Aggregation (ELGCA) module and a convolutional MLP. (c) A detailed view of the proposed ELGCA module performing the following key operations: (i) capturing local spatial context  using a 3x3 depth-wise convolution ($\bar X^i_{lo}$)
and (ii) global context aggregation ($A^i_{att}$) through a pooled-transpose (PT) attention operation,
and multi-channel feature aggregation using a 1x1 convolution ($Z^i$).
To enhance the efficiency of our PT attention,
we perform transposed attention ($G$) having linear complexity with the number of tokens, on pooled $Q^i$ and $K^i$ tokens (denoted as $\bar{Q}^i$, $\bar{K}^i$) on a sub-set of channels ($C/4$).
The aforementioned operations within our ELGCA module are performed in parallel on different groups (subsets) of channels obtained through channel splitting, leading to improved computational efficiency.

§ RELATED WORK

As mentioned earlier, the automatic CD has been widely studied over the past decade and various CD methods have been proposed <cit.>. Earlier, conventional methods strive to identify the changes via computing the image difference techniques <cit.>, in the bi-temporal images. To do so, these approaches typically generate a change map by identifying the pixels that have substantial differences and then employ a threshold to finally obtain the segmentation mask. Nevertheless, these approaches exhibit few benefits due to the incapabilities of the RGB features to discriminate between semantic and irrelevant changes. In addition, different classification approaches have been proposed which employ hand-crafted features <cit.>  and discriminate the change features from the temporal images via considering statistical modeling <cit.>, decision tree <cit.>,   change vector analysis <cit.>, support vector machine <cit.>, and clustering approaches <cit.>.

Later, deep learning approaches based on deep Convolutional Neural Networks (CNNs)   have shown significant performance improvement on numerous computer vision task such as segmentation <cit.>, object tracking <cit.>, action recognition <cit.> and many more.

The CNNs have shown their capability to better capture the underlying characteristics and are the key factor in utilizing them for the remote sensing change detection field <cit.>.
Various deep learning methods have also been explored in the field of remote sensing change detection <cit.>.
Ham et al. <cit.>  used a deconvolutional network as a semantic segmentation algorithm to detect unregistered buildings from UAVs (Unmanned Aerial Vehicles) images.
Zhao et al. <cit.> and Liu et al. <cit.> employed Mask-RCNN method for semantic segmentation of buildings and outdoor sports venues, respectively.  BiDateNet <cit.> proposed a  deep learning framework based on U-Net architecture to obtain a change map.
Guo et al. <cit.> configured a Fully Convolutional Network in the form of a Siamese Network to detect changes in continuous images. Khan et al. <cit.> used a CNN-based weakly supervised technique to directly detect and localize the change detection for the input image pairs.  CDNet <cit.>, FC-EF, FC-Siam-diff, and FC-Siam-Conc <cit.> employed boundary decisions to obtain the change map.

Although CNN methods exhibit better performance in the field, they strive to capture the long-range dependencies whereas transformers may struggle in capturing accurate localization due to lack of fine details <cit.>.
This motivated the researchers to combine the large receptive fields of transformers with the local contextual features of convolutions to obtain discriminative feature representations <cit.>.
Jiang et al. <cit.> introduce a co-attention module to enhance the features for better detection of change regions. The authors utilize the Siamese network in encoder-decoder architecture and use the VGG-16 network as a backbone feature extractor. Afterward, extracted features are enhanced by utilizing the co-attention module and input to the decoder for change map prediction.
Lie et al. <cit.> propose a hybrid transformer model that uses a Siamese CNN encoder to extract features and concatenate the encoded feature of two streams. Then, they capture the global contextual information using a transformer and use a cascading up-sampling decoder to predict the change map.
Chen et al. <cit.> propose a bi-temporal image transformer (BIT) network to model the long-range context information for CD task. The BIT utilizes the ResNet18 <cit.> to extract features from bi-temporal images followed by a transformer encoder to capture global context relationships in bi-temporal images. A decoder then projects back the encoded features to spatial space and obtains a binary change map.

Song et al. <cit.> introduce a hybrid model based on ResNet18 <cit.> backbone feature extractor and multiscale SwinTransformer <cit.> for enhancing the extracted feature maps.
Ke et al. <cit.> propose a hybrid transformer module to leverage the Siamese-based CNN backbone. Extracted features from the CNN backbone are enhanced using a hybrid transformer and a cascade decoder predicts the change map.
Bandara et al. <cit.> introduce a transformer-based encoder to capture better feature representation at multiple scales. The decoder fuses the multiscale features from the bi-temporal image and predicts the binary change mask.
Chen et al. <cit.> introduce a dual attention mechanism to capture the discriminative features for better detection of CD regions.
Fang et al. <cit.> propose an ensemble channel attention module (ECAM) to aggregate and refine features at multiple levels to get better CD results.
Zhang et al. <cit.> propose a deeply supervised difference discrimination network that takes the extracted bi-temporal features and utilizes a difference discrimination network to predict a change map.
Chen et al. <cit.> propose a basic and pyramid spatio-temporal attention module to obtain better discriminative features for CD task.

Contrary to the above, we propose an efficient global and local context aggregator that combines the features from  different context aggregators in a parallel manner to obtain rich feature representations while maintaining fewer parameters at multi-scale levels that can discriminate the complex change regions with better accuracy.

§ BASE FRAMEWORK

Our base framework is a transformer-based Siamese network consisting of an encoder, fusion modules, and a   decoder. The encoder comprises four stages, each incorporating multi-head self-attention layers followed by a convolutional MLP (multi-layer perceptron) network. This base framework takes pre- and post-change images as input and extracts corresponding feature pairs at different stages using a shared encoder.

At each stage, the extracted features undergo a linear layer transformation before being fed to a fusion module. The fusion module encodes the semantic changes between the image pairs by concatenating the features, denoted as $\hat X^i_{pre}$ and $\hat X^i_{post}$, across the channel dimension. Subsequently, a convolutional layer is employed to reduce the number of channels, followed by an activation layer.

The fused features from all stages are then passed to a decoder, where these features are merged and further processed through several convolutions and transposed convolution layers to obtain higher-resolution feature maps. Finally, these upsampled feature maps are utilized in a prediction layer to obtain the final change map.

§ METHOD

This section describes the overall architecture of the network and explains the ELGCA module in detail.

§.§ Overall Architecture

The overall network architecture of the proposed change detection (CD) framework, named as ELGC-Net, is illustrated in Fig. <ref>-(a). Similar to the baseline framework, the proposed CD framework consists of a Siamese encoder, fusion modules, and a   decoder. The Siamese encoder takes a pair of satellite images as input and outputs multi-scale feature maps $\hat X^i_{pre}$ and $\hat X^i_{post}$ at its four stages $i=1,2,3,4$. At each stage, the feature maps are initially down-sampled through a patch embedding layer and then passed through a series of encoder blocks. The structure of our encoder block is shown in Fig. <ref>-(b) which comprises our novel efficient local-global context aggregator (ELGCA) module and a convolutional MLP layer.

The focus of our design is the introduction of an efficient local-global context aggregator (ELGCA) module within the encoder block. Our ELGCA module aims to leverage the benefits of both convolution and self-attention <cit.> to aggregate local and global context, thereby improving the accuracy of predicted change maps while reducing the number of parameters and FLOPs.

At each stage,  similar to the base framework, the feature pairs from the encoder blocks are fused through   fusion modules performing simple linear projection,  feature concatenation, and a $1\times 1$ convolution operation to reduce the number of channels.
These fused features are subsequently passed to a   decoder for change map prediction.

Decoder Architecture:  Similar to the baseline, our ELGC-Net utilizes a decoder composed of several convolution and transpose convolution layers as illustrated in Fig. <ref>. The multi-scale fused features ($\hat X^{i}_{fused}$ where $i \in {1,2,3,4}$) from four stages are concatenated along channel dimension and fed to a $1 \times 1$ convolution layer. Afterwards, we apply the transpose convolution to increase the spatial resolution of the feature maps. Then, a residual block comprising of two $3 \times 3$ convolution layers is utilized to enhance the feature maps. The transpose convolution, cascaded with a residual block, is utilized twice to obtain the same spatial dimension as the input of the model. Finally, we apply a convolution layer to obtain prediction scores consisting of two channels where the first channel corresponds to no-change and the second channel corresponds to change class scores. The binary change map is obtained by using the argmax operation along the channel dimension.

Next, we present a detailed discussion of our ELGCA module.

ELGC-Net decoder takes fused features $\hat X^{i}_{fused}$ from the four stages (where $i = {1,2,3,4}$) and concatenates them along the channel dimension. Then, it utilizes $1 \times 1$ convolution to project features. Afterward, the upsampling is performed twice to obtain the same spatial dimension as the model input using cascaded transpose convolution and a residual block composed of two convolution layers.
Finally, a convolution layer is utilized to obtain the prediction scores having two channels.

§.§ Efficient Local-Global Context Aggregation (ELGCA) Module

Our ELGCA module strives for capturing local and global context information while reducing the computational complexity compared to self-attention through  careful design choices.
Motivated by the inception module <cit.> that  stacks multiple convolution operators of different  kernel sizes  in parallel, we split the channels into different groups and perform diverse operations on each channel split (see Fig. <ref>-(c)).
Specifically, we channel-wise split the input features and input them into two separate context aggregators to  obtain local  and global contextual information.

Let  $X^i \in \mathcal{R}^{H^i \times W^i \times C^i}$ be the input feature to the ELGCA module at  the $i^{th}$ stage, where ($H^i, W^i, C^i$) denote its height, width, and channels, respectively.  We first perform channel-wise splitting of  input feature $X^i$, resulting in  $ X^i_{gl}$,   $X^i_{lo}$ $\in \mathcal{R}^{H^i \times W^i \times \frac{C^i}{2}}$ from which  inputs to our  PT attention  and local context aggregator are obtained. i.e,

\begin{equation}
X^i_{gl}, X^i_{lo} = Split(X^i)
\end{equation}

Next, we describe the PT attention and local  context aggregation within our ELGCA module.

Pooled-transpose (PT)  Attention:
As mentioned earlier, the  standard self-attention has  quadratic computational complexity with respect to number of tokens, whereas the proposed PT attention has only linear complexity with respect to tokens.

We perform   a $1 \times 1$ convolution
on  $X^i_{gl}$  and split  it into $Z^i$, $Q^i$, $K^i$, and $V^i$  features, respectively. Here, $Q^i$, $K^i$, and $V^i$ are used as query, key, and value for our PT attention capturing global contextual information, whereas  $Z^i$ performs   multi-channel feature aggregation. The key steps of our PT attention are  the following:

At first, our PT attention   performs average- and max- pooling operations over the  $Q^i$ and $K^i$ features, respectively to obtain pooled features $\bar{Q}^i$ and
$\bar{K}^i$ that are robust to slight  variations.   Specifically, we perform  $3\times 3$  average (Avg) pooling with stride 2 on  $Q^i$ to capture  mean features $\bar {Q}^i$,  and employ a  $2\times 2$  max (Max) pooling with stride 2 on  $K^i$ to obtain  dominant features $\bar {K}^i$.

Secondly, we employ transposed attention ($G$), between    $\bar {Q}^i$,  $\bar {K}^i$, and $V^i$ embedding,  having linear complexity with respect to the number of tokens on a subset of feature channels ($C/4$). Our PT attention operations  to obtain  feature representation $A^i_{att}$   can be summarized as follows:

\begin{equation}
{A^i}_{att} = V^i  \times [\sigma (\bar{K^i}^T  \times \bar{Q}^i)],
\end{equation}

where $\times$ indicates  matrix multiplication   and $\sigma$ represents the Softmax operations, respectively. Here, $\bar {Q}^i$, $\bar{K}^i$  $\in \mathcal{R}^{\frac{H^i}{2}\frac{W^i}{2} \times \frac{C^i}{4}}$, and  ${V}^i$ $\in \mathcal{R}^{H^iW^i \times \frac{C^i}{4}}$.
Our PT attention leads to improved performance and generalization capabilities,  while simultaneously reducing the convergence time and computational cost.

Local Contextual Aggregation:
In order to capture spatially and channel-wise local contextual information, we perform  $3\times3$  depthwise convolution on input $X^i_{lo}$, resulting in feature representations $\bar X^i_{lo}$. This   depthwise convolution based local context aggregation effectively captures local contextual information with the minimum  number of model parameters and FLOPs, and it alleviates the requirement of positional embedding <cit.>.

Finally, the local and global context aggregated features  with various  receptive fields (such as $\bar X^i_{lo}$, $Z^i$, and $A^i_{att}$) are merged using concatenation operation to  obtain enriched local-global contextual aggregated feature ($\bar X^i$). In contrast to <cit.>, where context aggregators are employed in a sequential manner,  our ELGCA module exploits multiple context aggregators at different subsets of the channels in a parallel fashion which leads to reduced computational complexity.

In summary, our ELGC-Net CD framework employs the ELGCA module at different stages of the encoder block and generates pre- and post-change feature maps ($\hat X^i_{pre}$ and $\hat X^i_{post}$) from the two streams. These $\hat X^i_{pre}$ and $\hat X^i_{post}$ features at each stage are fused ($\hat X^i_{fuse}$) using a fusion module and then taken as input to a decoder for obtaining the final change map predictions.

State-of-the-art comparison on LEVIR-CD, DSIFN-CD, and CDD-CD. The results are reported in terms of the change class IoU, change class F1, and overall accuracy (OA) metrics. Our ELGC-Net efficiently combines the global and local contextual features and achieves superior performance in terms of all metrics over these datasets. Here, the best two results are highlighted in red and blue colors, respectively.

1|l|2*Method     1c|2*Params (M)      3c|LEVIR-CD     3c|DSIFN-CD     3c|CDD-CD

IoU     F1      OA
IoU     F1      OA
IoU     F1      OA

BIT  <cit.>
80.68     89.31      98.92

51.22     67.74     89.72
80.01     88.90     97.47

FC-Siam-diff <cit.>

75.92     86.31     98.67

48.44     65.26     89.06
54.57     70.61     94.95

FC-Siam-conc <cit.>

71.96     83.69     98.49

44.33     61.43     84.51
60.14     75.11     94.95

DTCDSCN <cit.>
78.05     87.67     98.77

48.46     65.29     88.14
85.34     92.09     98.16

SNUNet <cit.>
78.83     88.16     98.82

47.13     64.07     86.54
72.11     83.89     96.23

STANet <cit.>
77.40     87.30     98.66

45.12     62.19     85.00
72.22     84.12     96.13

H-TransCD <cit.>
81.92     90.06     99.00
–     –     –
–     –     –

IFNet   <cit.>
78.77     88.13      98.87
43.28     60.42     84.41
71.91     90.30     97.71

MSTDSNet   <cit.>
78.73     88.10      98.56
–     –     –
–     –     –

TransUNetCD <cit.>
83.67     91.11     –

–     –     –
94.50     97.17     –

ChangeFormer  <cit.>
82.48     90.40      99.04

53.26     69.50     90.56
81.53     89.83     97.68

ELGC-Net-LW (Ours)
82.36     90.33     99.03

57.92     73.35     91.57
93.48     96.63     99.21

ELGC-Net (Ours)
83.83     91.20       99.12

58.42      73.76      92.68
94.50      97.17      99.33

Qualitative results on the LEVIR-CD dataset. We present comparison with the best five existing change detection methods in literature, whose codebases are publicly available. The highlighted region shows that our method (ELGC-Net) is better at detecting the change regions as compared to  CNN-based including FC-Siam-diff <cit.>, STANet <cit.>, DTCDSCN <cit.> and transformer-based such as BIT <cit.> and ChangeFormer <cit.> methods.

§ EXPERIMENTS

§.§ Datasets and Evaluation Protocols

In our experiments, we use three public remote sensing change detection datasets to verify the effectiveness of the proposed network.

§.§.§ LEVIR-CD <cit.>
is a large-scale building change detection dataset. It consists of 637 high-resolution (0.5 m/pixel) image pairs of spatial size $1024 \times 1024$ taken from Google Earth. The time span of these image pairs ranges from 5 to 14 years. The dataset focuses on the changes related to building construction and demolition. In our experiments, we used the cropped version of the dataset having non-overlapping cropped patches of size $256 \times 256$. Following  <cit.>, we utilize the default data split of size 7120, 1024, and 2048 image pairs for train, validation and test sets respectively.

§.§.§ DSIFN-CD <cit.>
is another binary change detection dataset that is collected from Google Earth. The dataset contains 6 high-resolution bi-temporal images covering six cities of China. The dataset is available in cropped sub-image pairs of size $256 \times 256$ with train, validation, and test split of 14400, 1360, and 28 image pairs, respectively.

§.§.§ CDD-CD <cit.>
is a well-known remote sensing change detection dataset that contains 11 seasonal varying image pairs. Among these, 7 image pairs have a spatial size of $4725 \times 2700$ pixels while 4 image pairs are of size $1900 \times 1000$ pixels. The dataset has been used in the clipped form (<cit.>) of size $256 \times 256$ pixels and contains 10000, 3000, and 3000 image pairs in train, validation, and test sets respectively.

Evaluation Protocol:
Similar to the  <cit.>, we evaluate the change detection results  using  change class intersection over union (IoU), change class F1-score, and overall accuracy (OA) protocols on all the datasets.

Qualitative results on the DSIFN-CD dataset with a comparison presented with the best five existing change detection methods in the literature, whose codebases are publicly available. The highlighted region in red box shows that ELGC-Net is better to identify the change regions as compared to the CNN-based FC-Siam-diff <cit.>, STANet <cit.>, DTCDSCN <cit.> as well as transformer-based BIT <cit.> and ChangeFormer <cit.> methods.

§.§ Implementation Details

The encoder of the proposed ELGC-Net is composed of four stages having 3,3,4,3 encoder blocks and outputs feature maps that have channel dimensions of 64, 96, 128, and 256, respectively. The input to the encoder is a pair of co-registered images with the size of $256 \times 256 \times 3$. The decoder outputs the binary change map of the same spatial resolution.

Furthermore, random flip, scale crop, color jitter, and Gaussian blur data augmentations are used in the training process. The network is trained on 4 Nvidia A100 GPUs and model weights are randomly initialized. We set the learning rate equal to 3.1e-4 and train the model for 300 epochs. The value of the learning rate is selected empirically. We linearly decay the learning rate till the last epoch. Following <cit.>, we use AdamW optimizer having weight decay of 0.01 and beta values of (0.9, 0.999).

During inference, we freeze the model weights and pass the input image pair to the model. The network outputs the probability change map having two channels. A binary change map is obtained through the argmax operation along the channel dimension.

Loss Function:  During training, we utilize pixel-wise cross-entropy loss to measure the performance of the network. For binary change detection problem, cross-entropy loss function is mathematically expressed as:

\begin{equation}
Loss = -\frac{1}{N}\sum_{i=1}^{N} {y_i}  {log(\hat {y}_i)} + {(1 - y_i)}  {log(1 - \hat {y}_i)}
\end{equation}

where $y$ is the ground truth class label, $\hat y$ is the prediction probability and $N$ is the number of pixels.

Qualitative results on the CDD-CD dataset. We show the comparison with the best five existing change detection approaches in the literature, whose codebases are publicly available. The highlighted region in the red box shows that ELGC-Net is better to detect the change regions with clear boundaries between the pre- and post-change images.

Ablation study over the  LEVIR-CD. We present the impact of our contributions including local contextual aggregation, multi-channel aggregation, and PT attention (global contextual aggregation) by replacing the self-attention of the baseline. The local contextual aggregation using depthwise convolution (row 1) has the least parameters.
PT attention (row 3) presents a significant improvement over local contextual aggregation  (row 1) and multi-channel aggregation (row 2) in terms of IoU while having increased parameters.
Our final approach (row 7) aggregates all contextual information results in optimal IoU score while having comparable parameters
to local contextual aggregation  (row 1). The best results are shown in bold text.

3|c|Context Aggregators      2*Params (M)     2*IoU     2*F1     2*OA

Local Contextual Aggregation
Multi-channel Aggregation
PT Attention

$\checkmark$              10.43      81.69     89.93     98.98

$\checkmark$         10.72      81.24      89.65     98.97

$\checkmark$     11.32      82.55      90.44     99.02

$\checkmark$     $\checkmark$         10.50     82.18      90.21     99.00

$\checkmark$         $\checkmark$     10.65     83.13       90.79     99.08

$\checkmark$     $\checkmark$     11.02     82.72      90.54     99.04

$\checkmark$     $\checkmark$     $\checkmark$     10.57     83.83      91.20     99.12

Comparison of parameters, FLOPs, memory consumption, and change class IoU metric on LEVIR-CD. ELGC-Net has fewer parameters and FLOPs as compared to the baseline and ChangeFormer.  The best results are shown in bold text.

Method     Params (M)     FLOPs (M)     Memory (MBs)     IoU      F1     OA

ChangeFormer <cit.>     41.03     138803.4     989     82.48     90.40     99.04

Baseline     11.63     138557.9     2087     81.93     90.07     99.00

ELGC-Net (SA)†      10.57     127133.5     1540     83.10     90.77     99.08
ELGC-Net-LW (Ours)     6.78     19815      704     82.36     90.33     99.03

ELGC-Net (Ours)     10.57     123590.5     722.5     83.83     91.20     99.12

† represents PT attention is replaced with self-attention.

§.§ State-of-the-art Quantitative Comparison

We present a quantitative comparison of our approach with state-of-the-art (SOTA) methods over LEVIR-CD, DSIFN-CD, and CDD-CD datasets in Table <ref>.

Tab. <ref> shows that BIT <cit.>, H-TransCD <cit.>, ChangeFormer <cit.>, and TransUNetCD <cit.> perform better on LEVIR-CD dataset by achieving IoU score of 80.68%, 81.92%, 82.48% and 83.67% respectively.

Our ELGC-Net achieves superior performance in terms of all metrics against the compared methods while obtaining the IoU gain of 1.4% and 0.2% over ChangeFormer <cit.> and TransUNetCD <cit.>, respectively. Noted, our method surpasses the state-of-the-art methods without using a pre-trained backbone.
Furthermore, we introduce ELGC-Net-LW which comprises a lightweight decoder. The aim is to further reduce the trainable parameters as the decoder is performing operations on high-resolution feature maps.

ELGC-Net-LW replaces the transpose convolution with bilinear interpolation operation and utilizes a combination of depthwise and standard convolution operations in place of computationally expensive residual modules.
Despite having fewer trainable parameters, our lightweight ELGC-Net-LW also performs better compared to SOTA methods except for TransUNetCD and ChangeFormer which have slightly better IoU scores. However, TransUNetCD and ChangeFormer have significantly higher numbers of trainable parameters as listed in Tab. <ref>.

Similar to the LEVIR-CD, our method demonstrates favorable performance over the DSIFN-CD dataset as shown in Tab. <ref>.

We notice that our ELGC-Net achieves favorable performance against the state-of-the-art approaches.
It is worth mentioning that ELGC-Net-LW, having very few parameters and without using any pre-trained backbone, performs outstandingly on DSIFN-CD.

We also report the state-of-the-art comparison over the CDD-CD dataset in  Tab. <ref>.
We observe that CNN-based DTCDSCN <cit.> obtains a 85.34% IoU score. However, our method achieves a significant gain and provides a 94.50% score of IoU.
On the other hand, transformer-based TransUNetCD <cit.> provides the best performance on CDD-CD by achieving an IoU score of 94.50%. Our ELGC-Net is also obtaining the same score compared to the TransUNetCD <cit.>. Notably, our ELGC-Net does not use any pre-trained backbone feature extractor and still achieves the same results as the state-of-the-art method. This demonstrates the proficiency of our proposed efficient context aggregator to better capture the local and global contextual information for the CD task.

§.§ Qualitative Comparison

We present the qualitative results of ELGC-Net in Fig. <ref>, Fig. <ref>, and Fig. <ref> for LEVIR-CD, DSIFN-CD, and CDD-CD datasets, respectively. Fig. <ref> shows that

ELGC-Net is able to detect the small change (see first row) that is not detected by CNN-based FC-Siam-diff <cit.>, STANet <cit.>, DTCDSCN <cit.> as well as transformer-based BIT <cit.> and ChangeFormer <cit.> methods. Similarly, in the second and third rows, our method accurately detects the change
regions while ChangeFormer <cit.> and BIT <cit.> methods strive to identify the changes precisely.

In Fig. <ref>, our ELGC-Net is able to detect the true change regions while suppressing the false detections whereas   both CNN-based and transformer-based including  FC-Siam-diff <cit.>, STANet <cit.>, DTCDSCN <cit.>, BIT <cit.> and ChangeFormer <cit.> have more false change detections in addition to the true change regions. The reason is that our ELGC-Net has better capability to capture the discriminative features at local and global context levels and possesses the properties of both convolution and self-attention operations that complement each other for the detection of subtle and large change regions.

Fig. <ref> shows the performance comparison among  CNN-based  including FC-Siam-diff <cit.>, STANet <cit.>, DTCDSCN <cit.> and  transformer-based  such as BIT <cit.> and ChangeFormer <cit.>  with our ELGC-Net over the CDD-CD dataset. We see that the proposed ELGC-Net is capable of detecting subtle changes considerably better than the SOTA method (first row). Furthermore, the BIT <cit.> and ChangeFormer <cit.> approaches struggle to detect the change region accurately in the bottom-right of the second row which is precisely detected by our ELGC-Net.

§.§ Ablation Study

§.§.§ Ablation of contextual features in ELGCA module

We evaluate  the effectiveness of our contributions including  local contextual aggregation, multi-channel aggregation, and PT attention (global contextual aggregation) over the LEVIR-CD dataset. To do so, we iteratively apply different combinations of our contributions and present their parameters,

and IoU scores  as shown in Tab. <ref>.

Although the local contextual information using depthwise convolutions has minimum parameters (10.43M)

compared to other variants ( as illustrated in Tab. <ref> row 1), it has reduced  IoU score of 81.69%. Whereas PT attention has an IoU score of 82.55% while having increased parameters.

Nevertheless, our final network aggregating all the contextual information provides the best IoU score of 83.83%, with comparable  parameters

compared to row 1 as displayed in the last row in Tab. <ref>. It is evident that our final ELGCA module exhibits reduced parameters and better IoU score while leveraging the aggregation of different context aggregators in a parallel fashion.

§.§.§ Comparison of Parameters and FLOPs

We present the comparison of the parameters and FLOPs of our ELGC-Net with ChangeFormer <cit.> and baseline. Tab. <ref> shows that ELGC-Net has fewer parameters and FLOPs compared to state-of-the-art ChangeFormer <cit.> while it outperforms on all three datasets. We also demonstrate the effectiveness of our pooled-transpose (PT) attention
to have linear complexity with respect to tokens (
as shown in  Tab. <ref> row 3).
To do so, we replace the PT attention with standard self-attention  and observed that our PT attention exhibits fewer FLOPs while enhancing the IoU score. Furthermore,  the peak memory usage during the inference stage also shows that ELGC-Net consumes lesser memory.
Additionally, we present the results of ELGC-Net-LW that utilizes a light-weight decoder in Tab. <ref>. Although, the IoU score of ELGC-Net-LW is slightly reduced, however, the number of trainable parameters and FLOPs are considerably lower than the baseline and ChangeFormer which further demonstrates the effectiveness of the ELGCA module.

Comparison of different attention layers and ELGCA module on LEVIR-CD. The best results are shown in bold text.

Layers     Params (M)     FLOPs (M)     IoU     F1     OA

Self Attention (Baseline)     11.63     138557.9     81.93     90.07     99.00

Separable Attention <cit.>     11.32     124104.8     83.13     90.78     99.07

EfficientFormer Attention <cit.>     15.07     267518.1     81.64     89.89     98.99

ELGCA     10.57     123590.5     83.83     91.20     99.12

§.§.§ Comparison between other context aggregators and Our's ELGCA

In Tab. <ref>, we present the comparison of the ELGCA module with self-attention, separable attention <cit.>, and EfficientFormer attention <cit.> layer on LEVIR-CD. We replaced the ELGCA module in the proposed network with the self-attention, separable attention, and EfficientFormer attention layer, and train the network with the same hyperparameters. Tab. <ref> shows that self-attention and separable attention <cit.> have 11.63 and 11.32 million parameters respectively, and separable attention <cit.> performs better as compared to self-attention <cit.> and EfficientFormer attention <cit.> by achieving IoU score of 83.13%. However, the proposed ELGCA module performs remarkably better compared to these attention layers by obtaining an IoU score of 83.83% while having less number of trainable parameters and FLOPs.

Ablation study of different pooling layers in self-attention of ELGCA module on LEVIR-CD. Average and maximum pooling is applied to the query and key features of the self-attention before computing the matrix multiplication. An optimal combination of average and maximum pooling gives the best performance highlighted in bold.

2|c|Pooling applied to     2*IoU     2*F1     2*OA

Query Features     Key Features

-     -     82.80     90.59     99.05

Avg Pool     Avg Pool     83.30     90.89     99.08

Max Pool     Max Pool     83.00     90.71     99.06

Max Pool     Avg Pool     83.24     90.85     99.08

Avg Pool     Max Pool     83.83     91.20     99.12

Comparison of loss functions on LEVIR-CD applied during the training of our model. The best results are highlighted in bold text.

Loss Function     IoU     F1     OA

Focal Loss      79.27     88.44     98.84

mIoU Loss      77.51     87.32     98.75

Cross-Entropy Loss      83.83     91.20     99.12

§.§.§ Ablation of pooling layers in ELGCA

Fig. <ref>-(c) shows the architecture of the proposed ELGCA module. As discussed earlier in the ELGCA module,  First, average and maximum pooling operations are applied on the query and key features to obtain slightly robust features against minor variations. Secondly, we employ the transposed attention having linear complexity with respect to tokens, to reduce the computational cost of matrix multiplication. We conduct various experiments with different pooling layers and present the results in Tab. <ref>.
We notice that  comparable performance is achieved without applying pooling operations for key and query.
However, we notice from Tab. <ref> that the optimal combination of average and maximum pooling operations provides the best IoU score as 83.83% where we use average pooling operation for query feature maps and maximum pooling operation for key feature maps.

§.§.§  Ablation of losses applied during training of ELGCA

Finally, we perform an ablation study to show the impact of different losses utilized during the training of our proposed model.  From Tab. <ref>, it is notable that cross-entropy (CE) loss is more suitable to train our model compared to other losses.

§ CONCLUSION

In this work, we present an efficient global and local context aggregation network (ELGC-Net) for remote sensing change detection problem. ELGC-Net utilizes the local and global contextual representations in an efficient way to capture enhanced semantic information while reducing the number of trainable parameters.

Additionally, ELGC-Net does not require the use of a pre-trained backbone feature extractor. Extensive experimental results show the effectiveness of our approach. The proposed ELGC-Net outperforms on three challenging datasets and achieves state-of-the-art results. In addition, our lighter variant provides comparable performance while having reduced parameters and FLOPs.
The prospective future direction is to further increase the efficiency of the network and adapt the model for real-time change detection task on edge devices.


# Photonic Snake States in Two-Dimensional Frequency Combs

Salim B. Ivars  Institut Universitari de Matemàtica Pura i Aplicada,
Universitat Politècnica de València, 46022 València, Spain  Departament de
Física, Universitat Politècnica de Catalunya, 08222 Terrassa (Barcelona),
Spain  Yaroslav V. Kartashov  ICFO–Institut de Ciències Fotòniques, The
Barcelona Institute of Science and Technology, 08860 Castelldefels
(Barcelona), Spain  Institute of Spectroscopy, Russian Academy of Sciences,
Troitsk, Moscow, 108840, Russia  P. Fernández de Córdoba  Institut
Universitari de Matemàtica Pura i Aplicada, Universitat Politècnica de
València, 46022 València, Spain  J. Alberto Conejero  Institut Universitari de
Matemàtica Pura i Aplicada, Universitat Politècnica de València, 46022
València, Spain  Lluis Torner  ICFO–Institut de Ciències Fotòniques, The
Barcelona Institute of Science and Technology, 08860 Castelldefels
(Barcelona), Spain  Universitat Politècnica de Catalunya, 08034 Barcelona,
Spain  Carles Milián<EMAIL_ADDRESS>Institut Universitari de
Matemàtica Pura i Aplicada, Universitat Politècnica de València, 46022
València, Spain

###### Abstract

Taming the instabilities inherent to many nonlinear optical phenomena is of
paramount importance for modern photonics. In particular, the so-called snake
instability is universally known to severely distort localized wave stripes,
leading to the occurrence of transient, short-lived dynamical states that
eventually decay. The phenomenon is ubiquitous in nonlinear science, from
river meandering to superfluids, and to date it remains apparently
uncontrollable. However, here we show that optical snake instabilities can be
harnessed by a process that leads to the formation of stationary and robust
two-dimensional zigzag states. We find that such new type of nonlinear waves
exists in the hyperbolic regime of cylindrical micro-resonators and it
naturally corresponds to two-dimensional frequency combs featuring spectral
heterogeneity and intrinsic synchronization. We uncover the conditions of the
existence of such spatiotemporal photonic snakes and confirm their remarkable
robustness against perturbations. Our findings represent a new paradigm for
frequency comb generation, thus opening the door to a whole range of
applications in communications, metrology, and spectroscopy.

Figure 1:  | Cylindrical microresonator and multi-stable resonances. a, sketch of the driven (from the left) micro-cylinder in the photonic snake regime. b-g, maximum of the cavity background (single color) field amplitude, $|\psi_{0}(x,z)|$, versus laser-cavity detuning, $\delta$, showing selected nonlinear resonances for driving strengths b-d $h_{0}=1.4$, e-g $h_{0}=1.7$, and various values of the transverse pump width: b,e $\sigma_{z}=4$; c,f $\sigma_{z}=8$; d,g $\sigma_{z}=12$. The cavity background may be smooth [inset $(i)$] or may feature single- and multi-stripe quasi one-dimensional solitons distributed along $z$ [insets $(ii)-(viii)$]. The number of stripes on each portion of the resonances is colour coded in the legend of d. The resonance corresponding to the one-dimensional microring is also shown, for reference, by the dashed curves labeled as ’$1D$’ in b, e. Thick black dots on the resonances in b,c,e mark the onsets of instabilities and correspond to bifurcation points for spatiotemporal states, including snakes [see Fig.2]. Labels $SN_{1,2}$ in e mark the saddle node bifurcations (branch turning points) further discussed in Fig.3. Rectangles in f enclose the regions zoomed in Figs.4b and 4e. Background states are stable for $|\psi_{0}|<1$ and may be unstable otherwise, within the grey shaded areas. Insets showing spatial profiles are plotted over the area $x\in[-L/2,L/2[$ (covering the whole cavity length, $L=16$) and $z\in[-16,16]$.

Figure 2:  | Photonic snake families bifurcating from the single stripe background. a,b, amplitude vs detuning branches for all photonic snakes supported by the microcylinder considered in Fig.1 with $\sigma_{z}=4$ for a $h_{0}=1.4$ and b $h_{0}=1.7$. Thick (thin) traces denote stable (unstable) states. The nonlinear resonances (black curves) are those shown in Figs.1b and 1e, respectively. Photonic snakes exhibiting $N$ periods around the cavity circumference are denoted by $S_{N}$, and selected profiles, at the positions of the squares in b, are shown in insets $(i)-(viii)$. Panels $(ii),(iii)$ show real, imaginary parts of the snake $S_{1}$. Inset in b is a zoom over the region where the $S_{1}$ branch approaches the $S_{5}$ branch, and as a consequence, the corresponding snakes [inset ($ix$)] present a profile mixing both periodicities: 5 small zigzag periods within a 1 period envelope. Other branches exist (not shown), associated to $N$-dark soliton states, denoted by $D_{N}$, with typical profiles as shown in $(x)-(xii)$. c,d show the snake’s zigzagging angle [cf. inset ($iv$)] as a function of detuning, $\alpha(\delta)$, corresponding to all branches in a,b (branches are colour-matched). e shows the existence (light area) and stability (black area) domains in the $\\{\delta,\sigma_{z}\\}$ plane for the snake family $S_{2}$ with $h_{0}=1.7$. All insets are plotted over the area $x\in[-L/2,L/2[$ and $z\in[-8,8]$. Figure 3: | Onset of snakes. Growth rates vs detuning of the internal modes associated to the background states defining the cavity resonance in Fig. 1e [$\sigma_{z}=4$, $h_{0}=1.7$]. This figure is a parametric plot along the resonance, i.e., the abscissa corresponds to the values of detuning along the branches of the resonance, starting at $\delta=1.4$. The left panel ($\delta\in[1.4,SN_{1}\lesssim 3.5$) corresponds to the top branch, the central panel ($\delta\in[SN_{1},SN_{2}]$) corresponds to the middle branch, and the right panel ($\delta\in[SN_{2}\gtrsim 2.2,2.5]$) corresponds to the lower and stable branch. Labels $\chi_{\psi}$ denote the internal modes the growth of which leads to the state $\psi$ [$\psi$ are either snakes $S_{N}$ (black) or dark solitons $D_{N}$ (red)]. Insets at the top show the real part of typical snake (left) and neck (right) modes. Snake modes with $N=0$ induce, when excited, a global drift along $z$ while neck modes with $N=0$ correspond to the universally unstable mode that dotes the whole middle branch with the detrimental zero-wavelength instability. The zero growth rate points correspond to the bifurcation points shown in Figs. 1 and 2 from where the corresponding states (shown in Fig. 2) emerge. Insets are plotted over the area $x\in[-L/2,L/2[$ and $z\in[-8,8].$ Figure 4: Coupled photonic snakes. a, intensity of a 2-snake with 3 periods [denoted by $S_{3}^{(2)}$] for $h_{0}=1.7$ and $\sigma_{z}=8$. b, branches of its existence and stability shown by the solid and thick black curves, respectively. Branches with other colors correspond to the resonances in Fig. 1f. A square marker is around the location where the profile in a is taken ($\delta=1.72$). The stable branch emerges from the bifurcation point labeled with B (left most thick red dot). Inset: (top) 2-stripe state, (bottom) perturbation impingin the flexural instability that transforms the 2-stripe state into the $S_{3}^{(2)}$ in a. White arrows indicate the direction of the local drift [arrows in inset and in a are in exact correspondence]. c, Existence and stability chart in the $\\{\delta,\sigma_{z}\\}$-plane for $h_{0}=1.7$, built by analysing the stability of the stable branch in b and varying pump-width. Horizontal dashed line marks $\sigma_{z}=8$, corresponding to b. d-f are analogous to a-c for a triple snake. Other snaking states, albeit unstable: g,i show the internal modes of the (2,3)-stripe states leading to states in h,j. Additional bifurcations reshape states of the type j into k and l. All $(x,z)$ panels are plotted over the area $x\in[-L/2,L/2[$ and $z\in[-10,10].$ Figure 5: Two-dimensional heterogeneous combs. a-d, intensity, $|\psi(x,z)|^{2}$, for several snake profiles with different tilt, $\alpha$, periodicity, and detuning, $\delta$ (see labels) at $\sigma_{z}=4$ and $h_{0}=1.7$. e-h, Fourier transform along the $x$ direction of the fields in a-d, illustrating the two-dimensional comb structure vs $z$ around the central wave-numbers, $k_{x}$. The frequency combs exhibit a high degree of heterogeneity along the cylinder’s axis ($z$), as evident from insets (I)-(IV), showing the broad one-dimensional combs at different locations: $z=1.72$, 1.25, 0.63, 0, respectively. i-l, Two-dimensional Fourier transforms of the fields in a-d illustrating snakes in momentum space, which reveal the angular spread of the two-dimensional comb frequencies (see text). Figure axes: a-d, $x\in[-L/2,L/2[$ and $z\in[-5,5]$; e-h, $k_{x}\in[-20,20]$ and $z\in[-5,5]$; i-l, $k_{x}\in[-20,20]$ and $k_{z}\in[-20,20]$. Spectra in e-l are normalised to the strongest comb line excluding the pump at $k_{x}=0$.

Since their discovery almost half a century ago, the so-called flexural
instabilities zakh_rub , which lead to the unarrestable decay of quasi-one-
dimensional self-sustained states such as solitons zakhPR ; kivpel , have been
encountered in all areas of nonlinear sciences. Namely, they have been found
to occur in the dynamics of river meandering meander , classical fluids fluid
, fermion yefNat ; cetoli  and polariton claude  superfluids, Bose-Einstein
condensates bec , chemistry chem , and, importantly, nonlinear photonics mam ;
tik ; gorza04 ; gorza11 ; gorza06 . The initial stages of the instability $-$
referred to as the snake instability zakhPR $-$ evinces the reshaping of
straight states into zigzags caused by spontaneous and local transverse drifts
with alternating directions, which occur prior to the irreversible decay of
the corresponding states. In contrast to the use of the instability in several
systems, e.g., to trigger spontaneous vortex formation bec ; mam ; tik , its
controlled arrest to form robust multi-dimensional snaking states remains
hitherto unexpected.

In parallel, in recent years, microring cavities have emerged as an
outstanding platform for exploring the fundamental properties and the
applications of stabilised one-dimensional dissipative nonlinear waves, in
both the normal and the anomalous dispersion regimes, under a plethora of
potentially unfavourable situations, such as strong dispersive braschSci ;
vahNatCom , inelastic karpov ; VahalaStokes ; gaeta , and thermal stone ;
obrther  effects. Each and every new solitonic states found in such micron-
sized systems correspond to highly stable frequency combs that contribute to
the developing field of comb-related applications pasqPR ; cundiffRMP . As a
consequence, microring-like cavities are highly appealing systems to push the
frontiers of knowledge of nonlinear phenomena as well as boosting applications
in photonics. Nowadays, most solitonic states (stable combs) that have been
experimentally observed in microrings occur in essentially single-comb forming
devices, except for the few exceptions that may be found in counter-
propagating soliton experiments vahCount ; kippCount , in systems with two
microrings dimerKipp ; gaetfew , or in cavities with a few transverse modes
kippsdm , where each of the individual combs are formally one-dimensional. The
potential of two-dimensional frequency combs, which to date is essentially
unexplored experimentally, remains to be uncovered.

Here we discover the existence of robust spatiotemporal photonic snakes in
cylindrical micro-resonators with normal group velocity dispersion (GVD) and
in the presence of diffraction along the cylinder’s axis. Such previously
unknown states are a continuous two-dimensional ensemble of heterogeneous
combs which are inherently synchronised by the nonlinearity. Thus, they
represent a whole new paradigm for frequency comb formation.

Results

Model. The spatiotemporal snakes predicted here form in hollow cylindrical
Kerr micro-resonators [cf. Fig. 1a]. The dynamics of the intra-cavity field
can be accurately described by the two-dimensional generalisation of the
Lugiato-Lefever equation lugiato87 ,

$\partial_{t}\psi=\frac{i}{2}(-\partial_{x}^{2}+\partial_{z}^{2})\psi-(1+i\delta)\psi+i|\psi|^{2}\psi+ih_{0}e^{-z^{2}/\sigma_{z}^{2}},$
(1)

where the terms on the right-hand side of the equation account for normal GVD
along the propagation (and periodic) coordinate, diffraction along the
cylinder’s axis, losses, laser-cavity detuning, focusing nonlinearity, and the
axially localised pump, respectively (see methods). For the sake of
generality, we discuss photonic snakes first in the physical frame of Eq.1
[cf. Figs. 1-5], which represents a universal hyperbolic model governing a
plethora of nonlinear physical phenomena in optics op1 ; op2  and other areas
of physics kevrev . In order to anchor our predictions into realistic optical
frequency comb forming experimental realizations, we report on the snakes’
robustness under the specific collection of perturbations expected in
cylindrical micro-cavities [cf. Fig. 6], including Raman and thermal
nonlinearities, among others (see methods).

Nonlinear resonances. When driven into the strongly nonlinear regime, the here
considered micro-cylinders develop nonlinear resonances, associated to the
single colour locked states, $\psi_{0}$
($\partial_{t}\psi_{0}=\partial_{x}\psi_{0}=0$), exhibiting an unusually rich
multi-stability landscape, as shown in Figs.1b-g. For peak intensities
$|\psi_{0}|\gtrsim 1$, corresponding to the existence domain of bright
solitons with anomalous GVD BarashPRE96  (formally analogous to diffraction
along $z$), the background states become single and multi-stripe solitons [cf.
insets $(ii)$-$(viii)$] distributed on top of the smooth background [cf. inset
$(i)$], recalling the spatial solitons reported in pioneering works mclaPRL83
. The presence of stripes impinges an intricate morphology of the resonances
strongly deviating from that of the tri-valued (bistable) one-dimensional
counterparts (marked for reference by the light-grey curves in Figs.1b-g).
Resonances associated to wider or more intense pumps (i.e., greater
$\sigma_{z}$ and $h_{0}$, respectively) contain states hosting more stripes
and reshape to exhibit splitting (Figs. 1c,e), nesting (Figs. 1d,f,g), and
closed loops (Fig. 1f). As discussed below, the properties of resonances are
intimately linked to the existence of snake states.

Photonic snakes. The soliton stripes described above are prone to snake
instabilities. When these develop, stripes distort and acquire in our
dissipative system a periodic zigzag snaking profile, becoming inhomogeneous
along $x$ and, hence, polychromatic. The central result of this work is that
the transverse drifts induced by the snake instability are fully arrested
after they exerted certain distortion to stripes so that perfectly stationary
snakes form [see Figs. 2$(i)$-$(ix)$ for typical profiles]. Figures 2a, b
present the branches [$\max(|\psi(x,z)|)$ vs detuning] corresponding to the
complete set of snakes existing with narrow pump, $\sigma_{z}=4$, for
different driving strengths, $h_{0}$ (see caption). Snakes featuring from 1 to
5 periods in the microcavity circumference exist and those with 2 to 4 periods
are stable for some detuning intervals. We denote different branches of snakes
featuring $N$ periods on the microresonator circumference as $S_{N}$. Figures
2c,d show that snakes feature pronounced zigzagging angles, $\alpha$ [cf.
inset $(iv)$], over most of their existence regions, so that they have a
strong trend to form narrow pulses along $x$. We emphasize that the robust
snakes exist over a finite and generous region of the parameter space, as
shown in Figs. 2a-d. To further illustrate this, Fig. 2e shows the existence
domain on the $\\{\delta,\sigma_{z}\\}$ plane for the $S_{2}$ family, which is
found to be stable within the black area, unstable within the light shaded
area, and non-existent otherwise (white area). In this work, optimal
conditions for the snakes stability are encountered only for strong pump
localisation along $z$. Outside the stability domain, snakes are exposed to
oscillatory and exponential instabilities leading, respectively, to breathing
and decay. The precise combination of the driving amplitude and aspect ratio
of the pump field, $\sigma_{z}/L$, (lying in the range
$[\frac{1}{8},\frac{1}{3}]$ in Fig.2e) is essential for the formation of
robust snakes.

Photonic snake families emerge supercritically from the top branch of the
cavity resonance, each with a different detuning threshold (marked with black
dots in Fig.2), corresponding to the detuning values at which stripes become
unstable due to the growth of snake type perturbations or internal modes
internal2  (see methods for the stability analysis description). Figure 3
shows all relevant growth rates for the three branches (top, middle, and
bottom) of the right-most cavity resonance shown in Fig. 1e. Snake-type
perturbations ($\chi_{S}$, see top insets in Fig. 3) with different periods
acquire positive growth rates on the top and middle branches of the resonance.
The nullity of growth rates defines exactly the bifurcation points (black dots
in Fig. 2b) or loci from which snakes emerge from stripes. We note that the
symmetry properties of the stripes internal modes allows one to readily
predict the symmetries to be inherited by the emerging snake family.

Although our attention is focused on the snake instability, we point out that
the middle branch of the resonance in Fig. 3 also presents neck instabilities
zakh_rub , which often compete with the snake type skrPRL98 and manifests upon
the growth of the corresponding internal modes, denoted by
$\chi_{D_{1}}$-$\chi_{D_{3}}$ and shown in the right inset of Fig. 3. In this
work, all of the spatiotemporal states found with traces of neck
perturbations, namely, the dark solitons in Figs. 2($x$-$xii$) and the hybrid
snake-dark solitons in Fig. 4j,l, were highly unstable and thus their
properties are not discussed in details.

The photonic snake families become much richer when snake instabilities
develop on the multi-stripe states present on cavity resonances with large
pump beam widths, $\sigma_{z}\gtrsim 8$ (cf. Fig. 1), what leads to the
formation of coupled multi-snake states. Figure 4a shows an example of robust
two-snake state with three periods along the microcylinder’s circunference,
denoted by $S_{3}^{(2)}$. This double snake forms when the snake-type mode
$\chi_{S_{3}^{(2)}}$ (Fig. 4b bottom-right inset) grows on top of the two-
stripe state (Fig. 4b top-right inset). The full amplitude vs detuning branch
of the $S_{3}^{(2)}$ snake is shown in Fig. 4b by the black curve. Other
branches correspond to the cavity resonance (cf. Fig. 1f). Existence and
stability of the $S_{3}^{(2)}$ family with $h_{0}=1.7$ is shown vs pump width
and detuning over the relevant parameter space in Fig. 4c. Figures 4d-f show
analogous results for a triple-snake. The above two snake families correspond
to particular cases where their branches connect background states with
different number of stripes: the $S_{3}^{(2)}$ [$S_{3}^{(3)}$] family
bifurcates at lower $\delta$ from the $\mathrm{B}$-point located on a 2 (3)
stripe state and merges with 4-stripe state at a higher $\delta$ (red dots).

Even though we only found stable in-phase multi-snakes, the system supports a
vast collection of states with different morphology. Some of these are
illustrated in Fig. 4h and Figs. 4j,k,l which appear after the anti-phase
snake modes in Figs. 4g and 4i grow, respectively, on two- and three-stripe
states. This results in the formation of anti-phase snakes that appear alone
(Fig. 4h), mixed with dark solitons (Fig. 4j), together with ellipsoids (Fig.
4k), or displaying asymmetries along $z$ (Fig. 4l).

Heterogeneous 2D combs. The most remarkable feature of photonic snake states
is that their spectrum is heterogeneous along $z$, while being inherently
synchronized by the nonlinearity. Figure 5 shows a soliton stripe (Fig. 5a)
and snakes of different tilts and periods (Figs. 5b-d), together with their
frequency comb distribution along $z$ (Figs. 5e-h). Heterogeneity is evident
in Figs. 5g,h, showing Fourier spectra along the $z$-axis, and it is further
illustrated in insets $I-IV$ by showing different combs extracted at different
axial positions (marked by horizontal white lines in Figs. 5c,g).
Synchronisation and heterogeneity of these spectra is readily important for
metrology and spectroscopy picqueNP . The spreading of the combs along the
cylinder’s axis, $z$, naturally introduces the notion of heterogeneous two-
dimensional comb, which constitutes a generalisation of the widely reported
one-dimensional comb and a central result of this work.

Figures 5i-l show the snakes in Figs. 5a-d, respectively, in the two-
dimensional momentum space, illustrating the angular spread of the different
spectral components, potentially important for an efficient collection of the
combs by external tapers or waveguides. Because dispersion in cylinders is
typically much smaller in $x$ than in $z$ gorodjosab ; milianPRL18 , the
physical values corresponding to $k_{z}$ are in practice much smaller than
those corresponding to $k_{x}$, so that the propagation angles
$\theta=\arctan({k_{z}/k_{x}})$ remain in the order of a few degrees, at most
(see methods). On the same reason, the height of panels $e-h$ is of the order
of a few mm’s, so that the comb heterogeneity occurs along $z$ over larger
scale than typical taper fiber widths, which enables the efficient light
collection by, e.g., arrays of waveguides.

Figure 6: Effect of perturbations. a, cavity resonance and the branch for the
two-period snakes, $S_{2}$, affected by perturbation originating from higher
order dispersion, Raman scattering and thermal nonlinearity (black) [resonance
is solely affected by thermal effects, and hence the label hot]. The results
without perturbations are also shown for reference (gray, cf. Fig.2a). b-d
show, respectively, the light field, the Raman vibration, and the thermal
distribution for the stable snake at $\delta=2$. e, the corresponding two-
dimensional comb. f, the maximum amplitude of the light field during the
propagation of the stable and robust snake ($\delta=2$), for a total time of
$155$ ns [50,000 roundtrips], under initial noise and perturbations at each
step. Right-inset: complex eigenvalues, $\lambda$, obtained from the stability
analysis performed on the stationary snake used as input of the propagation,
taking into account all above mentioned perturbations. Instabilities
correspond to positive real parts of $\lambda$. g, analogous results as for f,
but for an oscillatory unstable snake ($\delta=2.24$). Small insets illustrate
the dynamical reshaping of the snake with time. Panels b-d and insets in g
span over the whole cavity circumference and over 4 mm along $z$.

Robustness. The photonic snakes and the two-dimensional combs exist in a
physical setting that is readily realizable. In particular, consider a hollow
micro-cylinder made of silica glass, with a radius $R=100\ \mu$m and wall
thickness $w=0.75\ \mu$m, pumped at $\lambda_{p}\approx 1.26\ \mu$m, with a
$\mathcal{Q}$-factor of $\sim 5\times 10^{6}$, and take into account the
expected specific effects introduced by chromatic dispersion, Raman and
thermal nonlinearities, and the localised coupling, in $x$ and $z$, between
the pump beam and the micro-cavity (see the full model in methods). Although
each of the above perturbations may be suppressed in various ways, we now show
that robust snakes occur in passive driven cylindrical micro-resonators.
Figure 6 presents an overview of our predictions in such scenarios, focusing
the attention on the snake family $S_{2}$ with $h_{0}=1.4$ and $\sigma_{z}=4$.
Figure 6a shows the thermal shift of the nonlinear resonance [black] with
respect to the unperturbed cold resonance [grey], as well as the $S_{2}$ snake
branch, which stable region is highlighted by the thick trace. An example of
robust snake is shown in Fig.6b, along with its associated Raman vibration
(Fig.6c) and thermal (Fig.6d) fields. The corresponding two-dimensional comb
is shown around the central lines in Fig.6e. We emphasize that in the context
of one-dimensional combs, solitons and related nonlinear waves have been
reported under the thermal nonlinearity iltherm , both experimentally stone
and theoretically menyukther . In one dimension, thermal effects mainly induce
a global shift in the laser-cavity detuning which, even if they introduce
their own instabilities iltherm ; menyukther , cannot compromise soliton
existence itself. However, in the two-dimensional case here considered,
thermal effects are inhomogenous along $z$ (cf. Fig.6d), and hence induce
thermal stress along $z$. As a consequence, the mere existence of the two-
dimensional snakes cannot be anticipated a priori, and neither can their
stability. The latter is formally predicted via linear stability analysis and
explicitly checked by long propagation runs (cf. Fig.6f) spanning over
$50,000$ cavity roundtrips (equivalent to 155 ns in our geometry). Propagation
of robust snakes under the whole plethora of realistic effects feature
constant peak amplitudes (as a consequence of the steady two-dimensional
profiles) despite the input random noise and step-to-step perturbations. On
the contrary, unstable snakes behave very differently. In particular, Fig.6g
shows the time evolution of an oscillatory unstable snake, featuring peak
amplitude and profile variations (see insets). The presence and absence of
instabilities, and the nature of them in the former case, are accurately
predicted by the stability analysis (see methods), which results are displayed
in insets of Figs.6f,g.

Discussion.

Stable photonic snakes present a trend to bifurcate super-critically from the
stable background upon the increase of detuning [cf. Figs. 2,4], even in the
presence of higher order effects [cf. Fig. 6], what strongly suggests that
they will be easily and deterministicaly excitable in realistic experiments by
the standard dynamical red-shift of the pump’s frequency. This simple
excitation mechanism of frequency combs in the normal GVD regime of passively
driven cavities, naturally provided here by the intrinsically two-dimensional
snakes, has remained hidden till date by the low dimensionality of microrings,
where relatively complex excitation methods were typically required. Indeed,
in one-dimensional microrings, solitonic combs in the normal GVD regime are
often excited by engineering a small anomalous GVD region around the pump’s
frequency, so modulational instability triggers comb formation. This is
achieved via mode coupling effects arising due to avoided crossings in multi-
modal microrings saavOE12 ; XueNP ; jang16 ; nazem21  or with the aid of an
auxiliary ring xuebis ; kim . A recent work demonstrates a much simpler and
straight-forward mechanism to excite (one-dimensional) dark solitons via self-
injection-locking VahalaSILDarlSol  and displaying the turnkey turnkey
operation. Photonic snake states may thus belong to the collection of simple
and deterministically excitable combs in cavities with dominant normal GVD,
which is of central importance for extending the formation of micro-cavity
combs to the long and short wavelength regions, far away from the telecom band
around $1.5\ \mu$m.

In addition, the two-dimensional combs reported here automatically enable the
possibility to host synchronised heterogeneous combs in a single device,
features that are of great importance in metrology and spectroscopy, and may
be difficult to combine in this specific manner with microrings. Indeed, comb
heterogeneity may be achieved in the one-dimensional context, e.g., via the
excitation of unbounded solitons copropagating along the same spatial channel
heteroNC . However, their different group velocities unavoidably yield de-
synchronisation. On the other hand, synchronisation was previously reported
for identical combs in distant microrings gaetSync . Remarkably, both features
were found simultaneously in the bi-modal stokes solitons VahalaStokes , where
heterogeneity appeared as a result of two copropagating and spectrally non-
overlapping combs. In the present work, differently with the above, and akin
to the two-dimensional geometry, the heterogeneity appears as a continuous
comb reshaping along the cylinder’s axis and with a fixed carrier frequency.
We note that the above remarkable findings heteroNC ; gaetSync ; VahalaStokes
were reported in the anomalous GVD region, where the physics is substantially
different than that in the normal GVD regime, subject of this work.

In summary, we have uncovered a fundamentally new mechanism to arrest and
control the ubiquitous and otherwise strong snake instability in cylindrical
Kerr micro-resonators and, as a result, the possibility to form complex robust
spatiotemporal snake states. Specifically, we have found that photonic snakes
naturally form perfectly synchronised heterogeneous comb ensembles in the
normal GVD (hyperbolic) regime of cylindrical micro-cavities. The phenomenon
represents a novel paradigm for the generation of optical frequency combs.

Online content

Methods, additional references, statements of data availability,
acknowledgements, details of author contributions and competing interests are
available in the online version of the paper.

## References

  * (1) V. E. Zakharov, A. M. Rubenchik, “Instability of waveguides and solitons in nonlinear media,” Sov. Phys. JETP 38, 494–500 (1974).
  * (2) E. A. Kuznetsov, A. M. Rubenchik, V. E. Zakharov,“Soliton stability in plasmas and hydrodynamics,” Phys. Rep. 142, 103–65 (1986).
  * (3) Y. S. Kivshar and D. E. Pelinovsky, “Self-focusing and transverse instabilities of solitary waves,” Phys. Rep. 331, 117-195 (2000).
  * (4) J. A. Constantin, T. Dunne, J. Ahmed, C. Legleiter, and E. D. Lazarus, “Sediment supply as a driver of river meandering and floodplain evolution in the Amazon Basin,” Nature Geoscience 7, 899-903 (2014).
  * (5) E. D. Brown, S. B. Buchsbaum, R. E. Hall, J. P. Penhune, K. F. Schmitt, K. M. Watson, and D. C. Wyatt, “Observations of a nonlinear solitary wave packet in the Kelvin wake of a ship,” J. Fluid Mech. 204, 263-293 (1989).
  * (6) T. Yefsah, A. T. Sommer, M. J. Ku, L. W. Cheuk, W. Ji, W. S. Bakr, and M. W. Zwierlein, “Heavy solitons in a fermionic superfluid,” Nature 499, 426-430 (2013).
  * (7) A. Cetoli, J. Brand, R. G. Scott, F. Dalfovo, and L. P. Pitaevskii, “Snake instability of dark solitons in fermionic superfluids,” Phys. Rev. A 88, 043639 (2013).
  * (8) F. Claude, S. V. Koniakhin, A. Maître, S. Pigeon, G. Lerario, D. D. Stupin, Q. Glorieux, E. Giacobino, D. Solnyshkov, G. Malpuech, and A. Bramati, “Taming the snake instabilities in a polariton superfluid,” Optica 7, 1660-1665 (2020).
  * (9) B. P. Anderson, P. C. Haljan, C. A. Regal, D. L. Feder, L. A. Collins, C. W. Clark, and E. A. Cornell, “Watching dark solitons decay into vortex rings in a Bose-Einstein condensate,” Phys. Rev. Lett. 86, 2926-2929 (2001).
  * (10) C. Luengviriya, U. Storb, G. Lindner, S. C. Müller, M. Bär, and M. J. B. Hauser, “Scroll Wave Instabilities in an Excitable Chemical Medium,” Phys. Rev. Lett. 100, 148302 (2008).
  * (11) A.V. Mamaev, M. Saffman, and A. A. Zozulya, “Propagation of dark stripe beams in nonlinear media: snake instability and creation of optical vortices,” Phys. Rev. Lett. 76, 2262-2265 (1996).
  * (12) V. Tikhonenko, J. Christou, B. Luther-Davies, and Y. S. Kivshar, “Observation of vortex solitons created by the instability of dark soliton stripes,” Opt. Lett. 21, 1129-1131 (1996).
  * (13) S.-P. Gorza, N. Roig, Ph. Emplit, and M. Haelterman, “Snake Instability of a Spatiotemporal Bright Soliton Stripe,” Phys. Rev. Lett. 92, 084101 (2004).
  * (14) S.-P. Gorza, Ph. Emplit, and M. Haelterman, “Observation of the snake instability of a spatially extended temporal bright soliton,” Opt. Lett. 31, 1280-1282 (2006).
  * (15) S.-P. Gorza, B. Deconinck, Ph. Emplit, T. Trogdon, and M. Haelterman, “Experimental Demonstration of the Oscillatory Snake Instability of the Bright Soliton of the $(2+1)D$ Hyperbolic Nonlinear Schrödinger Equation,” Phys. Rev. Lett. 106, 094101 (2011).
  * (16) V. Brasch, M. Geiselmann, T. Herr, G. Lihachev, M. H. P. Pfeiffer, M. L. Gorodetsky, and T. J. Kippenberg, “Photonic chip–based optical frequency comb using soliton Cherenkov radiation,” Science 351, 357–360 (2016).
  * (17) X. Yi, Q.-F. Yang, X. Zhang, K. Y. Yang, X. Li, and K. Vahala, “Single-mode dispersive waves and soliton microcomb dynamics,” Nat. Commun. 8, 1-9 (2017).
  * (18) M. Karpov, H. Guo, A. Kordts, V. Brasch, M. H. P. Pfeiffer, M. Zervas, M. Geiselmann, and T. J. Kippenberg, “Raman self-frequency shift of dissipative Kerr solitons in an optical microresonator,” Phys. Rev. Lett. 116, 103902 (2016).
  * (19) Q.-F. Yang, X. Yi, K. Y. Yang, and K. Vahala, “Stokes solitons in optical microcavities,” Nat. Phys. 13, 53–57 (2017).
  * (20) M. Yu, Y. Okawachi, R. Cheng, C. Wang, M. Zhang, A. L. Gaeta, and M. Lončar, “Raman lasing and soliton mode-locking in lithium niobate microresonators,” Light: Science & Applications 9, 1-7 (2020).
  * (21) J. R. Stone, T. C. Briles, T. E. Drake, D. T. Spencer, D. R. Carlson, S. A. Diddams, and S. B. Papp, “Thermal and nonlinear dissipative-soliton dynamics in Kerr-microresonator frequency combs,” Phys. Rev. Lett. 121, 063902 (2018).
  * (22) E. Obrzud, S. Lacomte, and T. Herr, “Temporal solitons in microresonators driven by optical pulses,” Nat. Photon. 11, 600-607 (2017).
  * (23) A. Pasquazi, M. Peccianti, L. Razzari, D. J. Moss, S. Coen, M. Erkintalo, Y. K. Chembo, T. Hansson, S. Wabnitz, P. Del’Haye, X. Xue, A. M. Weiner, and R. Morandotti, “Micro-combs: a novel generation of optical sources,” Phys. Rep. 729, 1–81 (2017).
  * (24) S. T. Cundiff and J. Ye, “Colloquium: Femtosecond optical frequency combs,” Rev. Mod. Phys. 75, 325-342 (2003).
  * (25) Q. F. Yang, X. Yi, K. Y. Yang, and K. Vahala, “Counter-propagating solitons in microresonators,” Nat. Photon. 11, 560-564 (2017).
  * (26) W. Weng, R. Bouchand, E. Lucas, and T. J. Kippenberg, “Polychromatic Cherenkov radiation induced group velocity symmetry breaking in counterpropagating dissipative Kerr solitons,” Phys. Rev. Lett. 123, 253902 (2019).
  * (27) A. Tikan, J. Riemensberger, K. Komagata, S. Hönl, M. Churaev, C. Skehan, H. Guo, R. N. Wang, J. Liu, P. Seidler, and T. J. Kippenberg, “Emergent nonlinear phenomena in a driven dissipative photonic dimer,” Nat. Phys. 17, 604–610 (2021).
  * (28) J. K. Jang, X. Ji, C. Joshi, Y. Okawachi, M. Lipson, and A. L. Gaeta, “Observation of Arnold Tongues in Coupled Soliton Kerr Frequency Combs,” Phys. Rev. Lett. 123, 153901 (2019).
  * (29) E. Lucas, G. Lihachev, R. Bouchand, N. G. Pavlov, A. S. Raja, M. Karpov, M. L. Gorodetsky, and T. J. Kippenberg,“Spatial multiplexing of soliton microcombs,” Nat. Photon. 12, 699–705 (2018).
  * (30) L. A. Lugiato and R. Lefever, “Spatial dissipative structures in passive optical systems,” Phys. Rev. Lett. 58, 2209-2211 (1987).
  * (31) C. Conti, S. Trillo, P. Di Trapani, G. Valiulis, A. Piskarskas, O. Jedrkiewicz, and J. Trull, “Nonlinear electromagnetic X waves,” Phys. Rev. Lett. 90, 170406 (2003).
  * (32) P. Di Trapani, G. Valiulis, A. Piskarskas, O. Jedrkiewicz, J. Trull, C. Conti, and S. Trillo, “Spontaneously generated X-shaped light bullets,” Phys. Rev. Lett. 91, 093904 (2003).
  * (33) L. A. Cisneros-Ake, R. Carretero-González, P. B. Kevrekidis, and B. A. Malomed, “Dynamics and stabilization of bright soliton stripes in the hyperbolic-dispersion nonlinear Schrödinger equation,” Commun. Nonlinear Sci. Numer. Simul. 74, 268-281 (2019).
  * (34) I. Barashenkov and Y. S. Smirnov, “Existence and stability chart for the ac-driven, damped nonlinear Schrödinger solitons,” Phys. Rev. E 54, 5707-5725 (1996).
  * (35) D. W. Mc Laughlin, J. V. Moloney, and A. C. Newell, “Solitary Waves as Fixed Points of Infinite-Dimensional Maps in an Optical Bistable Ring Cavity,” Phys. Rev. Lett. 51, 75–78 (1983).
  * (36) D. E. Pelinovsky, Y. S. Kivshar, and V. V. Afanasjev, “Internal modes of envelope solitons,” Physica D 116, 121-142 (1998).
  * (37) D. V. Skryabin and W. J. Firth, “Modulational instability of solitary waves in non-degenerate three wave mixing: the role of phase symmetries,” Phys. Rev. Lett. 81, 3379-3382 (1998).
  * (38) N. Picqué and T. W. Hänsch, “Frequency comb spectroscopy,” Nat. Photon. 13, 146–157 (2019).
  * (39) Y. A. Demchenko and M. L. Gorodetsky, “Analytical estimayes of eigenfrequencies, dispersion, and field distribution in whispering gallery resonators,” J. Opt. Soc. Am. B 30, 3056–3063 (2013).
  * (40) C. Milián, Y. V. Kartashov, D. V. Skryabin, and L. Torner, “Clusters of cavity solitons bounded by conical radiation,” Phys. Rev. Lett. 121, 103903 (2018).
  * (41) V. S. Ilchenko and M. L. Gorodetskii, “Thermal nonlinear effects in optical whispering gallery microresonators,” Laser Phys. 2, 1004-1009 (1992).
  * (42) A. Leshem, Z. Qi, T. F. Carruthers, C. R. Menyuk, and O. Gat, “Thermal instabilities, frequency-comb formation, and temporal oscillations in Kerr microresonators,” Phys. Rev. A 103, 013512 (2021).
  * (43) A. A. Savchenkov, A. B. Matsko, W. Liang, V. S. Ilchenko, D. Seidel, and L. Maleki, “Kerr frequency comb generation in overmoded resonators,” Opt. Express 20, 27290-27298 (2012).
  * (44) X. Xue, Y. Xuan, Y. Liu, P.-H. Wang, S. Chen, J. Wang, D. E. Leaird, M. Qi, and A. M. Weiner, “Mode-locked dark pulse Kerr combs in normal-dispersion microresonators,” Nat. Photon. 9, 594-600 (2015).
  * (45) J. K. Jang, Y. Okawachi, M. Yu, K. Luke, X. Ji, M. Lipson, and A. L. Gaeta, “Dynamics of mode-coupling-induced microresonator frequency combs in normal dispersion,” Opt. Express 24, 28794-28803 (2016).
  * (46) E. Nazemosadat, A. Fülöp, Ó. B. Helgason, P.-H. Wang, Y. Xuan, D. E. Leaird, M. Qi, E. Silvestre, A. M. Weiner, and V. Torres-Company, “Switching dynamics of dark-pulse Kerr frequency comb states in optical microresonators,” Phys. Rev. A 103, 013513 (2021).
  * (47) X. Xue, Y. Xuan, P.-H. Wang, Y. Liu, D. E. Leaird, M. Qi, and A. M. Weiner, “Normal-dispersion microcombs enabled by controllable mode interactions,” Laser Photon. Rev. 9, L23-L28 (2015).
  * (48) B. Y. Kim, Y. Okawachi, J. K. Jang, M. Yu, X. Ji, Y. Zhao, C. Joshi, M. Lipson, and A. L. Gaeta, “Turn-key, high-efciency Kerr comb source,” Opt. Lett. 44, 4475–4478 (2019).
  * (49) W. Jin, Q.-F. Yang, L. Chang, B. Shen, H. Wang, M. A. Leal, L. Wu, M. Gao, A. Feshali, M. Paniccia, K. J. Vahala, J. E. Bowers, “Hertz-linewidth semiconductor lasers using CMOS-ready ultra-high-Q microresonators,” Nat. Photon. 15, 346-353 (2021).
  * (50) B. Shen, L. Chang, J. Liu, H. Wang, Q. F. Yang, C. Xiang, R. N. Wang, J. He, T. Liu, W. Xie, J. Guo, D. Kinghorn, L. Wu, Q. Ji, T. J. Kippenberg, K. Vahala, and J. E. Bowers, “Integrated turnkey soliton microcombs,” Nature 582, 365-369 (2020).
  * (51) W. Weng, R. Bouchand, E. Lucas, E. Obrzud, T. Herr, and T. J. Kippenberg, “Heteronuclear soliton molecules in optical microresonators,” Nat. Commun. 11, 1-9 (2020).
  * (52) J. K. Jang, A. Klenner, X. Ji, Y. Okawachi, M. Lipson, and A. L. Gaeta, “Synchronization of coupled optical microresonators,” Nat. Photon. 12, 688–693 (2018).

Online Methods

Time evolution equations. The nonlinear dynamics of the electric field’s
envelope orbiting around a cylindrical microresonator is described by the
following system of coupled equations,

$\displaystyle
i\partial_{t}\psi=-ib_{1,1}\partial_{X}\psi+\frac{1}{2}(\partial_{X}^{2}-\partial_{z}^{2})\psi+\hat{D}_{hod}\psi-[i-\delta]\psi-$
$\displaystyle-[(1-f_{R})|\psi|^{2}-\Delta_{T}+Q]\psi-
h_{0}\exp\left(-\frac{z^{2}}{\sigma_{z}^{2}}\right)\xi(X),$ (2)
$\displaystyle\partial_{t}\Delta_{T}=-A\int_{0}^{L}|\psi|^{2}\frac{dX}{L}-B\Delta_{T},\
\ $ (3)
$\displaystyle\partial_{t}^{2}Q=-\frac{2\gamma_{R}\tau}{\gamma}\partial_{t}Q-\frac{\tau^{2}\Omega_{R}^{2}}{\gamma^{2}}\left[Q-f_{R}|\psi|^{2}\right],\
\ $ (4)
$\displaystyle\xi(X)=\frac{1}{N}\sum_{m=-\infty}^{+\infty}\exp\left(-\frac{[X+mL]^{2}}{\sigma_{X}^{2}}\right),$
(5)

where $\psi$, $\Delta_{T}$, $Q$ are the optical, thermal, and molecular
vibrational fields, respectively. The above model assumes that only one radial
mode family of the cylinder is at play, which may be unambiguously achieved by
considering a hollow cylinder with an thin wall width, e.g., $w=0.75\ \mu$m,
as we used in Fig.6. The cylinder’s dispersion is given by

$\displaystyle\hat{D}\equiv\sum_{q=0}^{+\infty}\sum_{p=q}^{+\infty}b_{q,p-q}(-i\partial_{z})^{p-q}(-i\partial_{X})^{q},$
(6) $\displaystyle
b_{q,p-q}\equiv\frac{B_{q,p-q}\gamma^{p/2-1}}{2^{p/2}|B_{0,2}|^{(p-q)/2}|B_{2,0}|^{q/2}},$
(7) $\displaystyle B_{q,p-q}\equiv\frac{\tau\omega^{(q,p-q)}}{q!(p-q)!(2\pi
R)^{p}},$ (8)
$\displaystyle\omega^{(q,p-q)}\equiv\left.\partial_{k_{x}}^{q}\partial_{k_{z}}^{p-q}\omega(k_{x},k_{z})\right|_{k_{x0},k_{z0}},$
(9)

where $\gamma$ is the normalised cavity loss, $R$ is the cylinder’s radius,
$\tau$ the roundtrip time, and $k_{x}$, $k_{z}$ the wavenumbers associated to
the the $X$, $z$ coordinates ($X$ is the frame at rest in the lab). The
dispersion terms with low $p,q$ indices account for: $p=0,q=0$, resonance
frequency ($\omega_{0}$); $p=1,q=1$, group velocity of the pump’s frequency
along $X$ ($b_{1,1}$); $p=2,q=2$, GVD; $p=2,q=0$, diffraction. The rest of
terms are all included in the higher order dispersion operator
$\hat{D}_{hod}\equiv\hat{D}-\omega_{0}\tau/\gamma+ib_{1,1}\partial_{X}-\frac{1}{2}(\partial_{X}^{2}-\partial_{z}^{2})$.
The even parity of the cylinder’s dispersion around $z$ (see, e.g., Gorjosab ;
joannoPRL ) and the fact that we expand around the $k_{z0}=0$ yields the
nullity of all coefficients with $p-q=1$. The relation between normalised and
physical coordinates is as follows: $X=X_{phys}/(2\pi
R)\sqrt{\gamma/(2|B_{2,0}|)}$, $z=Z_{phys}/(2\pi
R)\sqrt{\gamma/(2|B_{0,2}|)}$.

In our simulations, the width of the numerical window along $x$ was $L=16$,
which together with the choice of $\gamma=0.001$ sets $B_{2,0}\approx-2\times
10^{-6}$, attainable with a silica glass cylinder of $R=100\ \mu$m and wall-
width $w=0.75\ \mu$m at $\lambda_{p}\approx 1.26\ \mu$m, which features a
roundtrip time $\tau=3.1$ ps and a quality factor
$\mathcal{Q}=\omega_{p}\tau/\gamma\approx 4.9\times 10^{6}$ ($\omega_{p}=2\pi
c/\lambda_{p}$), reasonable for cylinders sumcyl . The main higher dispersion
terms are given by $B_{3,0}\approx 9.26\times 10^{-10}$,
$B_{4,0}\approx-2.83\times 10^{-13}$, $B_{0,2}\approx 1.21\times 10^{-4}$. The
normalised detuning is $\delta=(\omega_{p}-\omega_{0})\tau/\gamma$.

The pump beam, assumed of Gaussian profile, has an amplitude $h_{0}$ and a
width $\sigma_{z}$. Under realistic conditions (cf. Fig.6), the pump beam is
also localised along $X$, and this is accounted for via the function $\xi(X)$,
where $N$ is the normalisation factor such that $\max(\xi(X))=1$.

Thermal effects are introduced via the light-to-phonon energy conversion,
$A=10^{-2}$, and the corresponding cooling rate, $B=5\times 10^{-2}$,
considering realistic values menyuckther . Thermal detuning, $\Delta_{T}$, is
assumed not to depend on $X$ since $\tau\sim$ps is much smaller than the
temperature diffusion time scale, in the order of the ns. Last, we note that
in Eq.3 we omitted a term $\sim\mu\partial_{z}\Delta_{T}$, accounting for the
heat diffusion along $z$, as its relative importance to the other terms is of
the order of $\sim 10^{-6}$ (the thermal diffusion coefficient is $\mu\approx
7.25\times 10^{-7}m^{2}/s$ kat ).

Raman scattering is introduced via the molecular vibrational field boyd ,
previously implemented in microresonators milianPRA , with standard parameter
for glass given by blow : Raman to Kerr effective fraction, $f_{R}=0.18$;
inverse phonon life-time, $\gamma_{R}=1/32$ fs-1; natural phonon frequency,
$\omega_{R}=1/12.2$ fs-1; and
$\Omega_{r}\equiv[\gamma_{R}^{2}+\omega_{R}^{2}]^{1/2}$.

Eq.1 in the main text is a particular case of the system Eqns.2-4 when higher
order dispersion, Raman scattering, thermal detuning, and pump azimuthal’s
localisation are disregarded ($\hat{D}_{hod}=Q=\Delta_{T}=0$, $\xi(X)=1$). The
time evolution of the above system, Eqns.2-5, is simulated via the fourth
order Runge-Kutta method.

Computation of stationary snakes. Stationary solutions (snakes or else) are
obtained numerically from Eq.1 or Eqns. 2-5 with the Newton-Raphson method in
the frame comoving with the nonlinear state, where they readily satisfy
$\partial_{t}\psi=0$. While equation 1 is already expressed in such frame,
Eqns.2-5 (expressed in the lab frame) are rewritten into the comoving frame
after the substitution $x=X-(b_{1,1}+v)t$, where $v$ is a velocity shift
induced by the higher order effects, which is computed together with the
nonlinear solution. When computing stationary states, the $x$-localisation of
the pump is disregarded (we set $\xi(X)=1$).

Stability of snakes. Linear stability analysis is performed to all stationary
solutions, represented by the tuple $\\{\psi_{s}\,Q_{s},\Delta_{T_{s}}\\}$.
Each field in the stationary solution is prone to develop instabilities. In
the initial stages of such instabilities fields are regarded as
$\psi=\psi_{s}+ae^{i\lambda t}+b^{*}e^{-i\lambda^{*}t}$, $Q=Q_{s}+ce^{i\lambda
t}+c^{*}e^{-i\lambda^{*}t}$, $\Delta_{T}=\Delta_{T_{s}}+de^{i\lambda
t}+d^{*}e^{-i\lambda^{*}t}$ (Note $Q,\Delta_{T}$ are real fields), where
$\lambda$ are the complex eigenvalues of the Jacobian matrix, obtained after
substituion of the above decomposition into the system of Eqns. 2-4 (with
$\xi(X)=1$) and linearising in $\lambda$. The real parts of $\lambda$,
$\mathrm{Re}(\lambda)$, are the growth rates yielding instabilities when
$\mathrm{Re}(\lambda)>0$.

Data availability

The data that support the plots within this paper and other findings of this
study are available from the corresponding author upon reasonable request.

Code availability

The analysis codes will be made available on reasonable request.

Acknowledgements

JAC and CM acknowledge support from the Spanish government via the Grant
PID2021-124618NB-C21 funded by MCIN/AEI/ 10.13039/501100011033 and by “ERDF A
way of making Europe”, by the “European Union”. CM acknowledges support from
Generalitat Valenciana PROMETEO/2021/082. PFC acknowledges partial support
from the Spanish gouvernement via the project PID2021-128676OB-I00 (MICINN).
LT acknowledges support by CEX2019-000910-S [MCIN/AEI/10.13039/501100011033],
Fundació Cellex, Fundació Mir Puig, and Generalitat de Catalunya (CERCA). YVK
academic research has been supported by the research project FFUU-2021-0003 of
the Institute of Spectroscopy of the Russian Academy of Sciences.

Author contributions

SBI and CM carried out the numerical simulations. CM conceived the project.
All authors contributed significantly to this work, discussed the results, and
contributed into the manuscript preparation.

Competing interests

The authors declare no competing interests.

## References

  * (1) Y. A. Demchenko and M. L. Gorodetsky, “Analytical estimates of eigenfrequencies, dispersion, and field distribution in whispering gallery resonators,” J. Opt. Soc. Am. B 30, 3056–3063 (2013).
  * (2) M. Ibanescu, S. G. Johnson, D. Roundy, C. Luo, Y. Fink, and J. D. Joannopoulos, “Anomalous Dispersion Relations by Symmetry Breaking in Axially Uniform Waveguides,” Phys. Rev. Lett. 92, 063903 (2004).
  * (3) M. Sumetsky, “Mode localization and the Q-factor of a cylindrical microresonator,” Opt. Lett. 35, 2385–2387 (2010).
  * (4) A. Leshem, Z. Qi, T. F. Carruthers, C. R. Menyuk, and O. Gat, “Thermal instabilities, frequency-comb formation, and temporal oscillations in Kerr microresonators,” Phys. Rev. A 103, 013512 (2021).
  * (5) T. Katsura, “Thermal Diffusivity of Silica Glass at Pressures up to 9 GPa,” Phys. Chem. Min. 20, 201-208 (1993).
  * (6) R. Boyd, “Nonlinear Optics” (Academic Press, New York, 2008).
  * (7) C. Milián, A. V. Gorbach, M. Taki, A. V. Yulin, and D. V. Skryabin, “Solitons and frequency combs in silica microring resonators: Interplay of the Raman and higher-order dispersion effects,” Phys. Rev. A 92, 033851 (2015).
  * (8) K.J. Blow and D. Wood, “Theoretical description of transient stimulated Raman scattering in optical fibers,” IEEE J. Quantum Electron. 25, 2665 (1989).