# Spatial-Temporal Decoupling Contrastive Learning for Skeleton-based Human
Action Recognition

Shaojie Zhang1 Jianqin Yin1,2111Concat Author Yonghao Dang1 1Beijing
University of Posts and Telecommunications
2Queen Mary School Hainan, Beijing University of Posts and Telecommunications
{zsj, jqyin<EMAIL_ADDRESS>
###### Abstract

Skeleton-based action recognition is a central task of human-computer
interaction. Current methods apply the modeling paradigm of image recognition
to it directly. However, the skeleton sequences abstracted from the human body
is a sparse representation. The features extracted from the skeleton encoder
are spatiotemporal decoupled, which may confuse the semantics. To reduce the
coupling and improve the semantics of the global features, we propose a
framework (STD-CL) for skeleton-based action recognition. We first decouple
the spatial-specific and temporal-specific features from the spatiotemporal
features. Then we apply the attentive features to contrastive learning, which
pulls together the features from the positive pairs and pushes away the
feature embedding from the negative pairs. Moreover, the proposed training
strategy STD-CL can be incorporated into current methods. Without additional
compute consumption in the testing phase, our STD-CL with four various
backbones (HCN, 2S-AGCN, CTR-GCN, and Hyperformer) achieves improvement on
NTU60, NTU120, and NW-UCLA benchmarks. We will release our code at:
https://github.com/BUPTSJZhang/STD-CL.

##  1 Introduction

Human action recognition is an active task in the computer vision area. Due to
its wide range of applications in human-computer interaction, video analysis,
virtual reality, and so on, this task has been researched extensively in the
past decade Ren et al. (2020). In recent years, with the development of depth
sensors Cao et al. (2017) and human pose estimation algorithms Dang et al.
(2019), the skeleton sequence consisting of coordinates of human joints can be
easily acquired. The so-called skeleton-based action recognition utilizing
human joints has attracted much interest due to its robustness to background
clutter, illumination, and viewpoint changes Zhang et al. (2019).

Figure 1: Comparison of the image recognition framework (top) with the action
recognition framework (down) in both the perspective of humans and Deep Neural
Network (DNN).

Generally, most of the previous approaches Shi et al. (2019); Chen et al.
(2021a, b) in skeleton-based action recognition follow a paradigm as shown in
Figure 1. A skeleton encoder extracts the spatiotemporal co-occurrence
feature, and a softmax linear classifier is used to project the feature to the
category distributions. However, the global features may mix the
spatiotemporal information, which confuses the spatial- and temporal-specific
information reflecting different semantics at the joint level and frame level.
For example, when distinguishing between the action ”put on a hat” and ”take
off a hat”, which have more differences in the frame-level dimension,
temporal-specific information is more important. But for the action ”point to
something” and ”throw”, they have more spatial differences in the joint-level
dimension, so spatial-specific information is more important to distinguish
them.

The key to action recognition is matching the motion patterns to the action
categories. When our humans do this job, we can infer the actions by
individually matching the spatial context information or the temporal context
information to the actions we have seen before. Thus, the motion patterns
might be decoupled into spatial motion patterns and temporal motion patterns.
The same action performed by different persons may have similar
spatiotemporal-decoupling attentive features. This motivates us to model
cross-sequence spatiotemporal relations of human skeleton sequences.

In recent years, contrastive learning He et al. (2020); Chen et al. (2020) has
achieved great success in self-supervised representation learning, due to its
strong ability to study cross-sample relations. It can pull together the
features from positive pairs and push away the features from negative pairs.

Based on the above analysis, contrasting spatial-specific features and
temporal-specific features via contrastive learning can reduce coupling and
improve the semantic information in the global spatiotemporal features.
Therefore, we propose a Spatial-Temporal Decoupling Contrastive Learning (STD-
CL) framework for skeleton-based action recognition. In STD-CL, as shown in
Figure 1, we design a Spatial-Temporal Feature Decoupling(STFD) module to
decouple the global spatiotemporal feature to spatial attentive feature and
temporal attentive feature. We contrast two memory banks to store the
decoupled spatial and temporal features respectively. In the training phase,
we sample the positives and negatives from the memory banks according to the
labels, and apply them in contrastive learning. Moreover, there is no
additional computational consumption during the testing phase.

Overall, our main contributions can be summarized as follows:

  * •

To reduce the spatial and temporal coupling and improve the semantics of the
co-occurrence features extracted from the skeleton encoder, we propose a novel
STD-CL framework to learn more abundant global features from the sequence.

  * •

Motivated by the powerful capability of contrastive learning for modeling
cross-samples, we apply it for supervised skeleton-based action recognition.
STD-CL builds two memory banks to model the cross-sequence spatiotemporal
relations. Moreover, our STD-CL can be seamlessly incorporated into current
approaches.

  * •

Extensive experiments show that STD-CL achieves significant improvements based
on four state-of-the-art methods (HCN, 2S-AGCN CTR-GCN, and Hyperformer) on
three benchmark datasets (NTU60, NTU120, and NW-UCLA). Moreover, the STD-Cl is
only used in the training phase without additional overhead in the testing
phase.

Figure 2: The framework of our proposed Spatial-Temporal Decoupling
Contrastive Learning (STD-CL). In the testing phase, The skeleton encoder
extracts spatiotemporal feature $\boldsymbol{X}$ from the input skeleton
sequence $\boldsymbol{S}$. A global average pooling (GAP) layer is used to
aggregate the global information. A classifier is adopted to project the
global feature $\boldsymbol{g}$ to the category distribution. In the training
phase, the feature $\boldsymbol{X}$ is fed into the Spatial-Temporal Feature
Decoupling (STFD) module to decouple spatial attentive feature
$\boldsymbol{s}$ and temporal attentive feature $\boldsymbol{t}$. The details
of STFD are shown on the right. Two memory banks are built to store the
decoupled features respectively. We contrast the positive and negative
instance features sampled from the memory banks according to the labels.

##  2 Related Work

###  2.1 Skeleton-based Action Recognition

Skeleton-based action recognition aims to classify actions from sequences of
human keypoints. In early studies, RNN was a natural choice to handle the
skeleton sequence. HBRNN Du et al. (2015) applied a hierarchical RNN
architecture to model the long-range temporal dependency of the skeleton
sequences. STA-LSTM Song et al. (2017) proposed an end-to-end spatial and
temporal attention model to focus on discriminative joints of the skeleton
within each frame of the inputs. Inspired by the success of CNN in image
recognition, CNN-based methods have also been studied. HCN Li et al. (2018)
treated the joint dimension as channels to aggregate different levels of
spatiotemporal context information. TA-CNN Xu et al. (2022) adopted a pure CNN
architecture to model the irregular skeleton topology.

The human body can be abstracted as a graph structure, which treats the joint
as a node and the bone between two joints as an edge. Thus, GCN-based methods
were introduced to this task. ST-GCN Yan et al. (2018) adopted GCN on the
predefined spatial-temporal graphs to model the relations between joints.
2S-AGCN Shi et al. (2019) modeled the correlation between two joints given
corresponding features with the self-attention mechanism. CTR-GCN Chen et al.
(2021a) proposed a channel-wise topology refinement graph convolution to model
fine-grained relation. InfoGCN Chi et al. (2022) combined a learning objective
and an encoding method to break the information bottleneck.

With the popularity of transformers in computer vision, transformer-based
methods have also been investigated for this task. ST-TR Plizzari et al.
(2021) proposed a two-stream transformer architecture to model the spatial and
temporal dimensions respectively. Hyperformer Zhou et al. (2022) applied
transformer to model spatiotemporal features from the hyper-graph of the
sequences.

However, most of the previous methods adopted the training framework we have
discussed in section 1, confusing the spatiotemporal information and semantics
in the global features.

###  2.2 Contrastive Learning

Recently, contrastive learning has been widely used in self-supervised
representation learning. The core idea is to learn the discriminative features
by pulling the positive pairs and pushing away the negative pairs in the
embedded space. The positive pairs are features from transformed versions of a
sample with different augmentations, and the negative pairs are from different
samples. SimCLR Chen et al. (2020) proposed to adopt a huge batch size to
learn more useful representations. MoCo He et al. (2020) used a momentum
encoder with a memory bank to store more negative instances. Since then,
memory bank mechanisms have been the paradigm for constructing abundant
negative pairs to obtain discriminative features.

In the field of skeleton-based action recognition, contrastive learning was
usually adopted to learn the invariance from the skeleton sequences for
unsupervised learning. Moreover, several works are applying contrastive
learning to supervised action recognition. GAP Xiang et al. (2023) proposed a
framework for contrasting the features from visual and text modalities to
enhance the representation by using knowledge about actions and human body
parts. SkeletonGCL Huang et al. (2023) explicitly explored the rich cross-
sequence relations by using graph contrastive learning. However, SkeletonGCL
can only be incorporated into the GCN-based methods, which limits its
generalization. FR-Head Zhou et al. (2023) employed multi-level contrastive
learning in GCN-based methods to distinguish ambiguous actions.

Compared with the methods mentioned above, we contrast the spatial decoupled
representations and temporal decoupled representations instead of the global
pooled features or the graph embeddings. Thus, our STD-CL can be incorporated
into the previous methods, no matter what model-based, without introducing
data from other modalities.

##  3 Methodology

###  3.1 Perliminary

####  3.1.1 Notations

In our work, bold-face lower/upper-case letters denote vectors and tensors,
respectively, and light-face lower/upper-case letters mean scalars.

####  3.1.2 Background

Skeleton-based action recognition is to predict the categories from the input
sequences. We define the human skeleton sequence with $J$ keypoints and
$T_{0}$ frames in 3D space as $\boldsymbol{S}\in\mathbb{R}^{J\times
T_{0}\times 3}$. The sequence is fed into a skeleton encoder $\mathcal{F}$ to
extract the features $\boldsymbol{X}\in\mathbb{R}^{J\times T\times C}$. The
skeleton encoder aggregates the spatial information and the temporal
information to obtain discriminative spatiotemporal features. After fully
spatiotemporal aggregations, a global average pooling (GAP) layer is adapted
to summarize the global features. Finally, a softmax-based fully-connected
(FC) layer maps the global feature to a probability prediction
$\hat{\boldsymbol{y}}\in\mathbb{R}^{K}$ of K candidate categories. The process
can be defined as:

$\boldsymbol{X}=\mathcal{F}(\boldsymbol{S})$ (1)
$\boldsymbol{y}=FC(GAP(\boldsymbol{X}))$ (2)

In the process of optimization, a standard cross-entropy loss is applied to
supervise the prediction distribution and the ground truth $\boldsymbol{y}$ as
follows:

$\mathcal{L}_{CE}=-\sum_{i}\boldsymbol{y}_{i}log\hat{\boldsymbol{y}_{i}}$ (3)

###  3.2 Spatial Temporal Feature Decoupling

Human actions have the property of spatio-temporal coupling, the feature
$\boldsymbol{X}$ after extraction contains complex spatial and temporal
relations, which may confuse the semantics. Thus, we decouple the feature to
spatial attentive feature and temporal attentive feature, which are more
discriminative features.

The details of the proposed STFD are shown in (b) of Figure 2. The feature
$\boldsymbol{X}$ is fed into two parallel branches for feature decoupling. In
the first spatial feature decoupling (SFD) module, a temporal pooling layer is
adopted to obtain the spatial attentive feature
$\boldsymbol{X}_{s}\in\mathbb{R}^{T\times C}$. The temporal attentive feature
$\boldsymbol{X}_{t}\in\mathbb{R}^{J\times C}$ is obtained through a spatial
pooling layer. Then, two linear transformation functions $\theta$ and $\phi$
with a reduction rate of $r$ are used to transform the attentive features into
a neatly compact representation as follows:

$\left\\{\begin{matrix}\boldsymbol{X}_{s}^{*}=\theta(\boldsymbol{X}_{s})=\boldsymbol{X}_{s}\boldsymbol{W}_{\theta}\in\mathbb{R}^{J\times
C/r}\\\ \\\
\boldsymbol{X}_{t}^{*}=\phi(\boldsymbol{X}_{t})=\boldsymbol{X}_{t}\boldsymbol{W}_{\phi}\in\mathbb{R}^{T\times
C/r}\end{matrix}\right.$ (4)

where, $\boldsymbol{W}_{\theta},\boldsymbol{W}_{\phi}\in R^{C\times C/r}$ are
the weights of the spatial and temporal transformation function respectively.
After this, to retain the attentive information, we squeeze the spatial
dimension into the channel dimension in the SFD module, and the TFD module is
similar. In the end, two liner transformation functions $\zeta$ and $\sigma$
are used to transform the features to the latent space for further contrastive
learning, which can be formulated as:

$\left\\{\begin{matrix}\boldsymbol{s}=\zeta(\boldsymbol{s}^{*})=\boldsymbol{s}^{*}\boldsymbol{W}_{\zeta}\in\mathbb{R}^{D}\\\
\\\
\boldsymbol{t}=\sigma(\boldsymbol{t}^{*})=\boldsymbol{t}^{*}\boldsymbol{W}_{\sigma}\in\mathbb{R}^{D}\end{matrix}\right.$
(5)

where, $\boldsymbol{W}_{\zeta}\in\mathbb{R}^{JC/r\times D}$,
$\boldsymbol{W}_{\sigma}\in\mathbb{R}^{TC/r\times D}$ are the weights of final
projection function. In this way, we decouple the spatial attentive feature
$\boldsymbol{s}$ and the temporal attentive feature $\boldsymbol{t}$ from the
spatiotemporal feature $\boldsymbol{X}$.

###  3.3 Spatial-Temporal Decoupling Contrastive Learning

To improve the semantics of the action representations, we apply the decoupled
features to contrastive learning. The core is to reduce intra-class distances
and increase inter-class distances of the spatial and temporal attentive
features.

####  3.3.1 Memory Bank

To obtain abundant negative pairs, we contrast two memory banks
$\mathcal{M}_{spa}$ and $\mathcal{M}_{tem}$ to store the decoupled attentive
features. Moreover, to avoid introducing hyperparameters and encourage the
model to fully explore the cross-sequence context, we set the two memory banks
to store the attentive features of all samples. Thus, the memory banks
$\mathcal{M}_{spa},\mathcal{M}_{tem}\in\mathbb{R}^{L\times D}$. $L$ is the
number of the instance of the dataset. When an instance is fed into the model
to extract features, the stored features in the memory banks are updated
according to the index.

$\left\\{\begin{matrix}\mathcal{M}_{spa}[i]=\boldsymbol{s}_{i}\\\ \\\
\mathcal{M}_{tem}[i]=\boldsymbol{t}_{i}\end{matrix}\right.$ (6)

where $i$ is the index of the instance in the dataset. Each element in
$\mathcal{M}_{s}pa$ and $\mathcal{M}_{t}em$ indicates a decoupled feature
embedding from an instance. In the memory banks, the corresponding labels with
the decoupled features are stored for the feature selection. Particularly, the
decoupled features with the same category label as the current instance are
selected as positive pairs. And negative pairs are sampled from the memory
banks with other labels.

####  3.3.2 Training Objective

The positive pairs and negative pairs sampled from the two memory banks are
used to achieve spatial-temporal decoupled contrastive learning. To pull
together the positive pairs and push away the negative pairs, we define the
similarity between the two vectors with cosine distances as follows:

$\left\langle\boldsymbol{u},\boldsymbol{v}\right\rangle=\frac{\boldsymbol{u}\cdot\boldsymbol{v}^{T}}{||\boldsymbol{u}||\;||\boldsymbol{v}||}$
(7)

where $\boldsymbol{u}$, $\boldsymbol{v}$ are the features sampled from the
memory banks. Then the InfoNCE loss adopted for spatial-specific feature
contrasting and temporal-specific feature contrasting can be written as
follows:

$\mathcal{L}_{NCE}^{spa}=-\sum\limits_{\boldsymbol{s}^{+}\in\mathcal{M}^{+}_{spa}}log\frac{\left\langle\boldsymbol{s},\boldsymbol{s}^{+}\right\rangle/\tau}{\left\langle\boldsymbol{s},\boldsymbol{s}^{+}\right\rangle/\tau+\sum\limits_{\boldsymbol{s}^{-}\in\mathcal{M}^{-}_{spa}}\left\langle\boldsymbol{s},\boldsymbol{s}^{-}\right\rangle/\tau}$

(8)

$\mathcal{L}_{NCE}^{tem}=-\sum\limits_{\boldsymbol{t}^{+}\in\mathcal{M}^{+}_{tem}}log\frac{\left\langle\boldsymbol{t},\boldsymbol{t}^{+}\right\rangle/\tau}{\left\langle\boldsymbol{t},\boldsymbol{t}^{+}\right\rangle/\tau+\sum\limits_{\boldsymbol{t}^{-}\in\mathcal{M}^{-}_{tem}}\left\langle\boldsymbol{t},\boldsymbol{t}^{-}\right\rangle/\tau}$

(9)

where $\tau$ is a temperature hyperparameter in contrastive learning.
$\mathcal{M}^{+}_{spa}$ and $\mathcal{M}^{+}_{tem}$ are the sets of the
decoupled features from the positive sequences with the same label as the
current feature in the memory banks. $\mathcal{M}^{-}_{spa}$ and
$\mathcal{M}^{-}_{tem}$ are the sets of negative features with different
lalbes.

Finally, the total loss function of our proposed STD-CL is defined as follows:

$\mathcal{L}=\mathcal{L}_{CE}+\mathcal{L}_{NCE}^{spa}+\mathcal{L}_{NCE}^{tem}$
(10)

####  3.3.3 Sampling Strategy

To achieve efficient contrast across sequences, the strategy for sampling
features from memory banks for contrastive learning has a large impact on the
experimental performance. Therefore, in this paper, we follow the sampling
strategy in Huang et al. (2023), which combines the hard mining strategy and
random sampling strategy. The hard mining strategy is introduced to focus on
the hard samples, which means the $N^{+}_{H}$ hardest positive samples with
the lowest similarity and $N^{-}_{H}$ hardest negative samples with the
highest similarity in equation 7. Furthermore, the random sampling strategy is
adopted to maintain global random contrast. There are $N^{-}_{R}$ random
negative samples selected to contrast with the input instance. In total, the
sampling strategy adopted in our STD-CL takes contrast efficiency into account
and contributes to the hard examples.

##  4 Experiment

Table 1: Top-1 action aclssification accuracy comparison ($\%$) with the
state-of-the-art methods on NTU 60 and NTU120 datasets. The numbers in gray
indicate the results reported in their papers. Particularly, the 2S incidents
the ensemble results of joint and bone modality. And the 4S incidents the
ensemble results of joint, bone, joint motion and bone motion modality. *
denotes the ensemble results of the streams reported in their papers.

Dataset | NTU 60 | NTU 120
---|---|---
Setting | X-Sub | X-View | X-Sub | X-Set
Method/Modality | J | B | 2S | 4S | J | B | 2S | 4S | J | B | 2S | 4S | J | B | 2S | 4S
VA-LSTM Zhang et al. (2017) | - | - | 79.4* | - | - | - | 87.6* | - | - | - | - | - | - | - | - | -
AGC-LSTM Si et al. (2019) | - | - | 89.2* | - | - | - | 95.0* | - | - | - | - | - | - | - | - | -
ST-GCN Yan et al. (2018) | 81.5 | - | - | - | 88.3 | - | - | - | - | - | - | - | - | - | - | -
SGCN Zhang et al. (2020) | - | - | 89.0 | - | - | - | 94.5 | - | - | - | 79.2 | - | - | - | 81.5 | -
Shift-GCN Cheng et al. (2020b) | 87.8 | - | 89.7 | 90.7 | 95.1 | - | 96.0 | 96.5 | 80.9 | - | 85.3 | 85.9 | 83.2 | - | 86.6 | 87.6
DC-GCN+ADG Cheng et al. (2020a) | - | - | 90.8 | - | - | - | 96.6 | - | - | - | 86.5 | - | - | - | 88.1 | -
Dynamic GCN Ye et al. (2020) | - | - | - | 91.5 | - | - | - | 96.0 | - | - | - | 87.3 | - | - | - | 88.6
MS-G3D Liu et al. (2020) | 89.4 | 90.1 | 91.5 | - | 95.0 | 95.3 | 96.2 | - | - | - | 86.9 | - | - | - | 88.4 | -
DSTA Shi et al. (2020) | - | - | 91.5 | - | - | - | 96.4 | - | - | - | 86.6 | - | - | - | 88.0 | -
MST-GCN Chen et al. (2021b) | 89.0 | 89.5 | 91.1 | 91.5 | 95.1 | 95.2 | 96.4 | 96.6 | 82.8 | 84.8 | 87.0 | 87.5 | 84.5 | 86.3 | 88.3 | 88.8
ST-TR Plizzari et al. (2021) | - | - | 89.9 | - | - | - | 96.1 | - | - | - | 82.7 | - | - | - | 84.7 | -
TA-CNN Xu et al. (2022) | 88.9 | 89.2 | 91.0 | 91.5 | 94.5 | 94.1 | 95.7 | 95.9 | 84.0 | 85.1 | 87.8 | 88.2 | 85.3 | 86.3 | 89.0 | 89.6
EfficientGCN-B4 Song et al. (2022) | - | - | 91.7 | - | - | - | 95.7 | - | - | - | 88.3 | - | - | - | 89.1 | -
InfoGCN Chi et al. (2022) | 89.4 | 90.6 | 91.3 | 92.3 | 95.2 | 95.4 | 96.2 | 96.7 | 84.2 | 86.9 | 88.2 | 89.2 | 86.3 | 88.5 | 89.4 | 90.7
HCN Li et al. (2018) | - | - | 84.3* | - | - | - | 89.9* | - | - | - | - | - | - | - | - | -
HCN w/ STD-CL | - | - | -* | - | - | - | -* | - | - | - | - | - | - | - | - | -
2S-AGCN Shi et al. (2019) | 88.9 | 89.2 | 91.0 | 91.5 | 94.5 | 94.1 | 95.7 | 95.9 | 82.9 | 85.1 | 87.8 | 88.2 | 85.3 | 86.3 | 89.0 | 89.6
2S-AGCN w/STD-CL | 89.4 | 89.8 | 91.5 | 91.9 | 94.7 | 94.6 | 96.1 | 96.3 | 84.1 | 85.1 | 87.8 | 88.2 | 85.3 | 86.3 | 89.0 | 89.6
CTR-GCN Chen et al. (2021a) | 89.8 | 90.2 | 92.0 | 92.4 | 94.8 | 94.8 | 96.3 | 96.8 | 84.9 | 85.7 | 88.7 | 88.9 | 86.7 | 87.5 | 90.1 | 90.5
CTR-GCN w/STD-CL | 90.4 | 90.8 | 92.2 | 92.7 | 95.2 | 94.9 | 96.4 | 96.8 | 85.5 | 85.7 | 88.7 | 88.9 | 86.7 | 87.5 | 90.1 | 90.5
Hyperformer Zhou et al. (2022) | 90.7 | 91.2 | 92.4 | 92.8 | 95.1 | 95.2 | 96.2 | 96.5 | 86.6 | 88.1 | 89.3 | 89.9 | 88.0 | 89.0 | 90.6 | 91.2
Hyperformer w/STD-CL | - | - | - | - | - | - | - | - | - | - | - | - | - | - | - | -

Table 2: Top-1 action aclssification accuracy comparison ($\%$) with the
state-of-the-art methods on the NW-UCLA dataset. The numbers in gray indicate
the results reported in their papers. Particularly, the 2S incidents the
ensemble results of joint and bone modality. And the 4S incidents the ensemble
results of joint, bone, joint motion and bone motion modality. * denotes the
ensemble results of the streams reported in their papers.

Dataset | NW-UCLA
---|---
Method/Modality | J | B | 2S | 4S
AGC-LSTM Si et al. (2019) | - | - | 93.3* | -
DC-GCN+ADG Cheng et al. (2020a) | - | - | 95.3 | -
Shift-GCN Cheng et al. (2020b) | 92.5 | - | 94.2 | 94.6
2S-AGCN Shi et al. (2019) | 92.0 | 92.2 | 95.0 | 95.5
2S-AGCN w/STD-CL | 93.8 | 93.1 | 95.9 | 97.0
CTR-GCN Chen et al. (2021a) | 94.6 | 91.8 | 94.2 | 96.5
CTR-GCN w/STD-CL | 94.8 | 94.6 | 95.7 | 97.2
InfoGCN Chi et al. (2022) | 93.8 | 94.2 | 95.5 | 96.1
Info-GCN w/STD-CL | 94.2 | 95.5 | 95.7 | 96.3
Hyperformer Zhou et al. (2022) | 92.7 | 95.0 | 95.0 | 96.6
Hyperformer w/STD-CL | 95.3 | 95.3 | 96.3 | 97.0

###  4.1 Dataset Settings

NTU RGB+D. NTU RGB+D (NTU60) Shahroudy et al. (2016) is a large-scale dataset
for skeleton-based action recognition, which contains 56880 videos performed
by 40 volunteers. The action sequences can be categorized into 60 classes.
Following the authors of this dataset recommendation, we process this dataset
into two benchmarks: cross-subject(X-sub) and cross-view(X-view). In the
cross-subject setting, sequences of 20 subjects are for training, and the
sequences of the rest 20 subjects are for validation. In the cross-view
setting, skeleton sequences are split by camera views. Samples from two camera
views are used for training, and the rest are used for evaluation.

NTU RGB+D 120. NTU RGB+D 120 (NTU120) dataset Liu et al. (2019) adds 57367 new
skeleton sequences and 60 new action classes to the original NTU RGB+D 60
dataset. There are 32 various video configurations in it, each of which
depicts a different location and background. The authors offered the cross-
subject(X-sub) and cross-setup(X-set) as two benchmark evaluations. In the
cross-subject setting, sequences from 53 subjects are for training, and
sequences from the other 53 subjects are for testing. In the cross-setup
setting, skeleton sequences are split by setup ID. Samples from even set-up
IDs are used for training, and the odd setup IDs are used for evaluation.

Northwestern-UCLA. Northwestern-UCLA (NW-UCLA) dataset Wang et al. (2014)
contains 1494 video clips of 10 different actions captured from three Kinect
cameras. We follow the same evaluation protocol in Wang et al. (2014): the
first two cameras for training and the other for testing.

###  4.2 Implementation Details

To fully validate the effectiveness and generalizability of STGCL, we take
four different model-based approaches as baseline models. (1) CNN-based
method. We select HCN Li et al. (2018) as the CNN-based method to validate
STD-CL. We reproduce the experimental results using the released code and
conduct our STD-CL on it. And we follow the training recipes as the paper
describes. (2) GCN-based method. 2S-AGCN Shi et al. (2019) and CTR-GCN Chen et
al. (2021a), the widely used GCN-based model for skeleton-based action
recognition. For CTR-GCN, we follow their training strategies. Particularly,
for 2S-AGCN, we update the training strategies from CTR-GCN, which can improve
the performance significantly. (3) Transformer-based method. For transformer-
based methods, we choose Hyperformer as the baseline.

For the details of contrastive learning, the size of the established memory
banks is the length of the dataset. The dimensions of the decoupled features
$\boldsymbol{t}$ and $\boldsymbol{s}$ is 256. For all datasets, the number of
sampling instances $N^{+}_{H}$, $N^{-}_{H}$ and $N^{-}_{R}$ are set as 128,
512, and 512, respectively. To select the model with the best performance, we
experiment with temperatures $\tau$ of 0.5, 0.8, and 1.0. All experiments are
conducted on four GeForce GTX NVIDIA 3090 GPUs. More implementation details
are presented in the supplementary material.

###  4.3 Comparsion With The State-of-The-Art

We compare our method with previous SOTA methods in Tables 1 and 2. To fully
validate the effectiveness of the proposed STD-CL, we compare the experimental
results in a single modality and the ensemble results of multi-modality. $J$
means the joint stream, $B$ denotes the bone stream, $2S$ denotes the ensemble
results of the above two streams in default, and $4S$ incidents the ensemble
results of additional joint motion and bone motion streams. Particularly, due
to the design of their model architecture, the experimental results of VA-LSTM
Zhang et al. (2017), ACG-LSTM Si et al. (2019) and HCN Li et al. (2018) are
the ensemble results of two streams. Therefore, we classify the results into
$2S$.

As shown in Table 1, on NTU60 and NTU120,

As shown in Table 2, on NW-UCLA, STD-CL outperforms 2S-AGCN by 1.5%. It also
outperforms CTR-CCN by 0.7%. And it outperforms the recent work InfoGCN and
Hyperformer by 0.2% and 0.4%, respectively. Considering that the performance
on this dataset is already very high, such improvement brought by STD-CL is
significant.

###  4.4 Comparsion With Other Contrastive Learning Method

Table 3: Top-1 accuracy comparison ($\%$) with recent contrastive learning
methods in supervised skeleton-based action recognition on the cross-subject
benchmark of NTU120 dataset.

Method | NTU120 Xsub
---|---
Baseline Chen et al. (2021a) | 84.9
GAP Xiang et al. (2023) | 85.5
SkeletonGCL Huang et al. (2023) | -
FR Head Zhou et al. (2023) | 85.5
STD-CL | 85.5
STD-CL(Hyperformer Zhou et al. (2022)) | 86.6+

We compare our method with the recent contrastive learning methods in
supervised skeleton-based action recognition. The comparison is conducted on
the cross-subject of the NTU120 dataset. As shown in Table 3, our STD-CL
achieves competitive improvements. Furthermore, most previous methods are only
compatible with GCN-based methods. To take advantage of our method, we combine
the STD-CL with Hyperformer, which is a transformer-based method. Our method
outperforms previous contrastive learning methods by 1.1%, which is a
promising result. For more comparisons, please refer to the supplementary
materials.

###  4.5 Ablation Study

In this section, we conduct experiments to evaluate the different experimental
settings on the cross-subject benchmark of the NTU120 dataset to verify the
design of the proposed STD-CL. The experimental results of ablation studies
are shown in Table 4.

Hyperparameter in STDF module.

Contrastive Strategies.

Performance on Ambiguous Actions.

Training consumption.

(a) Impact of the channel reduction rate in STFD. Method | $r$ | Acc(%)
---|---|---
CTR-GCN | - | 84.9
4 | -
8 | 85.5
16 | -

(b) Impact of the dimension of the decoupled features. Method | $D$ | Acc(%)
---|---|---
CTR-GCN | - | 84.9
128 | -
256 | 85.5
512 | -

(c) Contrastive strategies Method | Acc(%)
---|---
Baseline (w/o contrast) | 84.9
w/SD-CL | -
w/TD-CL | -
w/G-CL | -
w/STD-CL w/ G-CL | -
STD-CL | 85.5

(d) Accuracy on different level actions. Setting | Acc(%)
---|---
Method/Level | Hard | Medium | Easy
2S-AGCN | 58.9 | 82.0 | 95.0
2S-AGCN w/ours | - | - | -
CTR-AGCN | 59.6 | 82.4 | 95.1
CTR-AGCN w/ours | - | - | -
Hyperformer | 58.9 | 82.0 | 95.0
Hyperformer w/ours | - | - | -

(e) The comparison of training consumption. Method | Memory-Usage (G) | Time (hours)
---|---|---
HCN | - | -
HCN w/ours | - | -
2S-AGCN | - | -
2S-AGCN w/ours | - | -
CTR-GCN | - | -
CTR-AGCN w/ours | - | -
Hyperformer | - | -
Hyperformer w/ours | - | -

Table 4: Ablation study with different experimental settings of STD-CL on
NTU120. The Acc incidents the accuracy of the action recognition. SD-CL means
applying spatial decoupled features to contrastive learning, and TD-CL means
applying temporal decoupled features to contrastive learning. G-CL means
applying the global features, which is the output of the global mean pooling
layer, to contrastive learning.

###  4.6 Discussions

Figure 3: Action class with accuracy differences higher than 4% between CTR-
GCN method in the cross-subject benchmark of NTU120 dataset.

##  5 Conclusion

In this paper, we propose a novel spatial-temporal decoupling contrastive
learning (STD-CL) framework for skeleton-based action recognition, which can
be combined with most of the previous skeleton encoders. We employ contrastive
learning to reduce intra-class distances and increase inter-class distances of
the spatial-specific features and temporal-specific features. In this way, the
spatiotemporal coupling in the global features decreases, and the semantics
increase, which is beneficial to action recognition. The extensive
experimental results with four various skeleton encoders on three datastes
demonstrate the effectiveness of the proposed framework.

## Acknowledgments

This work was supported partly by the National Natural Science Foundation of
China (Grant No. 62173045, 62273054), partly by the Fundamental Research Funds
for the Central Universities (Grant No. 2020XD-A04-3), and the Natural Science
Foundation of Hainan Province (Grant No. 622RC675).

## References

  * Cao et al. [2017] Zhe Cao, Tomas Simon, Shih-En Wei, and Yaser Sheikh.  Realtime multi-person 2d pose estimation using part affinity fields.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 7291–7299, 2017.
  * Chen et al. [2020] Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton.  A simple framework for contrastive learning of visual representations.  In Proceedings of the International Conference on Machine Learning, pages 1597–1607. PMLR, 2020.
  * Chen et al. [2021a] Yuxin Chen, Ziqi Zhang, Chunfeng Yuan, Bing Li, Ying Deng, and Weiming Hu.  Channel-wise topology refinement graph convolution for skeleton-based action recognition.  In Proceedings of the IEEE International Conference on Computer Vision, pages 13359–13368, 2021.
  * Chen et al. [2021b] Zhan Chen, Sicheng Li, Bing Yang, Qinghan Li, and Hong Liu.  Multi-scale spatial temporal graph convolutional network for skeleton-based action recognition.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 1113–1122, 2021.
  * Cheng et al. [2020a] Ke Cheng, Yifan Zhang, Congqi Cao, Lei Shi, Jian Cheng, and Hanqing Lu.  Decoupling gcn with dropgraph module for skeleton-based action recognition.  In Proceedings of the European Conference on Computer Vision, pages 536–553. Springer, 2020.
  * Cheng et al. [2020b] Ke Cheng, Yifan Zhang, Xiangyu He, Weihan Chen, Jian Cheng, and Hanqing Lu.  Skeleton-based action recognition with shift graph convolutional network.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 183–192, 2020.
  * Chi et al. [2022] Hyung-gun Chi, Myoung Hoon Ha, Seunggeun Chi, Sang Wan Lee, Qixing Huang, and Karthik Ramani.  Infogcn: Representation learning for human skeleton-based action recognition.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 20186–20196, 2022.
  * Dang et al. [2019] Qi Dang, Jianqin Yin, Bin Wang, and Wenqing Zheng.  Deep learning based 2d human pose estimation: A survey.  Tsinghua Science and Technology, 24(6):663–676, 2019.
  * Du et al. [2015] Yong Du, Wei Wang, and Liang Wang.  Hierarchical recurrent neural network for skeleton based action recognition.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1110–1118, 2015.
  * He et al. [2020] Kaiming He, Haoqi Fan, Yuxin Wu, Saining Xie, and Ross Girshick.  Momentum contrast for unsupervised visual representation learning.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 9729–9738, 2020.
  * Huang et al. [2023] Xiaohu Huang, Hao Zhou, Bin Feng, Xinggang Wang, Wenyu Liu, Jian Wang, Haocheng Feng, Junyu Han, Errui Ding, and Jingdong Wang.  Graph contrastive learning for skeleton-based action recognition.  2023\.
  * Li et al. [2018] Chao Li, Qiaoyong Zhong, Di Xie, and Shiliang Pu.  Co-occurrence feature learning from skeleton data for action recognition and detection with hierarchical aggregation.  roceedings of the International Joint Conference on Artificial Intelligence., 2018.
  * Liu et al. [2019] Jun Liu, Amir Shahroudy, Mauricio Perez, Gang Wang, Ling-Yu Duan, and Alex C Kot.  Ntu rgb+ d 120: A large-scale benchmark for 3d human activity understanding.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 42(10):2684–2701, 2019.
  * Liu et al. [2020] Ziyu Liu, Hongwen Zhang, Zhenghao Chen, Zhiyong Wang, and Wanli Ouyang.  Disentangling and unifying graph convolutions for skeleton-based action recognition.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 143–152, 2020.
  * Plizzari et al. [2021] Chiara Plizzari, Marco Cannici, and Matteo Matteucci.  Spatial temporal transformer network for skeleton-based action recognition.  In Pattern Recognition. ICPR International Workshops and Challenges: Virtual Event, January 10–15, 2021, pages 694–701. Springer, 2021.
  * Ren et al. [2020] Bin Ren, Mengyuan Liu, Runwei Ding, and Hong Liu.  A survey on 3d skeleton-based action recognition using learning method.  arXiv preprint arXiv:2002.05907, 2020.
  * Shahroudy et al. [2016] Amir Shahroudy, Jun Liu, Tian-Tsong Ng, and Gang Wang.  Ntu rgb+ d: A large scale dataset for 3d human activity analysis.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1010–1019, 2016.
  * Shi et al. [2019] Lei Shi, Yifan Zhang, Jian Cheng, and Hanqing Lu.  Two-stream adaptive graph convolutional networks for skeleton-based action recognition.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 12026–12035, 2019.
  * Shi et al. [2020] Lei Shi, Yifan Zhang, Jian Cheng, and Hanqing Lu.  Decoupled spatial-temporal attention network for skeleton-based action-gesture recognition.  In Proceedings of the Asian Conference on Computer Vision, 2020.
  * Si et al. [2019] Chenyang Si, Wentao Chen, Wei Wang, Liang Wang, and Tieniu Tan.  An attention enhanced graph convolutional lstm network for skeleton-based action recognition.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1227–1236, 2019.
  * Song et al. [2017] Sijie Song, Cuiling Lan, Junliang Xing, Wenjun Zeng, and Jiaying Liu.  An end-to-end spatio-temporal attention model for human action recognition from skeleton data.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 31, 2017.
  * Song et al. [2022] Yi-Fan Song, Zhang Zhang, Caifeng Shan, and Liang Wang.  Constructing stronger and faster baselines for skeleton-based action recognition.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 45(2):1474–1488, 2022.
  * Wang et al. [2014] Jiang Wang, Xiaohan Nie, Yin Xia, Ying Wu, and Song-Chun Zhu.  Cross-view action modeling, learning and recognition.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 2649–2656, 2014.
  * Xiang et al. [2023] Wangmeng Xiang, Chao Li, Yuxuan Zhou, Biao Wang, and Lei Zhang.  Generative action description prompts for skeleton-based action recognition.  In Proceedings of the IEEE International Conference on Computer Vision, pages 10276–10285, 2023.
  * Xu et al. [2022] Kailin Xu, Fanfan Ye, Qiaoyong Zhong, and Di Xie.  Topology-aware convolutional neural network for efficient skeleton-based action recognition.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pages 2866–2874, 2022.
  * Yan et al. [2018] Sijie Yan, Yuanjun Xiong, and Dahua Lin.  Spatial temporal graph convolutional networks for skeleton-based action recognition.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 32, 2018.
  * Ye et al. [2020] Fanfan Ye, Shiliang Pu, Qiaoyong Zhong, Chao Li, Di Xie, and Huiming Tang.  Dynamic gcn: Context-enriched topology learning for skeleton-based action recognition.  In Proceedings of the ACM International Conference on Multimedia, pages 55–63, 2020.
  * Zhang et al. [2017] Pengfei Zhang, Cuiling Lan, Junliang Xing, Wenjun Zeng, Jianru Xue, and Nanning Zheng.  View adaptive recurrent neural networks for high performance human action recognition from skeleton data.  In Proceedings of the IEEE International Conference on Computer Vision, pages 2117–2126, 2017.
  * Zhang et al. [2019] Hong-Bo Zhang, Yi-Xiang Zhang, Bineng Zhong, Qing Lei, Lijie Yang, Ji-Xiang Du, and Duan-Sheng Chen.  A comprehensive survey of vision-based human action recognition methods.  Sensors, 19(5):1005, 2019.
  * Zhang et al. [2020] Pengfei Zhang, Cuiling Lan, Wenjun Zeng, Junliang Xing, Jianru Xue, and Nanning Zheng.  Semantics-guided neural networks for efficient skeleton-based human action recognition.  In proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1112–1121, 2020.
  * Zhou et al. [2022] Yuxuan Zhou, Zhi-Qi Cheng, Chao Li, Yanwen Fang, Yifeng Geng, Xuansong Xie, and Margret Keuper.  Hypergraph transformer for skeleton-based action recognition.  arXiv preprint arXiv:2211.09590, 2022.
  * Zhou et al. [2023] Huanyu Zhou, Qingjie Liu, and Yunhong Wang.  Learning discriminative representations for skeleton-based action recognition.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 10608–10617, 2023.


# Highly nonlinear dipolar exciton-polaritons in bilayer MoS2

Biswajit Datta  Department of Physics, City College of New York, New York, NY,
USA<EMAIL_ADDRESS><EMAIL_ADDRESS>Mandeep Khatoniar  Department
of Physics, City College of New York, New York, NY, USA  Department of
Physics, Graduate Center of the City University of New York (CUNY), New York,
NY, USA  Prathmesh Deshmukh  Department of Physics, City College of New York,
New York, NY, USA  Department of Physics, Graduate Center of the City
University of New York (CUNY), New York, NY, USA  Félix Thouin  Department of
Engineering Physics, École Polytechnique de Montréal, Montréal, Quebec, Canada
Rezlind Bushati  Department of Physics, City College of New York, New York,
NY, USA  Department of Physics, Graduate Center of the City University of New
York (CUNY), New York, NY, USA  Simone De Liberato  School of Physics and
Astronomy, University of Southampton, Southampton, UK  Stephane Kena Cohen
Department of Engineering Physics, École Polytechnique de Montréal, Montréal,
Quebec, Canada  Vinod M. Menon  Department of Physics, City College of New
York, New York, NY, USA  Department of Physics, Graduate Center of the City
University of New York (CUNY), New York, NY, USA<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Realizing nonlinear optical response in the low photon density limit in solid-
state systems has been a long-standing challenge. Semiconductor microcavities
in the strong coupling regime hosting exciton-polaritons have emerged as
attractive candidates in this context. However, the weak interaction between
these quasiparticles has been a hurdle in this quest. Dipolar excitons provide
an attractive strategy to overcome this limitation but are often hindered by
their weak oscillator strength. The interlayer dipolar excitons in naturally
occurring homobilayer MoS2 alleviates this issue owing to their hybridization
of interlayer charge transfer exciton and intralayer B exciton. Here we
demonstrate the formation of dipolar exciton polaritons in bilayer MoS2
resulting in unprecedented nonlinear interaction strengths. A ten-fold
increase in nonlinearity is observed for the interlayer dipolar excitons
compared to the conventional A excitons. These highly nonlinear dipolar
polaritons will likely be a frontrunner in the quest for solid-state quantum
nonlinear devices.

## Introduction

Photons are becoming ubiquitous in emerging quantum technologies like quantum
communication and metrology due to the ability to propagate long distances
while being robust against decoherence 1, 2. To further extend the span of
their utilities, platforms to achieve and implement strongly interacting
photons in solid-state systems are highly desirable 3. Conventional materials
do not exhibit nonlinear response at power levels associated with single
photons. In this context, remarkable advances have been made in cold atom
systems to realize interacting photons at the single-particle level 4. A
strong contender for the generation of strongly interacting photons in the
solid-state are exciton-polaritons formed via non-perturbative coupling of
cavity photons with excitonic resonances. Although they can be modeled as non-
interacting quantum fluids at low densities, beyond a critical density,
saturation and short-range exchange interactions become significant and give
rise to various phenomena like the appearance of spontaneous coherence,
parametric down-conversion, and superfluidity 5. Strong spatial confinement
along with Coulomb interactions give rise to even stronger correlations
between polaritons that can no longer be described using a mean-field theory.
Such interactions give rise to non-Poissonian statistics of laser
transmission, dubbed as polariton blockade. Preliminary evidence of such non-
classical correlation was recently observed in confined polariton systems in
GaAs 6, 7 but a small ratio of interaction to dissipation rates resulted in
only a weak violation of classical correlations.

Transition metal dichalcogenides (TMDCs) have garnered much attention for
their exceptional optoelectronic properties and have been used to demonstrate
a wide array of fundamental phenomena and technological applications 8.
Furthermore, their two-dimensional (2D) nature allows for the formation of
heterojunctions or homojunctions with arbitrary twist angles resulting in
emergent properties 9, 10, 11, 12. Their strong binding energy and large
oscillator strength render them capable of forming polaritons at room
temperature 13, which has also been shown to retain the intriguing properties
of the 2D excitons 14, 15, 16. Polariton interactions in TMDCs under different
configurations are being studied extensively. Realization of Fermi polarons
17, trion polaritons 18, 19, and excited-state Rydberg exciton-polaritons 20
have shown great potential in their abilities to harness strong polariton
interactions. Moiré exciton-polaritons were recently demonstrated using
heterobilayers of WS2/MoSe2 where the electronic confinement potential arising
from the twisted heterostructure was shown to enhance the nonlinearity 9.
However, all the platforms above rely on short-range exchange interactions or
phase-space filling that pose a bottleneck in realizing few polariton non-
linearity under current experimental capabilities.

In this context, a very attractive possibility is the use of spatially
indirect interlayer excitons (IE) in TMDC heterostructures 21 that can have a
permanent dipole moment and hence support highly interacting dipolar
polaritons. Formation of polaritons with excitons that possess a permanent
dipole moment has been shown to enhance the polariton interactions both in
resonant 22, 23 and non-resonant excitation schemes 24. However, their
inherently low oscillator strength creates an impediment in reaching the
strong coupling regime without hybridization with a direct intralayer exciton.
Furthermore, the large in-built interfacial electric fields make electrical
tuning of the energies of the IE far more difficult. Bilayer MoS2 provides a
highly attractive platform to solve these issues 25, 26, 27, 28. The IE in
naturally stacked 2H bilayer MoS2 remarkably has both an out-of-plane (static)
dipole moment and an in-plane (oscillating) dipole moment. Due to the in-plane
dipole moment, the IE in MoS2 bilayer has an oscillator strength of
approximately 36% of that of the intralayer A exciton along with strong
absorption that is visible up to room temperatures. Moreover, a strong
response to DC electric fields has been demonstrated in these systems, thus
confirming their dipolar nature 25, 26, 28. In this work, we achieve a strong
coupling of microcavity photons with the IEs (along with intralayer A and B
excitons) in bilayer MoS2. The IE polariton shows 10 fold enhancement of the
polariton nonlinearity compared to the intralayer A exciton owing to its
dipolar and phase space filling contributions. Such enhanced nonlinear
response makes them appealing to realize strongly interacting polaritons in
condensed matter systems. In addition, the ease of fabrication and realization
of multi-polariton species in this system makes it a practical and
fundamentally interesting material for studying polariton physics.

## Results

Fig. 1a shows a schematic of the bands that form interlayer exciton (IE) in
bilayer MoS2. IE1 and IE2 excitons in which electrons are in layer 1, and
layer 2 respectively, are energetically degenerate at zero external bias while
the hole is delocalized among both layers. The black arrows in Fig. 1a
indicate the transitions that form the IE. This peculiar charge distribution
is at the heart of producing both in-plane and out-of-plane dipole moments of
IE. Fig. 1b shows the charge distribution in real space for IE1 and IE2. Owing
to the spatial separation of the electron and the hole, it acquires a
permanent dipole moment, the nature of which has been a topic of much interest
recently 25, 27, 28, 26. IEs in MoS2 homobilayer can be thought of as an
admixture of B intralayer exciton with an optically dark but electric field
tunable charge-transfer exciton, which is typically found in TMDC
heterobilayers 29. As a result, IEs in MoS2 homobilayer acquire both a strong
oscillator strength and electric field tunability 29. Fig. 1c shows the white
light differential reflection of the bilayer sample at zero power limit at 7
K. The peaks at 1.9365 eV, 2.000 eV, and 2.1076 eV correspond to the A, IE,
and B excitons, respectively. Interestingly, the absorption of the IE is
prominent even at room temperature, which provides a straightforward method to
identify the MoS2 bilayer flakes after mechanical exfoliation (see SI Fig.S3).

We probe the nonlinear optical response of the IE and A exciton using resonant
pump experiments (see Methods). The relative shift of A exciton and IE so
obtained are plotted in Fig. 1d against the estimated density of A exciton and
IE, respectively (see Methods for the details of the exciton density
calculation). Both the excitons blueshift with increased excitation density,
with the IE shift being more pronounced. Exchange interaction is the dominant
mechanism responsible for the blue shift for the A exciton, while dipolar
interactions dominate the interaction in the IE response as exchange
interaction is likely to be small given the larger spatial separation of
electron and holes 30. To better understand the exciton dynamics, we have also
performed ultrafast pump-probe spectroscopy on the bilayers (see SI section 10
for the pump-probe data). In this case, the spectra show a redshift of the
resonances accompanied by a nearly identical transient response regardless of
which exciton is first pumped (A, IE or B). This can be understood if carriers
rapidly relax out of the K valley to the band minima on a scale faster than
the $\sim$ 200 fs time resolution. This supports our hypothesis that the A and
IE transitions are nearly homogeneously broadened with a lifetime dictated by
the intervalley scattering rate. This agrees well with recent coherent
spectroscopy and numerical calculations performed for bilayers of MoSe2 31.
The above approach is used to estimate the exciton density in Fig. 1d.

Exciton polaritons were realized in the bilayer MoS2 by embedding it in a
microcavity. Fig. 2a shows the schematic of the structure used in this
experiment. It consists of a bilayer MoS2 encapsulated by thin (20 nm)
hexagonal boron nitride (hBN), which is sandwiched between a bottom dielectric
distributed Bragg reflector (DBR) mirror and a top silver mirror. See SI Note
1 and Methods for the fabrication details and cavity structure. The cavity is
designed such that the bilayer MoS2 flake sits close to an electric field
anti-node, allowing us to observe five dispersive modes associated with the
different polaritonic states as shown in Fig. 2b. These distinct polariton
modes arise due to the hybridization of the cavity photon mode with the
various excitonic states present in the bilayer system. We name the polariton
branches as pol-1 through pol-5, with pol-1 corresponding to the lowest energy
and pol-5 the highest. From the Fig. 2b, we see that the 1s state of A, the
IE, B exciton and 2s Rydberg state of A exciton all strongly couple to the
same cavity mode at 7 K. We fit the data with a five-coupled oscillator model
where the energy of all the four excitons, the cavity mode, their Rabi
splitting, and the effective refractive index of the system are treated as fit
parameters. The fit results in exciton energies $E_{A1s}=1.9323\pm 0.0003$ eV,
$E_{IE}=2.0014\pm 0.0001$ eV, $E_{A2s}=2.078\pm 0.002$ eV, $E_{B}=2.111\pm
0.001$ eV and the Rabi splittings $\Omega_{A1s}=40.4\pm 0.3$ meV,
$\Omega_{IL}=21.4\pm 0.1$ meV, $\Omega_{A2s}=13\pm 0.5$ meV, $\Omega_{B}=51\pm
0.4$ meV. The exciton energies obtained through the fit agree well with the
experimentally determined exciton energies with slight shifts due to the
strain and changes in the dielectric environment in the cavity. In our device,
the bare cavity mode is positively detuned from the A1s exciton by
$\delta_{C-A_{X}}$ = 17 meV. See SI Note 2 for details on the coupled
oscillator model. Fig. 2c shows the line cuts of the full k-space data at
three different k∥ where the cavity mode is resonant with the A exciton
(yellow), IE (red) and B exciton (blue), respectively. We notice that the Rabi
splitting of the IE polariton remains observable even at 77 K (see SI Fig.S3).

We measured the white light reflectivity at the IE and A resonance and monitor
the energy shifts of the polariton branches, see SI Note 5 and SI Note 8 for
the details of the experiment. Fig. 3a shows the energy of the upper and lower
branch of the IE polariton as a function of polariton density at the inplane
wavevector corresponding to the zero-detuning $k_{\parallel}$ (where the
cavity mode and exciton energy are degenerate). Two Lorentzian fits at each
density is used to obtain the energy of the upper and lower polaritons. SI
Note 4 discusses the details of the polariton density calculation and SI Note
5 provides the raw data of the density dependent differential reflection of IE
polariton. We can see from Fig. 3 that the lower branch of the IE polariton
moves with pump power more than the upper branch. We also observe a
simultaneous increase of the the magnitude of zero-detuning k|| along with the
reduction of Rabi splitting. This suggests the presence of both exciton energy
renormalization and saturation effect. SI Note 7 shows the density-dependent
Rabi splitting of IE polariton. Taking the derivative of the ELP with respect
to density, we calculate the strength of the nonlinearity, gLP. Since the rate
of blue shift saturates at high density (olive curve in Fig. 3a),
g${}^{IE}_{LP}$ reduces with polariton density as seen in Fig. 3b. Note that
at the lowest powers accessible in our measurement, the IE lower polariton
branch already shows a power-dependent blueshift. On the other hand, within
the accessible range of the polariton density in our experiment,
g${}^{A}_{LP}$ remains nearly a constant (see SI Note 8 for the blueshift data
of the lower polariton branch of A exciton). The obtained strength of the non-
linearity in the low density limit for IE and A exciton polaritons are
g$\mathrm{{}^{IE}}_{\mathrm{LP}}$ $\sim(100\pm 2)$ µeV µm2 and
g$\mathrm{{}^{A}}_{\mathrm{LP}}$ $\sim(10\pm 0.2)$ µeV µm2 respectively
corresponding to 10 fold increase in the nonlinear response for the IE
polaritons.

## Discussion

We note that the g$\mathrm{{}^{A}}_{\mathrm{LP}}$ value is comparable to the
previous report for A exciton-polariton in monolayer TMDCs 20, 32. This allows
us to compare the nonlinearity of the IE in bilayer MoS2 with that of Rydberg
A2s state reported in ref.20 – we find that the nonlinear response of IE is
even larger than the A2s state by a factor of $\sim$2.7. More importantly, we
would like to highlight a few critical differences between these two studies.
The approximately 4 times enhancement in the gLP for A2s state in ref.20 is
attributed to the enhanced phase-space filling due to a higher Bohr radius of
2s compared to 1s. However, this enhancement in nonlinearity comes at the cost
of reduced oscillator strength. This is overcome in ref.20 by stacking three
layers of monolayer TMDC separated by thin hBN. There is, however, no clear
experimental route towards using higher-order Rydberg states beyond $n>2s$ in
TMDCs to get even higher enhancement as the oscillator strength almost
vanishes for higher-order Rydberg states. Additionally, the higher-order
Rydberg states are also more susceptible to disorder due to their larger Bohr
radius. On the other hand, the spectral signature of interlayer exciton in
bilayer MoS2 is visible even at room temperature and was shown to have a
significant response to applied out of plane electric field 25. This provides
a platform with the possibility of creating highly interacting dipolar
polaritons controlled by external electric fields.

Nonlinear polariton interactions in the present case have contributions from
exciton-exciton Coulomb interaction (which includes both dipolar and exchange
terms) and phase space filling effect. For intralayer excitons with a small
Bohr radius only the exchange term and phase space filling effect is important
30. On the other hand, for interlayer excitons in traditional quantum well
systems, the direct dipole-dipole interaction dominates, and the exchange
interaction is neglected due to the separation between charges in the two
layers 33. In contrast in MoS2 homobilayers, the hole is delocalized across
the two layers, and hence there is a finite overlap with the electron
wavefunction leading to small but non-zero exchange interactions. Thus in the
present experiments, in addition to the saturation effects, both dipolar and
exchange interaction contribute to the overall nonlinear response, with the
latter being less significant owing to the finite separation between the
electron and hole.

To discern the contributions to the polaritonic nonlinearity due to dipolar
and exchange interaction of excitons ($\Delta\mathrm{E_{exc-exc}}$) versus
phase space filling ($\Delta\mathrm{E_{sat}}$) we carry out the following
analysis. Fig. 4a shows the Hopfield coefficients of all the excitonic
components and the cavity photon as a function of $k_{\parallel}$. This is
used to obtain the shift in polariton energy, $\Delta
E^{\mathrm{IE}}_{\mathrm{LP}}$ as a function of Hopfield cavity photon
fraction ($C$), see Fig. 4b. The interaction strength due to phase space
filling can be written in terms of
$g_{\mathrm{LP}}=4g_{\mathrm{SAT}}|C||X|^{3}$, where X is the Hopfield
coefficient for the exciton. Although this is a multi-exciton system, the
contribution of the other excitons close to the zero-detuning $k_{\parallel}$
of the IE is small. Thus, we can write
$g^{IE}_{\mathrm{LP}}=4g_{\mathrm{SAT}}|C|(1-|C|^{2})^{3/2}$ which is a non-
monotonic function of the cavity Hopfield coefficient as shown in the Fig. 4b
inset (red curve). For the case of exchange and dipolar interaction, the
interaction strength scales as
$g_{\mathrm{LP}}=g_{\mathrm{XX}}|X|^{4}=g_{\mathrm{XX}}(1-|C|^{2})^{2}$ 23, 33
which monotonically decreases with increasing $|C|$ (cyan curve in Fig. 4b
inset). We find that the measured nonlinear response of the IE polariton is
the highest at an in-plane wave vector larger than zero-detuning
$k_{\parallel}$, as seen in Fig. 4b. We also notice that $\Delta
E^{\mathrm{IE}}_{\mathrm{LP}}$ changes non-monotonically with the cavity
Hopfield coefficient. The tilted parabolic shape of the measured $\Delta
E^{\mathrm{IE}}_{\mathrm{LP}}$ as a function of $|C|$ indicates a roughly
equal contributions from both $\Delta\mathrm{E_{exc-exc}}$ and
$\Delta\mathrm{E_{sat}}$.

Fig. 4c shows a schematic of the polariton energies at two different powers
for $\Delta\mathrm{E_{sat}}>\Delta\mathrm{E_{exc-exc}}$. Due to both the blue
shift of the exciton and saturation effect, the lower branch moves with
increasing power more than the upper branch. Analyzing the blueshift of the IE
lower polariton branch and redshift of the IE upper polariton branch at the
zero-detuning $k_{\parallel}$ we calculate the magnitude of
$\Delta\mathrm{E_{exc-exc}}$ and $\Delta\mathrm{E_{sat}}$ in our system. At
the zero-detuning $k_{\parallel}$ we write the energy shift of IE lower branch
as $\Delta E^{\mathrm{IE}}_{\mathrm{LP}}$ = $\Delta\mathrm{E_{exc-exc}}$ \+
$\Delta\mathrm{E_{sat}}$ and the energy shift of IE upper branch as $\Delta
E^{\mathrm{IE}}_{\mathrm{UP}}$ = $\Delta\mathrm{E_{exc-
exc}}-\Delta\mathrm{E_{sat}}$. From these equations we obtain
$\Delta\mathrm{E_{exc-exc}}$ and $\Delta\mathrm{E_{sat}}$ as a function of
polariton density, see SI Note 9. The ratio of $\Delta$Esat/$\Delta$Eexc-exc
as a function of polariton density is shown in Fig. 4d, indicating that the
magnitude of the dipole mediated exciton-exciton interaction nonlinearity and
saturation nonlinearity are comparable especially at low polariton densities.
Further increase of dipole-dipole interaction by creating imbalance between up
and down dipoles (via applied electric field or providing strain) will allow
to access regime where dipolar interactions dominate.

We now discuss the origin of the net dipolar interaction arising in our system
at the zero electric field. The two ground states of the interlayer exciton
are quasi degenerate at the zero electric field and should not possess net
dipole moment in an ideal limit. However, as our experiments indicate, dipolar
interactions play a crucial role in the 10 fold enhancement in interaction
strength we observe for the interlayer exciton-polaritons. This raises
questions about the origin of the dipole moment. Structural asymmetry, strain,
residual electric fields, and disorder are the likely reasons for breaking the
symmetry and resulting in a net dipole moment. Given that the bilayer MoS2 is
sandwiched between two hBN layers of similar thickness, one can rule out the
role of structural asymmetry. While local strain and built-in fields could
have highly localized effects, we observed similar enhancement over the entire
bilayer sample, and hence one can conclude that local strain or fields are
also not the key contributors. Finally, if one looks at disorder arising from
defects present in the 2D TMDCs, the typical defect densities observed in
mechanically exfoliated monolayer MoS2 is on the order of $10^{12}$-$10^{13}$
$\mathrm{cm}^{-2}$ as has been shown before 34. Using the lower bound of
$10^{12}$ $\mathrm{cm}^{-2}$ we can estimate the lower bound for dipolar
interaction arising from local symmetry breaking between the two layers of the
bilayer system. Using a Poisson distribution of defects, we estimate 36% of
the carriers in one layer see a defect in its vicinity, thereby introducing a
net dipole moment. This effect has been hypothesized to be the reason behind
the observation of a sizable second harmonic generation in bilayer MoS2 even
at zero perpendicular DC electric field 35. This local imbalance of the up and
down dipole moment will increase the dipole-dipole interaction observed in our
experiment.

In summary, we demonstrate extremely large optical nonlinearity arising from
IE polaritons in homobilayer MoS2. The reported nonlinearity is an order of
magnitude larger than that reported for the A excitons in TMDCs, bringing the
IE polaritons closer to the polariton blockade regime. The enhanced nonlinear
response of the exciton-polaritons is attributed to the dipolar interaction of
the interlayer excitons and the saturation nonlinearity with the former
starting to dominate in the weak pump limit. The homobilayer MoS2 system
presents itself as an attractive platform to realize IE polaritons without the
need for precise twist angle control while having the potential to reach
elevated operational temperatures. Using an electric field to induce
preferential dipole moment along with lateral photon confinement in the
microcavity is expected to push the bilayer MoS2 polariton system into the
quantum nonlinear regime.

## Methods

## Fabrication details

The DBR used in our experiment is made of 8 pairs of SiO2 (104.9 nm) and TiO2
(64.7 nm) layers and were deposited by radio frequency sputtering on an
intrinsic Si chip. MoS2 and hBN were exfoliated from bulk crystals (from
2Dsemiconductor Inc.) using blue tape (Nitto) and scotch tape respectively.
MoS2 was exfoliated onto PDMS substrate and hBN was exfoliated onto a 300 nm
SiO2/Si substrate. Bilayer MoS2 flakes were identified by reflection
spectroscopy – the additional dip at $\sim$ 639 nm in the reflection signal at
room temperature is the hallmark of bilayer MoS2. hBN/bilayer MoS2/hBN
heterostructure stacking and transfer were done using the well-known
polypropylene carbonate transfer technique 36. The chosen top and bottom hBN
layers were of similar thickness $\sim$ 20 nm. The final stack was then
transferred onto the DBR at a temperature of 150°C. The chip was kept in
chloroform for 12 hours to remove the polypropylene carbonate residue. 240 nm
polymethylmethacrylate (PMMA) (495 A4 from Michrochem) was spin-coated to form
a 3$\lambda$/2 cavity. Finally, a silver layer (40 nm) was deposited by
electron-beam evaporation for the top mirror of the microcavity. Optical
images of the device and the bare cavity response are given in SI Note 1.

## Optical measurement details

We recorded the angle-resolved reflection spectra using the Fourier space
imaging technique. A broadband halogen light source was used for reflection
measurements. A supercontinuum pulsed light source (NKT Photonics, repetition
rate 80 MHz, pulse duration 20 ps) was used to study the power dependence of
the bare excitons and polaritons. Appropriate long pass and short pass filters
were used in the input to excite a narrow band around the excitons/polaritons
under study. For the IE polariton, the excitation bandwidth of the
supercontinuum laser was chosen to be $\sim$40 meV around the IE energy in the
power-dependent measurements. The set-up was coupled with a Princeton
Instruments monochromator with a CCD camera. A 50X, 0.65 numerical aperture
objective was used for all the measurements at 7 K in a closed cycle Montana
cryostat. The polariton dispersion was revealed by imaging the back aperture
of the microscope objective (Fourier plane) onto the camera. The spot size of
the laser was 1 µm2 on a uniform area of the sample of dimension $\sim$500
µm2.

## Exciton density calculation in the absence of the top mirror

For the data presented in Fig.1d, we excite A, IE and B exciton simultaneously
at various optical fluences. We have used a supercontinuum laser with a 20 ps
pulse width and broadband excitation (1.907 eV to 2.246 eV) scheme for this
measurement. At all fluences, each spectrum is fitted with the sum of three
Lorentzians to track the evolution of each feature with excitation fluence
(see SI Note 6). To calculate the density of individual excitons, we do the
following.

We first calculate the absorbed power density by the bilayer MoS2

$P_{\mathrm{absorbed}}=\frac{\int_{E_{min}}^{E_{max}}P_{laser}(E)L(E-E_{X})dE}{A_{beam}}$
(1)

Here $P_{laser}(E)$ is the intensity profile of the incident supercontinuum
laser, $L(E-E_{X})$ is the absorbance of X exciton around its energy $E_{X}$,
and $A_{beam}$ is the size of the focused laser beam on the sample. We run the
integration within 74 meV (53 meV) band around the A exciton (IE) peak energy,
respectively. Beyond these bands, the reflection contrast of A exciton and IE
are not detectable. The absorbed photon density from the absorbed power
density is calculated taking into account the life time of these excitons (45
fs for A exciton 53 fs for IE). We assume all of the absorbed photons form
excitons.

## Transient reflectivity measurements

Transient reflection spectroscopy measurements presented in SI Note 10 were
obtained using an ultrafast transient absorption spectrometer (Harpia, Light
Conversion). An ultrafast optical parametric oscillator (Orpheus-F, Light
Converison) provided the pump pulses and white light continuum generation in a
sapphire crystal provided the probe pulse. These were filtered using a pinhole
inside a telescope and collinearly sent to a home-built microscope equipped
with a 50x objective (NA 0.45). The sample was held at 3.5 K using a closed-
cycle cryostat (Montana Instruments) overlooked by the microscope. The output
lens of the filtering telescope was adjusted to obtain pump and probe spots as
large as the bilayer flake. Before its detection by a spectrometer, the pump
was removed from the light reflected by the sample using a polarizer so that
only the probe beam was detected.

## Error in estimating density of polaritons

Currently, the density error bars consider the error in measuring power,
determining the energy and linewidth of the polariton. However, one primary
source of the error is the lifetime of the polaritons. To address this issue
we have done pump-probe experiment which validates the use of the linewidth of
the polaritons to determine the polariton lifetime. Results of our pump-probe
experiment is included in the SI note 10. Another likely source of error is
estimating efficient photons’ loading into a cavity. This is because the
conventional input-output theory does not account for lossy mirrors, as in our
experiment. Also photon loading into microcavities from free space is not
efficient due to mode mismatch between the free space laser modes and bi-
orthonormal modes in lossy cavities 37. However, the relative enhancement of
the nonlinearity of IE with respect to A remains unaffected due to the above
errors as they exist for both the IE and A exciton polariton.

## Acknowledgements

We acknowledge support from ARO through grants W911NF-17-1-0312 and
W911NF-22-1-0091; NSF DMR-2011738, and AFOSR-FA2386-21-1-4087. The low
temperature spectroscopy facility was supported by NSF MRI grant -
DMR-1726573. S.D.L. acknowledges support from a Royal Society Research
Fellowship and from the Philip Leverhulme Prize of the Leverhulme Trust.

## Contributions

V.M.M., M.K. and B.D. conceived the experiments. B.D. fabricated the device.
R.B. helped in the fabrication. B.D., M.K. and P.D. performed the
measurements. F.T. did the time resolved pump probe measurements. B.D. did the
data analysis with inputs from M.K., P.D. and V.M.M. S.D.L. and S.K.C. gave
inputs in the data interpretation. B.D, M.K. and V.M.M wrote the manuscript
with comments from all authors.

Figure 1:  Excitons in MoS2 bilayer and white light absorption. (a) Schematic
of the band structure (around the K point) of the bilayer MoS2 showing (black
arrows) the participation of the different bands in the IE formation. The blue
and red colors represent the up and down spin, respectively. L1 and L2 refer
to the layer 1 and layer 2 of the bilayer MoS2, respectively. (b) Schematic
showing the two possible charge configurations in the absence of bias voltage
between the two layers. In both cases, the electron is localized in one of the
layers, and the hole is distributed among both the layers. (c) White-light
differential reflection at 7 K showing three dips corresponding to A, IE, and
B exciton, respectively. (d) Relative shift of excitonic features observed in
instantaneously-probed pulsed reflectivity with increasing exciton density.
Pulse duration = 20 ps. The K valley lifetime of A exciton (45 fs) and IE (53
fs) estimated from their respective linewidths are used to estimate the
exciton density from the pump fluence. The error bars represent the variance
of the fitted parameters caused by noise in the measurement.

Figure 2:  Device schematic and polariton dispersion. (a) Schematic of the
device. The bilayer MoS2 is sandwiched between two hexagonal boron nitride
flakes of almost similar thickness $\sim$20 nm. The entire stack is
transferred by the polypropylene carbonate (PPC) method onto a SiO2 terminated
SiO2-TiO2 DBR. As the cavity spacer, a layer of 240 nm thick PMMA was spin-
coated on the stack. Finally, 40 nm of Ag was deposited via electron beam
evaporation to form a $3\lambda/2$ cavity, where $\lambda$ is the cavity
resonance wavelength. (b) White light differential reflection showing all
branches of the polaritons at 7 K. 1s and 2s states of A exciton, IE, and B
exciton all are strongly coupled to the same cavity mode. This data was
measured by illuminating the sample with incoherent halogen light with a power
of few nW at the sample. The yellow dashed lines are the polariton eigenvalues
obtained from fit using the five coupled oscillator model. The brown dashed
lines denote the energy of the excitons and the cavity mode resulting from the
fitting. Here $k_{0}=\frac{2\pi}{\lambda_{c}}$; $\lambda_{c}$ is the
wavelength of the uncoupled cavity mode at $k_{\parallel}=0$. (c) Line cuts at
different $k_{\parallel}$ of the color plot shown in panel-b. It clearly shows
that each exciton strongly couples, resulting in two polariton modes that lie
above and below the excitonic resonances.

Figure 3:  Comparison of nonlinearity of the IE polariton with A exciton
polariton. (a) Energy of the lower and upper branch of the IE polariton as a
function of the polariton density at zero-detuning $k_{\parallel}$. The
movement of the lower branch is more than the upper branch due to the combined
effect of exciton-exciton interaction and saturation. The error bars in energy
represent the uncertainty in determining the peak of the Lorentzian fit
(polariton energy) to the reflection data. (b)
g${}^{\mathrm{IE}}_{\mathrm{LP}}$ (g${}^{\mathrm{A}}_{\mathrm{LP}}$) as a
function of the IE (A) lower branch polariton density at zero-detuning
$k_{\parallel}$. gLP is defined as the local slope of the polariton energy vs.
polariton density graph. g${}^{\mathrm{A}}_{\mathrm{LP}}$ and
g${}^{\mathrm{IE}}_{\mathrm{LP}}$ are calculated from the lower branch of A
polariton and IE polariton, respectively. Note that the gLP of IE polariton
$\sim(100\pm 2)$ µeV µm2 is almost 10 times larger than the gLP of A exciton
polariton $\sim(10\pm 0.2)$ µeV µm2. All the power-dependent nonlinear
measurements were carried out using a pulsed supercontinuum laser (20 ps
pulsewidth) with proper bandpass filter in the input to excite only one
polariton species. The polariton density corresponds to the specific branch
being excited. The density error bars take into account the error in measuring
power and determining the energy of the polariton and linewidth from the
Lorentzian fits. The error bars in gLP consider the error in determining the
density and the energy of the polariton.

Figure 4:  Nonlinear response of the IE polariton. (a) The Hopfield
coefficients of the lower IE polariton showing the contributions of cavity
photon and the four excitons as a function of $k_{\parallel}$. (b)
Experimentally measured blueshift of the IE lower branch as a function of
$k_{\parallel}$ at a polariton density $\sim 95$ µm-2. Here $k_{\parallel}$ is
converted to cavity Hopfield coefficient using the previous plot. The
blueshift of the lower branch of IE ($\Delta E^{\mathrm{IE}}_{\mathrm{LP}}$)
is highest close to the zero-detuning $k_{\parallel}$, but the shape looks
like a tilted parabola. This can be explained by taking both exciton-exciton
interaction and saturation effect into account. Inset shows the theoretically
expected dependence of exciton-exciton interaction and saturation effect on
$|C|$. It can be seen that superimposition of these two graphs can explain the
tilted parabola shape of $\Delta E^{\mathrm{IE}}_{\mathrm{LP}}$ in our data.
(c) Cartoon showing the IE polariton dispersion at two different excitation
powers where the combined effect of exciton-exciton interaction and saturation
effect is present. Here, X1 and X2 are the exciton energies at power P1 and
P2, respectively. The zero-detuning $k_{\parallel}$ increases with increasing
power since the exciton blueshifts. Due to the combined effect of exciton
blueshift and saturation, the lower IE polariton branch moves more than upper
IE polariton branch. (d) Ratio of estimated saturation nonlinearity and
exciton-exciton interaction nonlinearity as a function the IE lower branch
polariton density. The error bars in energy represent the uncertainty in
determining the peak of the Lorentzian fit (polariton energy) to the
reflection data.

## References

  * 1 Gisin, N. & Thew, R.  Quantum communication.  _Nature photonics_ 1, 165–171 (2007).
  * 2 Giovannetti, V., Lloyd, S. & Maccone, L.  Quantum-enhanced measurements: beating the standard quantum limit.  _Science_ 306, 1330–1336 (2004).
  * 3 Chang, D. E., Vuletić, V. & Lukin, M. D.  Quantum nonlinear optics—photon by photon.  _Nature Photonics_ 8, 685–694 (2014).
  * 4 Lukin, M. D. _et al._ Dipole blockade and quantum information processing in mesoscopic atomic ensembles.  _Phys. Rev. Lett._ 87, 037901 (2001).  URL https://link.aps.org/doi/10.1103/PhysRevLett.87.037901.
  * 5 Carusotto, I. & Ciuti, C.  Quantum fluids of light.  _Reviews of Modern Physics_ 85, 299–366 (2013).  1205.6500.
  * 6 Muñoz-Matutano, G. _et al._ Emergence of quantum correlations from interacting fibre-cavity polaritons.  _Nature Materials_ 18, 213–218 (2019).  URL https://doi.org/10.1038/s41563-019-0281-z.
  * 7 Delteil, A. _et al._ Towards polariton blockade of confined exciton–polaritons.  _Nature Materials_ 18, 219–222 (2019).  URL https://doi.org/10.1038/s41563-019-0282-y.
  * 8 Mak, K. F. & Shan, J.  Photonics and optoelectronics of 2d semiconductor transition metal dichalcogenides.  _Nature Photonics_ 10, 216–226 (2016).
  * 9 Zhang, L. _et al._ Van der waals heterostructure polaritons with moiré-induced nonlinearity.  _Nature_ 591, 61–65 (2021).
  * 10 Alexeev, E. M. _et al._ Resonantly hybridized excitons in moiré superlattices in van der waals heterostructures.  _Nature_ 567, 81–86 (2019).
  * 11 Zhang, L. _et al._ Twist-angle dependence of moiré excitons in ws 2/mose 2 heterobilayers.  _Nature communications_ 11, 1–8 (2020).
  * 12 Tang, Y. _et al._ Tuning layer-hybridized moiré excitons by the quantum-confined stark effect.  _Nature Nanotechnology_ 16, 52–57 (2021).
  * 13 Liu, X. _et al._ Strong light-matter coupling in two-dimensional atomic crystals.  _Nature Photonics_ 9, 30–34 (2014).  URL http://dx.doi.org/10.1038/nphoton.2014.304.  1406.4826.
  * 14 Sun, Z. _et al._ Optical control of room-temperature valley polaritons.  _Nature Photonics_ 11, 491–496 (2017).  URL https://doi.org/10.1038/nphoton.2017.121.
  * 15 Chen, Y.-J., Cain, J. D., Stanev, T. K., Dravid, V. P. & Stern, N. P.  Valley-polarized exciton–polaritons in a monolayer semiconductor.  _Nature Photonics_ 11, 431–435 (2017).  URL https://doi.org/10.1038/nphoton.2017.86.
  * 16 Dufferwiel, S. _et al._ Valley-addressable polaritons in atomically thin semiconductors.  _Nature Photonics_ 11, 497–501 (2017).  URL https://doi.org/10.1038/nphoton.2017.125.
  * 17 Tan, L. B. _et al._ Interacting polaron-polaritons.  _Phys. Rev. X_ 10, 021011 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevX.10.021011.
  * 18 Emmanuele, R. P. A. _et al._ Highly nonlinear trion-polaritons in a monolayer semiconductor.  _Nature Communications_ 11, 3589 (2020).  URL https://doi.org/10.1038/s41467-020-17340-z.
  * 19 Kravtsov, V. _et al._ Nonlinear polaritons in a monolayer semiconductor coupled to optical bound states in the continuum.  _Light: Science & Applications_ 9, 1–8 (2020).
  * 20 Gu, J. _et al._ Enhanced nonlinear interaction of polaritons via excitonic rydberg states in monolayer wse2.  _Nature Communications_ 12, 2269 (2021).  URL https://doi.org/10.1038/s41467-021-22537-x.
  * 21 Rivera, P. _et al._ Observation of long-lived interlayer excitons in monolayer mose 2–wse 2 heterostructures.  _Nature communications_ 6, 1–6 (2015).
  * 22 Suárez-Forero, D. G. _et al._ Enhancement of parametric effects in polariton waveguides induced by dipolar interactions.  _Phys. Rev. Lett._ 126, 137401 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevLett.126.137401.
  * 23 Togan, E., Lim, H.-T., Faelt, S., Wegscheider, W. & Imamoglu, A.  Enhanced interactions between dipolar polaritons.  _Phys. Rev. Lett._ 121, 227402 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevLett.121.227402.
  * 24 Rosenberg, I. _et al._ Strongly interacting dipolar-polaritons.  _Science Advances_ 4 (2018).  URL https://advances.sciencemag.org/content/4/10/eaat8880.  https://advances.sciencemag.org/content/4/10/eaat8880.full.pdf.
  * 25 Leisgang, N. _et al._ Giant stark splitting of an exciton in bilayer mos 2\.  _Nature Nanotechnology_ 15, 901–907 (2020).
  * 26 Lorchat, E. _et al._ Excitons in bilayer mos 2 displaying a colossal electric field splitting and tunable magnetic response.  _Physical Review Letters_ 126, 037401 (2021).
  * 27 Gerber, I. C. _et al._ Interlayer excitons in bilayer mos 2 with strong oscillator strength up to room temperature.  _Physical Review B_ 99, 035443 (2019).
  * 28 Peimyoo, N. _et al._ Electrical tuning of optically active interlayer excitons in bilayer mos2.  _Nature Nanotechnology_ 1–6 (2021).
  * 29 Deilmann, T. & Thygesen, K. S.  Interlayer excitons with large optical amplitudes in layered van der waals materials.  _Nano letters_ 18, 2984–2989 (2018).
  * 30 Erkensten, D., Brem, S. & Malic, E.  Exciton-exciton interaction in transition metal dichalcogenide monolayers and van der waals heterostructures.  _Physical Review B_ 103, 045426 (2021).
  * 31 Helmrich, S. _et al._ Phonon-assisted intervalley scattering determines ultrafast exciton dynamics in ${\mathrm{mose}}_{2}$ bilayers.  _Phys. Rev. Lett._ 127, 157403 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevLett.127.157403.
  * 32 Stepanov, P. _et al._ Exciton-exciton interaction beyond the hydrogenic picture in a mose 2 monolayer in the strong light-matter coupling regime.  _Physical Review Letters_ 126, 167401 (2021).
  * 33 Byrnes, T., Recher, P. & Yamamoto, Y.  Mott transitions of exciton polaritons and indirect excitons in a periodic potential.  _Physical Review B_ 81, 205312 (2010).
  * 34 Hong, J. _et al._ Exploring atomic defects in molybdenum disulphide monolayers.  _Nature communications_ 6, 1–8 (2015).
  * 35 Shree, S. _et al._ Interlayer exciton mediated second harmonic generation in bilayer mos2.  _Nature communications_ 12, 1–9 (2021).
  * 36 Wang, L. _et al._ One-dimensional electrical contact to a two-dimensional material.  _Science_ 342, 614–617 (2013).
  * 37 Carusotto, I. & Ciuti, C.  Quantum fluids of light.  _Reviews of Modern Physics_ 85, 299 (2013).

Supplementary Materials: Highly nonlinear dipolar exciton-polaritons in
bilayer MoS2

##  I Device fabrication

Figure S1:  Optical micrograph of the device. (a) Optical micrograph of the
hBN/MoS2 bilayer/hBN heterostructure on a distributed Bragg reflector (DBR)
before top silver mirror deposition. The dashed line marks the uniform area of
the heterostructure where measurements are done. Figure S2:  PMMA bare cavity
mode. (a) Dispersion of the PMMA bare cavity mode measured away from the
hBN/bilayer MoS2/hBN heterostructure. The bare cavity mode gets 0.23 eV red
shifted on the sample area due to the hBN and bilayer MoS2. (b) Line cut of
the bare cavity mode at $k_{||}=0$. The red curve is the fitted Lorentzian
with a FWHM of 15.5 meV.

##  II Details of coupled oscillator model

The five polariton branches in our experiment can be described with the eigen
modes of the following five coupled oscillator model.

$\left(\begin{array}[]{ccccc}\text{$E_{cav}$}&\text{$\Omega_{A}/2$}&\text{$\Omega_{IE}/2$}&\text{$\Omega_{A2s}/2$}&\text{$\Omega_{B}/2$}\\\
\text{$\Omega_{A}/2$}&\text{$E_{A}$}&0&0&0\\\
\text{$\Omega_{IE}/2$}&0&\text{$E_{IL}$}&0&0\\\
\text{$\Omega_{A2s}/2$}&0&0&\text{$E_{A2s}$}&0\\\
\text{$\Omega_{B}/2$}&0&0&0&\text{$E_{B}$}\end{array}\right)$

where
$E_{\mathrm{cav}}=\frac{c\hbar}{e\mathrm{n_{c}}}\sqrt{\mathrm{k_{v}}^{2}+k_{||}^{2}}=\frac{c\hbar}{e\mathrm{n_{c}}}\sqrt{\mathrm{k_{v}}^{2}+\left(10^{7}\frac{k_{||}}{k_{0}}\right)^{2}}$,
$k_{v}=\frac{2\pi}{\lambda_{c}}n_{c}$, $k_{0}=\frac{2\pi}{\lambda_{c}}\approx
10^{7}$ nm-1, e is the electronic charge. Here the wavelength is in meter and
energy is in eV. The subscripts cav, A, IE, A2s, and B refer to the cavity
mode, A exciton, imterlayer exciton, 2s state of A exciton, and B exciton
respectively. $\Omega_{X}$ and $E_{X}$ refer to the Rabi splitting and energy
of the X exciton respectively.

The fitting yields $\lambda_{c}=637.5\pm 0.1$ nm, $n_{c}=1.450\pm 0.003$,
$E_{A1s}=1.9323\pm 0.0003$ eV, $E_{IL}=2.0014\pm 0.0001$ eV, $E_{A2s}=2.078\pm
0.002$ eV, $E_{B}=2.111\pm 0.001$ eV and the Rabi splittings
$\Omega_{A1s}=40.4\pm 0.3$ meV, $\Omega_{IL}=21.4\pm 0.1$ meV,
$\Omega_{A2s}=13\pm 0.5$ meV, $\Omega_{B}=51\pm 0.4$ meV.

##  III Reflection measurements on bare exciton and exciton-polariton data at
elevated temperature

Figure S3:  White light reflection and exciton-polariton at elevated
temperatures. (a) The red (black) color plot shows the differential reflection
of bilayer (monolayer) MoS2 on PDMS substrate. The small peak in between the A
and B exciton in bilayer is the signature of the IE. (b) Strong coupling of
the IE and A exciton at 77 K in bilayer MoS2.

##  IV Details of the polariton density calculation

Polariton density is calculated from the Gross-Pitaevskii (GP) equation. The
GP equation at wavevector k can be written as

$i\hbar\frac{\partial\psi_{LP}(k,t)}{\partial
t}=\left[\epsilon_{LP}(k)-\frac{i\hbar\gamma_{LP}}{2}\right]\psi_{LP}(k,t)+\hbar
F_{P}(k,t)$ (S1)

Here $\epsilon_{LP}(k)=\hbar\omega_{LP}(k)$ is the energy of the lower
polariton of an exciton and $\gamma_{LP}$ is the line width of the polariton.
$F_{P}(k,t)$ is the coherent driving. The pulse width of our supercontinuum
laser is 20 ps which is much larger than the polariton lifetime (30-50 fs).
This allows the pumping to be treated as continuous wave which is resonant
with our cavity. We can express the pumping term as
$F_{P}(k,t)=F_{P}(k)e^{-i\omega_{P}t}$ and
$\psi_{LP}(k,t)=\tilde{\psi}_{LP}(k)e^{-i\omega_{P}t}$ where $\omega_{P}$ is
the pump frequency. This simplifies the above equation to

$\left[\hbar\omega_{P}-\epsilon_{LP}(k)+\frac{i\hbar\gamma_{LP}}{2}\right]\tilde{\psi}_{LP}(k)=\hbar
F_{P}(k)$ (S2)

The $F_{P}(k)$ term can be written from the input-output relation.

$F_{P}(k)=C(k)\sqrt{\frac{\eta\cdot P_{\text{int }}(k)}{\hbar\omega_{P}}}$
(S3)

$\eta=\frac{\left|t_{\text{top mirror}}\right|^{2}}{\tau_{\text{trip }}}$ (S4)

Here C(k) is the photon Hopfield coefficient of the polariton branch.
$P_{int}(k)$ is the incident power on the top mirror and $\eta$ is the
coupling coefficient. $t_{\text{top mirror}}$ is the transmission of top
mirror, and $\tau_{trip}$ is the photon trip time in the cavity. This gives
$\eta=4.1\times 10^{12}$ s-1. Equation S2 can be written as

$\begin{gathered}\tilde{\psi}_{LP}(k,\omega_{P})=\frac{\hbar
F_{P}(k)}{\hbar\omega_{P}-\epsilon_{LP}(k)+\frac{i\hbar\gamma_{LP}}{2}}\end{gathered}$
(S5)

Polariton density at the wavevector k and frequency $\omega_{P}$ is

$\begin{gathered}\left|\psi_{LP}(k,\omega_{P},t)\right|^{2}=\left|\tilde{\psi}_{LP}(k,\omega_{P})\right|^{2}=\frac{|C|^{2}\frac{\eta\cdot
P_{\text{int
}}(k)}{\hbar\omega_{P}}}{\left(\omega_{P}-\omega_{LP}(k)\right)^{2}+\left(\frac{\gamma_{LP}}{2}\right)^{2}}\end{gathered}$
(S6)

Since the polaritons have a finite width in energy, the total polariton
density at the wavevector k can be found out by integrating over the all
frequency range

$\left|\tilde{\psi}_{LP}(k)\right|^{2}=\int\frac{|C|^{2}\frac{\eta\cdot\xi(k,\omega_{P})}{\hbar\omega_{P}}}{\left(\omega_{P}-\omega_{LP}(k)\right)^{2}+\left(\frac{\gamma_{LP}}{2}\right)^{2}}d\omega_{P}$
(S7)

Here, $\xi(k,\omega_{P})$ denotes the incident power density
(power/frequency/wavevector).

The integral above can be numerically calculated but we can make an
approximation to get a close form solution for the density. We note that the
Lorentzian part of the integrand for $\gamma_{LP}\rightarrow 0$ is a Dirac
delta function

$\delta(\omega_{P}-\omega_{LP})=\lim_{\gamma_{LP}\rightarrow
0}\frac{1}{\pi}\frac{\frac{\gamma_{LP}}{2}}{(\omega_{P}-\omega_{LP})^{2}+(\frac{\gamma_{LP}}{2})^{2}}$
(S8)

Using the above formula and taking the power density $\xi(k,\omega_{P})$
constant over the narrow range of excitation we get

$\left|\tilde{\psi}_{LP}(k)\right|^{2}\approx\frac{2\pi|C(k)|^{2}\eta\cdot\xi}{\hbar\omega_{LP}(k)\gamma_{LP}(k)}$
(S9)

Here, $\gamma_{LP}(k)$ is determined from the experimental data.

Total real space polariton density can be found out by summing
$\left|\tilde{\psi}_{LP}(k)\right|^{2}$ in the k-space for the experimental
wavevector range

$\begin{split}\left|\tilde{\psi}_{LP}\right|^{2}&=\sum_{k=k_{min}}^{k_{max}}\left|\tilde{\psi}_{LP}(k)\right|^{2}\\\
&=\frac{2\pi\eta\cdot\xi}{\hbar}\sum_{k=k_{min}}^{k_{max}}\frac{|C(k)|^{2}}{\omega_{LP}(k)\gamma_{LP}(k)}\end{split}$
(S10)

If $P_{0}$ is the measured real space peak power density (power/energy) of the
pulsed supercontinuum laser excitation then $\xi=P_{0}/N$, where N is the
number of k point in between $k_{min}$ and $k_{max}$.

The above equation can also be written as

$\begin{split}\left|\tilde{\psi}_{LP}\right|^{2}&=\sum_{k=k_{min}}^{k_{max}}\frac{2\pi}{\omega_{\mathrm{LP}}(k)\,\tau_{\mathrm{trip}}}\times\frac{|C(k)|^{2}|t_{\mathrm{top\,mirror}}|^{2}\,R}{\gamma_{\mathrm{LP}}(k)}\end{split}$
(S11)

Here, $\frac{2\pi}{\omega_{\mathrm{LP}}\tau_{trip}}\sim 1$ and $R$ is the
incident number of photons/sec.

##  V Raw data for the nonlinearity of the IE polariton

Figure S4:  Raw data for the nonlinearity of the IE polariton at 7 K. (a)
K-space differential reflection data showing the dispersion of the IE lower
branch (pol-2) and IE upper branch (pol-3) at the lowest density (1.9 µm-2).
The excitation bandwidth of the supercontinuum laser for the IE is set 40 meV
by using a low-pass and a high-pass filter. (b) Line cut of the k-space data
at zero detuning $k_{\parallel}$ showing the differential reflection of the IE
lower branch (pol-2) and IE upper branch (pol-3) at different polariton
density. Note that since the IE exciton blueshifts with increasing density,
the zero detuning $k_{\parallel}$ also increases with increasing density. All
the power dependent nonlinear measurements (for both samples in the cavity and
outside the cavity) were carried out using a pulsed supercontinuum laser (20
ps pulsewidth) with proper bandpass filter in the input to excite only one
polariton species.

Fig. S4a shows the differential reflection of the IE polaritons excited (with
a 20 ps pulsewidth supercontinuum laser) in a narrow band used to study the
nonlinearity of the IE polariton. Fig. S4b shows the differential reflection
at the zero detuning $k_{\parallel}$ for IE polariton. The upper polariton
branch red shifts and lower polariton branch blue shifts with increasing
density reducing the Rabi splitting. We also notice that the two branches move
asymmetrically – the lower branch moves more than the upper branch at zero
detuning $k_{\parallel}$. This is because both exciton-exciton interaction and
saturation give rise to blue shift for the lower branch but they produce
opposite shift for the upper branch – as described in the main text.

##  VI Fits to determine the blueshift of the bare excitons

Figure S5:  Power dependent reflectance measurement to measure the blue shift.
Three Lorentzian centered at A exciton, IE and B exciton respectively are
fitted to the data at each power. The blue shift of the Lorentzians centered
at A exciton and IE are plotted in the main text as a function of estimated
exciton density. The method of exciton density estimation is described in
Methods.

##  VII Rabi splitting of interlayer exciton polariton as a function of
polariton density

Figure S6:  Rabi splitting of IE polariton at the zero detuning
$k_{\parallel}$ as a function of polariton density at 7 K. The blue line is
the fit to the usual formula
$\Omega=\frac{\Omega_{0}}{\sqrt{1+\frac{n}{n_{\text{sat}}}}}$ for free 2D
excitons.

We plot the Rabi splitting of the IE polariton at the zero detuning
$k_{\parallel}$ as a function of density, see Fig. S6. We mote that the usual
formula for free 2D excitons
$\Omega=\frac{\Omega_{0}}{\sqrt{1+\frac{n}{n_{\text{sat}}}}}$ does not fit our
data well. The inadequacy of the free 2D exciton formula to describe the
density dependence of the Rabi splitting was also observed in a recent study 1
on the Moiré exciton-polariton in a hetrobilayer TMDC.

##  VIII Energy of the lower branch of A exciton polariton as a function of
polariton density

Figure S7:  Energy of the lower branch of A exciton polariton as a function of
polariton density at 7 K from which the $g_{LP}^{A}$ is calculated and plotted
in Fig.3b of the main manuscript.

##  IX Separation of the saturation nonlinearity and exciton-exciton
interaction nonlinearity for IE polariton

Figure S8:  Separation of the saturation nonlinearity and exciton-exciton
interaction nonlinearity for IE polariton. Calculated $\Delta\mathrm{E_{exc-
exc}}$ and $\Delta\mathrm{E_{sat}}$ at the zero detuning $k_{\parallel}$ from
the experimental data as a function polariton density.

##  X Probing non-linearities in the reflectance spectrum of bilayer MoS2 on
two distinct time scales

Figure S9:  Probing non-linearities in the reflectance spectrum of bilayer
MoS2 after the excitons scatter from K-valley. (a) Transient reflectivity
spectra probed as a function of white-light probe energy and pump-probe delay.
(b) Transient reflectivity spectra probed 1.5 ps after the resonant excitation
of the A (full), B (dashed) and IE (dot-dashed) excitons. The absorbed photon
density is about 4x104 photons/µm${}^{\textrm{2}}$, 2x104
photons/µm${}^{\textrm{2}}$ and 4x104 photons/µm${}^{\textrm{2}}$
respectively.

Fig. S9a shows transient reflectivity measured on the bilayer film. The sample
is first excited by an ultrashort pump pulse (200 fs) resonant with one of the
exciton transitions. The relative change in reflectivity induced by the pump
is then probed as a function of time delay using a broadband pulse. The
transient reflectivity spectra obtained when pumping the A, B and IE
resonances at a pump-probe delay of 1.5 ps are plotted in Fig. S9b. Regardless
of the exciton transition that is initially excited, all the excitonic
resonances redshift and saturate (which corresponds to the observed
differential lineshapes). Time traces of the transient reflectivity are shown
in Fig. S10.

In contrast, we notice that all the excitonic features blueshift when the
sample is probed with a single broadband pulse lasting 20 ps (main text Fig.
1d). The opposite shifts observed using these two probes provide insight in
the excitation regimes induced in the bilayers. When the probe is delayed from
the pump all excitonic features are shifted to lower energies (Fig. S9a). This
is consistent with the expected contribution of relaxed carriers in these
materials 2, 3, 4. Indeed, the pump-probe delay used in this experiment is
much greater than the K-valley exciton lifetime (45 fs for A exciton and and
53 fs for IE). In this case, the probe interacts with the sample when only
relaxed carriers remain.

When the reflectivity is probed while the material is still being pumped, as
is the case in the single pulse reflectivity experiment, the dominant non-
linear effect is instead a shift of all features to higher energies. Indeed,
when the pulse interacts with the sample, A, B and IE excitons are
continuously created while their constituents quickly relax out of the K
valley. This quasi-continuous pumping leads to the steady-state presence of A,
IE and B excitons amidst a growing background of relaxed carriers until the
pulse no longer interacts with the sample. Moreover, the blueshift arising
from excitonic repulsion dominates over the (smaller) redshift from the
background of relaxed carriers. Given the long pulses used in our single pulse
reflectivity, compared to the lifetime of the excitons ($\sim$ 50 fs), there
are about 20 ps/50 fs = 400 times more relaxed carriers than excitons. The
strength of the K-valley exciton-exciton interactions are therefore much
higher than that due to the interaction with the background carriers.

Figure S10:  Probing non-linearities in the reflectance spectrum of bilayer
MoS2 after the excitons scatter from K-valley. Transient reflectivity spectra
probed as a function of white-light probe energy and pump-probe delay. The
legend corresponds to the probe exciton energy where the time trace is taken.

## References

  * 1 Zhang, L. _et al._ Van der waals heterostructure polaritons with moiré-induced nonlinearity.  _Nature_ 591, 61–65 (2021).
  * 2 Schmidt, R. _et al._ Ultrafast coulomb-induced intervalley coupling in atomically thin ws2.  _Nano letters_ 16, 2945–2950 (2016).
  * 3 Steinhoff, A., Rosner, M., Jahnke, F., Wehling, T. O. & Gies, C.  Influence of excited carriers on the optical and electronic properties of mos2.  _Nano letters_ 14, 3743–3748 (2014).
  * 4 Sie, E. J. _et al._ Observation of exciton redshift–blueshift crossover in monolayer ws2.  _Nano letters_ 17, 4210–4216 (2017).