# Effectiveness of the Recent Advances
in Capsule Networks

Nidhin Harilal∗
Department of Computer Science
University of Colorado, Boulder
<EMAIL_ADDRESS>
Rohan Patil∗
Department of Computer Science
University of California, San Diego
<EMAIL_ADDRESS>

###### Abstract

Convolutional neural networks (CNNs) have revolutionized the field of deep
neural networks. However, recent research has shown that CNNs fail to
generalize under various conditions and hence the idea of capsules was
introduced in 2011, though the real surge of research started from 2017. In
this paper, we present an overview of the recent advances in capsule
architecture and routing mechanisms. In addition, we find that the relative
focus in recent literature is on modifying routing procedure or architecture
as a whole but the study of other finer components, specifically, squash
function is wanting. Thus, we also present some new insights regarding the
effect of squash functions in performance of the capsule networks. Finally, we
conclude by discussing and proposing possible opportunities in the field of
capsule networks.

* These authors contributed equally to this work

##  1 Introduction

Over the last few years, neural networks have made remarkable progress in
various tasks ranging from vision tasks like image recognition and object
segmentation to machine translation tasks. The availability of huge amounts of
data has made it possible for Neural Networks to excel in different areas of
Computer Vision. The considerable success lies in the fact that CNNs can
automatically extract high-level features from images, which are much more
powerful than human-designed features.

Capsule Networks (CapsNet) were introduced by Hinton et al. [12, 34], which
addressed the significant limitations of CNNs and showed superior performance
on the MNIST [22] dataset. Presently, capsule networks are regarded as one of
the most promising breakthroughs in deep learning. The primary reason behind
this is that capsules’ idea provides a much more promising way of dealing with
different variations in images, including position, scale, orientation, and
lighting, than the currently employed methods in the neural networks
community.

Despite capsule networks showing an increased positive impact on various
vision tasks [17, 2, 3], the lack of architectural knowledge has limited
researchers to exploit the full potential of this new field. Therefore, this
paper aims to provide insights behind the working of capsule networks and
critically review the latest advances in this field. First, we briefly discuss
the reason behind the introduction of capsule networks, followed by a
description of its components. Then, we provide an analysis describing aspects
and limitations of various components of different capsule networks, followed
by an analysis of augmenting squash functions. Lastly, we conclude our
analysis by describing opportunities for further research in this field.

##  2 Convolutional Neural Networks (CNNs)

Convolutional neural networks (CNN) are feed-forward neural networks that can
extract features from data in a hierarchical structure. The architecture of
CNN is inspired by visual perception [16]. CNN architectures date back decades
[21], consisting of convolutional layers that can extract high-level features
from images. CNNs detect these features in images and learn how to recognize
objects with this information. Layers near the start detect simpler features
like edges, and as the layers get deeper, they detect more complex features
like eyes, noses, or an entire face in case of face recognition. It then uses
all of these features, which it has learned to make a final prediction. Deep
CNNs have provided a significant contribution in computer vision tasks such as
image classification [18, 9]. They have also been successfully applied to
other computer vision fields, such as object detection [32, 4, 25], face
recognition [29], etc.

#### Limitations:

CNNs perform exceptionally great when they are inferenced over images that
resemble the dataset [7]. Despite being successful, CNN performs poorly when
it receives the same image with a different viewpoint [1, 26]. Convolving
kernel across an image ensures invariance, but it doesn’t certify
equivariance. Max-pooling [30] was introduced to further aid in creating the
positional invariance. On close scrutiny, the pooling operation stacked with a
convolutional layer will only detect the features but not preserve any spatial
relationships between the detected features. The difference could be explained
with a CNN detecting a human face. An average human face will have a pair of
eyes, nose, and mouth; however, an image in which these parts are present
doesn’t qualify it to be an image of a human face. Therefore, a system just
detecting certain features in the image without having information about their
spatial arrangement could fail in many scenarios. In short, The pooling
operation, along with the convolution, was supposed to introduce positional,
orientational, and proportional invariances but rather became a cause for them
[26]. Including more data or using methods like data, augmentation is used to
tackle this problem, which ensured that the model was trained on as many
viewpoints/ orientations as possible. However, this is a very crude way of
handling this.

##  3 Capsule Networks (CapsNet)

Hinton et al. [12] proposed the first capsule networks, which, unlike
conventional CNNs, were designed to encode and preserve the underlying spatial
information between its learned features. The basic idea behind the
introduction of capsules by Hinton et al. [12] was to create a neural network
capable of performing inverse graphics [11]. Computer graphics deals with
generating a visual image from some internal hierarchical representation of
geometric data. This internal representation consists of matrices that
represent the relative positions and orientation of the geometrical objects.
This representation is then converted to an image that is finally rendered on
the screen. Hinton and his colleagues in their paper [12] argued that humans
deconstruct a hierarchical representation of the visual information received
through eyes and use this representation for recognition. From a pure machine
learning perspective, this means that the network should be able to
deconstruct a scene into co-related parts, which can be hierarchically
represented [19]. To achieve this, Hinton et al. [12] proposed to augment the
conventional idea of neural network architecture to reflect the idea of
several entities. Each entity encapsulates a certain number of neurons and
learns to recognize an implicitly defined visual entity over a limited domain
of viewing conditions and deformations.

###  3.1 Evolution of Capsules

In the first proposed capsule networks [12], the output of each capsule
consisted of the probability that a specific feature exists along with a set
of instantiating parameters111Instantiating parameters [12] may include and
encode the underlying pose, lighting, and deformation of a particular feature
relative to the other features that are detected by the capsules in the
image.. The goal of Hinton et al. [12] was not recognizing or classifying the
objects in an image but rather to force the outputs of a capsule to recognize
the pose in an image. A significant limitation of this first implementation of
capsule networks was that it required an additional external set of inputs
that specified how the image had been transformed for each of the entities to
work. Although these transformations could be learned in principle, the more
significant challenge was to devise a way to instruct each of the capsules to
discover the underlying hierarchical relationship between the transformations
in a complete end-to-end train setting without explicitly using additional
inputs.

To overcome this challenge, Sabour et al. [34] proposed a new definition for
capsules. Each capsule consisted of a series of convolutional layers followed
by a fully-connected layer. Unlike the output from conventional artificial
neurons, which are scalar, the output from their proposed capsules is
represented as an activity vector. The length of this vector encoded the
probability of detecting a feature, and the state of the detected feature,
also referred to as instantiating parameters, is encoded as the direction of
these vectors. To represent the capsule network as vectors whose length
represents probabilities, they used a non-linear "squashing" function
described in Equation 1.

$v_{j}=\frac{{||s_{j}||}_{2}^{2}}{1+{||s_{j}||}_{2}^{2}}\frac{s_{j}}{{||s_{j}||}_{2}}$
(1)

where $v_{j}$ is the vector output of $j^{th}$ capsule and $s_{j}$ is its
total input. The above-defined non-linearity ensures that the short vector
gets shrunk to zero length and long vectors to slightly below 1. This is the
basic definition of capsules, and most further improvements in capsule
networks utilized this definition of capsules.

###  3.2 Routing mechanism

The routing mechanism is one of the key factors in making the capsule networks
work [34]. Routing controls the transfer of information between the capsules.
Higher-level capsules get the input from the lower level capsules. The first
and the last capsule layers are referred to as the Primary capsule and class
capsule layer, respectively. Feature extraction from the image is done by the
convolutional layers present inside each capsule, and the output is fed
according to the computed squash vectors. Figure 1 first routing method for
the capsule network was proposed by Sabour et al. [34].

Figure 1: First Routing based Capsule structure (Dynamic Routing) by Sabour et
al. [12]

In this, each capsule $i$ (where $1\leq i\leq N$) in layer $l$ has an activity
vector $\mathbf{u_{i}}\in\mathbb{R}$ to encode the spatial information of
features in the form of instantiating parameters using the activation
function. This output vector $\mathbf{u_{i}}$ of the $i^{th}$ lower-level
capsule is then fed into all capsules in the neighbouring high-level
$(l+1)^{th}$ layer. The total input to a capsule $s_{j}$ is a weighted sum of
"prediction vectors" $\hat{u}_{j|i}$ it received from the lower-level
capsules. For all layers of capsules except the primary capsules, The $j^{th}$
capsule at layer $l+1$ will receive $u_{i}$ and find its product with a
corresponding weight matrix $\mathbb{W}_{ij}$ which it is learns during
training.

$\mathbf{s}_{j}=\sum_{i}c_{ij}\hat{u}_{j|i},\;\;\hat{u}_{j|i}=\mathbb{W}_{ij}\mathbf{u}_{i}$
(2)

Coupling coefficients ($c_{ij}$) ensures that the output of capsule at level
$l$ is linked to that in layer $l+1$. These coefficients are determined by a
"routing softmax" whose initial logits $b_{ij}$ are the log prior
probabilities that capsule $i$ should be coupled to capsule $j$.

$c_{ij}=\frac{\exp{(b_{ij})}}{\sum_{k}\exp{(b_{ik})}}$ (3)

These coefficients are then updated iteratively based on the agreement between
the current output $v_{j}$ of each capsule, $j$, in the layer above and the
prediction $\hat{u}_{j|i}$ made by capsule $i$. In short, lower level capsules
are responsible for sending its input to certain higher-level capsules based
on the computed coefficients such that it “agrees” in establishing a relation
between the input and the task at hand. This is the essence of the dynamic
routing method for capsule networks and this mechanism of "agreement" or
"voting" between the capsules is referred as routing-by-agreement [34].

##  4 Advances in Routing Mechanism

The routing method is the distinguishing factor of capsule networks as it
provides the mechanism to learn the hierarchial representations, essentially
carving out a parse tree. The limitations of dynamic routing for generalizing
over complex data have led researchers to introduce other routing methods.

#### EM Routing:

It is an Expectation-Maximisation routing procedure [13]. The two main ideas
of this method are: first, the relationship between entity and viewpoint can
be encoded as a $4\times 4$ matrix, and second, the coefficients will be
iteratively updated so that a cluster of similar votes activate the next layer
capsules. According to Hinton et al. [13], this method is similar to
attention; the clustering is focusing on specific features, but here the
direction of attention is reversed in the sense that the competition is
between higher layer capsules for votes from lower-level capsules. Another
aspect is the robustness of adversarial attacks. Capsule networks with EM
Routing show better robustness than CNNs, but interestingly perform similar to
a CNN when adversarial examples are generated in a black-box manner using a
CNN.

#### KL Routing:

Wang et al. [36] present a routing procedure as an optimization process with a
KL regularization between current and previous state coupling coefficients.
For convenience, we call it KL Routing throughout the paper. In their setting,
they eliminate the squash function in the iteration part and calculate the
squash vector (calculating the squashed vector is different from the original)
at the end of iterations. The critical part is that they derive their provided
robust mathematical modeling of routing as an agglomerative fuzzy K-Means
algorithm [24].

#### Self Routing:

It is a Mixture-of-Experts [28] approach towards routing [8]. The process is
non-iterative, thus reducing the computation time. Moreover, it shows better
accuracy and robustness towards adversarial attacks as compared to EM routing.
The authors also claim that the $4\times 4$ matrix multiplication constraint
used in EM routing is too restrictive for learning a good representation of
complex data with multi-level routing layers. There are two learnable weight
matrices $W^{route}$ and $W^{pose}$, which are used for calculating the
interaction coefficients and prediction vectors in Self Routing.

Table 1: Comparing Routing Procedures.

Routing | Inspiration | Squash function | Iterative
---|---|---|---
EM Routing | Gaussian Mixture Model | None | Yes
KL Routing | Agglomerative Fuzzy K-Means algorithm | $v_{j}=\frac{||\sum_{i}c_{ij}o_{j|i}||}{1+max_{k}||\sum_{i}c_{ik}o_{k|i}||}s_{j}$ | Yes
Self Routing | Mixture-of-Experts | $v_{j}=\frac{{||s_{j}||}_{2}^{2}}{1+{||s_{j}||}_{2}^{2}}\frac{s_{j}}{{||s_{j}||}_{2}}$ | No

Here $o_{j|i}=\frac{1}{{||W_{ij}||}_{\mathcal{F}}}W_{ij}\mu_{i}$

##  5 Capsule Network Variants and Modifications

This section will highlight a few modifications and developments on the
original capsule network, which are not limited to routing, summarised in
Table 2

#### CapProNet

In CapProNet, each capsule is associated a subspace of the input vector [40].
Formally, each capsule learns a weight matrix $W$, the columns of which form
the subspace’s basis. The input vector is projected on these subspaces before
the routing procedure. The authors claim that the advantage of capsule
subspaces is that it is possible to find the best direction along which a
specific feature’s information about the position, orientation, scales, and
intrinsic characteristics can be represented better.

#### RS-CapsNet

The major motivation of RS-CapsNet [38] is from ResNet [10] and Squeeze-and-
Excitation Network [14]. Their work on routing is majorly to mitigate the
issues of stacked capsule layers, which do not perform well [37]. In RS-
CapsNet, a linear combination of the capsule outputs is used for the next
layer’s input. According to the authors, since the primary capsules are
reshaped from convolutional feature maps, redundant background information is
present. Taking a linear combination of capsule outputs will reduce the
background noise as well as the network will have increased non-linearity.

#### Sparse Unsupervised Capsules (SUPCAPS)

The work of Rawlinson et al. [31] shows how capsules can be used in an
unsupervised setting. Furthermore, their work shows the importance of using
sparsity in unsupervised learning. In their approach to creating a generator,
they initially removed the mask and margin loss function, but the result was
that the network started to act as an autoencoder. Furthermore, individual
learning capsules were unable to learn about specific features. To mitigate
this issue, they added dropout as it was observed that the coefficients of
capsules were similar, implying that the capsules were unable to carve a parse
tree. Adding dropout restored the ability to route in capsules. The discovery
that latent capsules lose their ability to specialize in their identity
matches the observation that stacked capsule layers have decreased performance
[37].

Table 2: Capsule Network Varinats and Modifications Variants |  Modification |  Contributions
---|---|---
CapProNet |  Each capsule learns a subspace and uses projection of input vector. |  Better representation of information. Increased speed.
RS-CapsNet |  Splice convolutional feature map to create ResNet type architecture. |  Allows stacking of capsule and increased non-linearity.
SUPCAPS |  Modified to use capsules as generator in unsupervised setting and added dropout. |  Study on loss of capsule learning capabilities and importance of dropout for capsules.

##  6 Discussion

The focus on routing has led people to search for various types of routing
procedures. The results on CIFAR10 [39] in [8] shows that EM routing performs
worse than that of the dynamic routing. In contrast, EM routing on SmallNORB
[23] gives similar results as that of Dynamic Routing, but the self-routing
procedure beats them all.

For KL Routing, experiments were done on a simple unsupervised perceptual
grouping task [6], where KL Routing performs better than dynamic routing.
There is no comparison between EM Routing and KL Routing, but given the type
of task, it appears that EM routing will also perform well. Given the
performance of EM routing on SmallNORB and Self Routing beating both Dynamic
and EM Routing, Self Routing will likely perform well on such tasks.
Furthermore, capsule networks give better performance when used in an ensemble
[37], supporting the Mixture-of-Experts approach towards routing.

It has been shown that residual networks behave like ensembles of networks
with different depths [35] and the authors for Self Routing mentions
integration with residual networks as future work. In contrast, the routing
procedure of RS-CapsNet slices the convolutional feature maps to construct
capsules, and the routing is done on these slices. This procedure can also be
construed as doing Self Routing on complete input given to multiple capsules.
Moreover, the slicing of input and then taking a linear combination of each is
the same as projecting the complete input vector onto different vector
subspaces. This approach is similar to CapProNet. Zhang et al. [40] apply the
subspace capsules only as of the last layer with other networks like ResNet
and DenseNet; their future work included the integration of the same in all
layers. It is interesting to note that RS-CapsNet performs better in
ensembles, so the same may hold for capsules networks with Self Routing.

In EM routing, we previously mentioned that black box attacks gave performance
similar to CNNs, which likely indicates that the methods used in white box
methods [5, 20] for generating adversarial examples are not effective against
capsules or at least are not smart enough. The work is done by Michels et al.
[27] shows that it is possible to fool capsule nets with sophisticated
attacks, though they did modify the algorithms, and for simpler attacks where
CNNs fail, capsule network does show robustness to attacks. Furthermore, their
work only consider the capsule network proposed by Sabour et al. [34] but not
the other variants, so it is not possible to generalize. Similarly, given the
similarities between Self Routing and RS-CapsNet, it is possible to also be
robust towards adversarial attacks, though this needs to be confirmed.

Rawlinson et al. [31] expects that sparse unsupervised capsules (SUPCAPS) will
also show robustness similar to capsules with EM Routing because of the
unsupervised training regime. Though this is just a speculation, it is based
on the improved generalization when dropouts were added. It appears that the
capsules need to be given a chance to explore other possibilities for learning
the underlying hierarchical relationship given the same inputs, and it seems
that dropouts are helping in this process.

Autoencoder based capsule structure by Hinton et al. [12] required to give
information about the augmentation in image. In the work of Sabour et al.
[34], this problem appears to have been solved by using the margin loss and
masking, but not specifically due to structuring of capsule. This is
observable from the fact that when margin loss and mask is removed, the
network starts to act as an autoencoder [31]. The failure of stacked capsules
[37] can be used to interpret that capsules in the hidden layer are not
learning the latent information of features. Possibly, this is because
capsules require margin loss to force them to explore different paths and thus
capsules require the special supervision of a loss directly.

##  7 Squash Function

The squash function plays the role of adding non-linearity in the routing
process. An attempt to modify the squash function was done in [37, 36, 15],
but there has not been much work. In the previously mentioned routings and
variants, equation 1 is used or squashing itself is removed as in the case of
EM Routing. We did not find a proper study of squash functions and hence did
some experiments, though they were limited to not changing the direction of
the vector. We propose the following variants of squash function:

$v_{j}=\frac{{||s_{j}||}_{m}^{m}}{1+{||s_{j}||}_{m}^{m}}\frac{s_{j}}{{||s_{j}||}_{2}}\,\,\,\,where\,\,m\in\mathbf{I}^{+}$
(4)
$v_{j}=\frac{{||s_{j}||}_{\infty}}{1+{||s_{j}||}_{\infty}}\frac{s_{j}}{{||s_{j}||}_{2}}$
(5)

For notation, we define the squash function represented by equation 4 as
$S_{m}$ and equation 5 as $S_{\infty}$. Our focus is on the relative
performance of $S_{1}$, $S_{2}$, $S_{3}$, $S_{4}$, $S_{5}$, $S_{10}$ and
$S_{\infty}$, for which we train the capsule network [34] on MNIST [22] and
CIFAR10 with five-fold validation.

Figure 2: Performance comparison of different norm-based squash functions. Table 3:  Results on MNIST and CIFAR10  | MNIST | CIFAR10
---|---|---
Function | Accuracy (%) | Loss | Accuracy (%) | Loss
$S_{1}$ | $99.408\pm 0.095$ | $29.179\pm 0.96$ | 66.108 $\pm$ 2.09 | 84.697 $\pm$ 1.48
$S_{2}$ | 99.512 $\pm$ 0.028 | 27.495 $\pm$ 0.49 | 74.525 $\pm$ 0.23 | 74.038 $\pm$ 1.43
$S_{3}$ | 99.417 $\pm$0.021 | 28.636 $\pm$ 0.58 | 75.068 $\pm$ 0.87 | 71.873 $\pm$ 1.59
$S_{4}$ | 99.582 $\pm$ 0.037 | 27.261$\pm$0.31 | 73.727 $\pm$ 1.19 | 72.593 $\pm$ 1.62
$S_{5}$ | 99.542 $\pm$ 0.037 | 27.971 $\pm$ 0.38 | 73.506 $\pm$ 0.32 | 77.401 $\pm$ 1.41
$S_{10}$ | 13.521 $\pm$ 0.021 | 98.473 $\pm$ 0.19 | 11.269 $\pm$ 0.27 | 121.713 $\pm$ 0.85
$S_{\infty}$ | 10.897 $\pm$ 0.013 | 114.513 $\pm$ 0.15 | 10.681 $\pm$ 0.21 | 132.612 $\pm$ 0.79

We observe that for $S_{1}$, the train accuracy remain close to that of test
but the overall performance is worse than $S_{2}$. In contrast, for $m>2$, the
train accuracy starts to shoot while the test accuracy is close to $S_{2}$.
The likely reason for this is that with increasing $m$, votes with little
increase in magnitude can give very high agreement. However, when the
magnitude drops below unity, the vote is highly penalized. As
$m\rightarrow\infty$, the final magnitude as the function of original
magnitude will have a jump from $0$ to $1$ at value unity.

For $S_{\infty}$ we find that the capsules miserably fail to generalize. This
is likely because $S_{\infty}$ fails to penalize disagreements as well as to
give high weight to agreements. The same logic can be extended towards the
family $S_{m}$. But this logic is limited in the sense that if we highly
penalize disagreements, then once a capsule has decided a route in training,
it will fail to discover others.

We observe that for $S_{2}$ and $S_{5}$, the training accuracy shoots up in
the initial epochs itself and the similar can be seen in test accuracy also,
but this is not the case for $S_{3}$ and $S_{4}$. This shows that the value of
$m$ is deciding how much previously learned knowledge the network will
preserve, but no clear relation is visible. For $S_{2}$, we expect that the
network will moderately penalize disagreements while the case of $S_{5}$ is
the other extreme. The behaviour of $S_{3}$ and $S_{4}$ is somewhat in the
middle and the slow learning is likely because the network does not easily
form or forget an opinion in the initial stages.

##  8 Conclusions

It is evident that capsules require extreme supervision for specializing them,
thus limiting the possibility of achieving better performance from deep
capsule networks. Dropouts, residuals and subspaces provide a way to solve
this problem. We observe that there needs to be proper comparison between
various routing methods on multiple types of datasets as certain routing
methods work good only for specific tasks. Also, there is little work on the
squash function which primarily brings non-linearity in capsules. Our work on
squash function shows how greatly it can affect the overall learning process
of capsule networks, thus, an analysis on finer components of these networks
are worth an investigation.

##  9 Future Work

We believe that the on key problem in capsule networks regarding the
specializing comes up because the training of the complete network is done in
an end-to-end manner. In contrast, it is possible that capsules are added
slowly to a capsule layer. Further, we think that incorporating techniques
similar to Rawlinson et. al. [31], it is possible to force the network to
explore other routes for learning the underlying hierarchical relationships
with the newly added capsules. Moreover, using subspaces can allow forcing
capsules to learn specific features, though it might be necessary to add
constraints on learning of subspaces to ensure that all capsules learn a
different one. Also, incorporating routing methods like Self Routing can
mitigate the issue of previous training that did not include the vote of the
newly added capsule. This can be corresponded with the concept of incremental
learning that is done in CNNs [33]. Incorporating techniques like incremental
learning in capsule networks is just like putting clay slowly to add new
pathways to the previously learned representations between capsules.

## References

  * [1] Amr Bakry, Mohamed Elhoseiny, Tarek El-Gaaly, and Ahmed M. Elgammal.  Digging deep into the layers of cnns: In search of how cnns achieve view invariance.  In Yoshua Bengio and Yann LeCun, editors, 4th International Conference on Learning Representations, ICLR 2016, San Juan, Puerto Rico, May 2-4, 2016, Conference Track Proceedings, 2016.
  * [2] Fei Deng, Shengliang Pu, Xuehong Chen, Yusheng Shi, Ting Yuan, and Shengyan Pu.  Hyperspectral image classification with capsule network using limited training samples.  Sensors, 18(9):3153, 2018.
  * [3] Kevin Duarte, Yogesh Rawat, and Mubarak Shah.  Videocapsulenet: A simplified network for action detection.  In Advances in Neural Information Processing Systems, pages 7610–7619, 2018.
  * [4] Ross Girshick, Jeff Donahue, Trevor Darrell, and Jitendra Malik.  Rich feature hierarchies for accurate object detection and semantic segmentation.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 580–587, 2014.
  * [5] Ian J Goodfellow, Jonathon Shlens, and Christian Szegedy.  Explaining and harnessing adversarial examples.  arXiv preprint arXiv:1412.6572, 2014.
  * [6] Klaus Greff, Sjoerd Van Steenkiste, and Jürgen Schmidhuber.  Neural expectation maximization.  In Advances in Neural Information Processing Systems, pages 6691–6701, 2017.
  * [7] Jiuxiang Gu, Zhenhua Wang, Jason Kuen, Lianyang Ma, Amir Shahroudy, Bing Shuai, Ting Liu, Xingxing Wang, Gang Wang, Jianfei Cai, and Tsuhan Chen.  Recent advances in convolutional neural networks.  Pattern Recognition, 77:354 – 377, 2018.
  * [8] Taeyoung Hahn, Myeongjang Pyeon, and Gunhee Kim.  Self-routing capsule networks.  In Advances in Neural Information Processing Systems, pages 7658–7667, 2019.
  * [9] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  In Proceedings of the IEEE international conference on computer vision, pages 1026–1034, 2015.
  * [10] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [11] Geoffrey Hinton, Alex Krizhevsky, Navdeep Jaitly, Tijmen Tieleman, and Yichuan Tang.  Does the brain do inverse graphics.  In Brain and Cognitive Sciences Fall Colloquium, volume 2, 2012\.
  * [12] Geoffrey E Hinton, Alex Krizhevsky, and Sida D Wang.  Transforming auto-encoders.  In International conference on artificial neural networks, pages 44–51. Springer, 2011.
  * [13] Geoffrey E Hinton, Sara Sabour, and Nicholas Frosst.  Matrix capsules with em routing.  In International Conference on Learning Representations, 2018.
  * [14] Jie Hu, Li Shen, and Gang Sun.  Squeeze-and-excitation networks.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 7132–7141, 2018.
  * [15] Wenkai Huang and Fobao Zhou.  Da-capsnet: dual attention mechanism capsule network.  Scientific Reports, 10(1):1–13, 2020.
  * [16] D. Hubel and T. Wiesel.  Receptive fields, binocular interaction and functional architecture in the cat’s visual cortex.  The Journal of Physiology, 160, 1962.
  * [17] Ayush Jaiswal, Wael AbdAlmageed, Yue Wu, and Premkumar Natarajan.  Capsulegan: Generative adversarial capsule network.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 0–0, 2018.
  * [18] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  In Advances in neural information processing systems, pages 1097–1105, 2012.
  * [19] Tejas D. Kulkarni, Will Whitney, Pushmeet Kohli, and Joshua B. Tenenbaum.  Deep convolutional inverse graphics network.  CoRR, abs/1503.03167, 2015.
  * [20] Alexey Kurakin, Ian Goodfellow, and Samy Bengio.  Adversarial examples in the physical world.  arXiv preprint arXiv:1607.02533, 2016.
  * [21] Y. LeCun, B. Boser, J. S. Denker, D. Henderson, R. E. Howard, W. Hubbard, and L. D. Jackel.  Backpropagation applied to handwritten zip code recognition.  Neural Comput., 1(4):541–551, December 1989.
  * [22] Yann LeCun, Corinna Cortes, and Christopher JC Burges.  The mnist database of handwritten digits, 1998.  URL http://yann. lecun. com/exdb/mnist, 10(34):14, 1998.
  * [23] Yann LeCun, Fu Jie Huang, and Leon Bottou.  Learning methods for generic object recognition with invariance to pose and lighting.  In Proceedings of the 2004 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2004. CVPR 2004., volume 2, pages II–104. IEEE, 2004.
  * [24] Mark Junjie Li, Michael K Ng, Yiu-ming Cheung, and Joshua Zhexue Huang.  Agglomerative fuzzy k-means clustering algorithm with selection of number of clusters.  IEEE transactions on knowledge and data engineering, 20(11):1519–1534, 2008.
  * [25] Tsung-Yi Lin, Piotr Dollár, Ross Girshick, Kaiming He, Bharath Hariharan, and Serge Belongie.  Feature pyramid networks for object detection.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2117–2125, 2017.
  * [26] Rosanne Liu, Joel Lehman, Piero Molino, Felipe Petroski Such, Eric Frank, Alex Sergeev, and Jason Yosinski.  An intriguing failing of convolutional neural networks and the coordconv solution, 2018.
  * [27] Felix Michels, Tobias Uelwer, Eric Upschulte, and Stefan Harmeling.  On the vulnerability of capsule networks to adversarial attacks.  arXiv preprint arXiv:1906.03612, 2019.
  * [28] Steven J Nowlan and Geoffrey E Hinton.  Evaluation of adaptive mixtures of competing experts.  In Advances in neural information processing systems, pages 774–780, 1991.
  * [29] Omkar M Parkhi, Andrea Vedaldi, and Andrew Zisserman.  Deep face recognition.  British Machine Vision Association, 2015.
  * [30] Marc’Aurelio Ranzato, Fu Jie Huang, Y-Lan Boureau, and Yann LeCun.  Unsupervised learning of invariant feature hierarchies with applications to object recognition.  In 2007 IEEE conference on computer vision and pattern recognition, pages 1–8. IEEE, 2007.
  * [31] David Rawlinson, Abdelrahman Ahmed, and Gideon Kowadlo.  Sparse unsupervised capsules generalize better.  arXiv preprint arXiv:1804.06094, 2018.
  * [32] Joseph Redmon, Santosh Divvala, Ross Girshick, and Ali Farhadi.  You only look once: Unified, real-time object detection.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 779–788, 2016.
  * [33] Deboleena Roy, Priyadarshini Panda, and Kaushik Roy.  Tree-cnn: a hierarchical deep convolutional neural network for incremental learning.  Neural Networks, 121:148–160, 2020.
  * [34] Sara Sabour, Nicholas Frosst, and Geoffrey E Hinton.  Dynamic routing between capsules.  In Advances in neural information processing systems, pages 3856–3866, 2017.
  * [35] Andreas Veit, Michael J Wilber, and Serge Belongie.  Residual networks behave like ensembles of relatively shallow networks.  In Advances in neural information processing systems, pages 550–558, 2016.
  * [36] Dilin Wang and Qiang Liu.  An optimization view on dynamic routing between capsules, 2018.
  * [37] Edgar Xi, Selina Bing, and Yang Jin.  Capsule network performance on complex data.  arXiv preprint arXiv:1712.03480, 2017.
  * [38] Shuai Yang, Feifei Lee, Ran Miao, Jiawei Cai, Lu Chen, Wei Yao, Koji Kotani, and Qiu Chen.  Rs-capsnet: An advanced capsule network.  IEEE Access, 8:85007–85018, 2020.
  * [39] Matthew D Zeiler and Rob Fergus.  Stochastic pooling for regularization of deep convolutional neural networks.  arXiv preprint arXiv:1301.3557, 2013.
  * [40] Liheng Zhang, Marzieh Edraki, and Guo-Jun Qi.  Cappronet: Deep feature learning via orthogonal projections onto capsule subspaces.  In Advances in Neural Information Processing Systems, pages 5814–5823, 2018.


Limits on Leptonic TeV Emission from the Cygnus Cocoon with Swift-XRT

David Guevel

0000-0002-0870-2328]David Guevel
Department of Physics, Wisconsin IceCube Particle Astrophysics Center, University of Wisconsin, Madison, WI, USA

0000-0001-5186-5950]Andrew Beardmore
School of Physics & Astronomy, University of Leicester,
University Road,
Leicester, LE1 7RH, United Kingdom

0000-0001-5624-2613]Kim L. Page
School of Physics & Astronomy, University of Leicester,
University Road,
Leicester, LE1 7RH, United Kingdom

Amy Lien
University of Tampa,
401 W Kennedy Blvd,
Tampa, FL, USA

0000-0002-5387-8138]Ke Fang
Department of Physics, Wisconsin IceCube Particle Astrophysics Center, University of Wisconsin, Madison, WI, USA

0000-0001-7523-570X]Luigi Tibaldo
IRAP, Université de Toulouse, CNRS, CNES, UPS, 9 avenue Colonel Roche, 31028 Toulouse, Cedex 4, France

Sabrina Casanova
Institute of Nuclear Physics Polish Academy of Sciences, PL-31342 IFJ-PAN, Krakow, Poland

Petra Huentemeyer
Department of Physics, Michigan Technological University, Houghton, MI, USA

$\gamma$-ray observations of the Cygnus Cocoon, an extended source surrounding the Cygnus X star-forming region, suggest the presence of a cosmic ray accelerator reaching energies up to a few PeV.
The very-high-energy (VHE; 0.1-100 TeV) $\gamma$-ray emission may be explained by the interaction of cosmic-ray hadrons with matter inside the Cocoon, but an origin of inverse Compton radiation by relativistic electrons cannot be ruled out.
Inverse Compton $\gamma$-rays at VHE are accompanied by synchrotron radiation peaked in X-rays. Hence,
X-ray observations may probe the electron population and magnetic field of the source.
We observed eleven fields in or near the Cygnus Cocoon with the Neil Gehrels Swift Observatory's X-Ray Telescope (Swift-XRT) totaling 110 ksec.
We fit the fields to a Galactic and extra-galactic background model and performed a log-likelihood ratio test for an additional diffuse component.
We found no significant additional emission and established upper limits in each field.
By assuming that the X-ray intensity traces the TeV intensity and follows an $dN/dE\propto E^{-2.5}$ spectrum, we obtained a 90% upper limit of $F_X < 8.7\times 10^{-11}\rm~erg\,cm^{-2}\,s^{-1}$ or $< 5.2\times 10^{-11}\rm~erg\,cm^{-2}\,s^{-1}$ on the X-ray flux of the entire Cygnus Cocoon between 2 and 10 keV depending on the choice of hydrogen column density model.
This suggests that no more than one quarter of the $\gamma$-ray flux at 1 TeV is produced by inverse Compton scattering, when assuming an equipartition magnetic field of $\sim 20\,\mu$G.

§ INTRODUCTION

The cosmic ray spectral energy distribution features a spectral break, commonly called “the knee," at a few PeV.
The knee is believed to mark the transition from Galactic to extra-galactic cosmic ray origin.
The accelerators that produce the highest energy Galactic cosmic rays, known as PeVatrons, have not been conclusively identified.
The Cygnus Cocoon, from which $\gamma$-rays above 100 TeV have been detected <cit.>, is a prime PeVatron candidate.
The Cygnus Cocoon contains the Cyg OB2 stellar association and NGC 6910 stellar cluster, as well as $\gamma$-Cygni Supernova Remnant but extends beyond them.
The winds of massive stellar clusters inside the Cocoon may potentially accelerate particles to PeV energies  <cit.>.

GeV and TeV $\gamma$-rays have been seen by the Fermi Large Area Telescope (LAT;  ) up to 100 GeV, the High Altitude Water Cherenkov Observatory (HAWC; ) beyond 100 TeV <cit.>, and by the Large High Altitude Air Shower Observatory (LHAASO) up to 1.4 PeV <cit.>.
The GeV emission detected by the Fermi-LAT indicates an cosmic-ray energy density that is 50% larger than the ISM with a harder spectrum <cit.>, suggesting that cosmic ray particles are freshly accelerated at the Cocoon. The GeV emission may be explained by either leptonic or hadronic models.
The HAWC Collaboration found a change in the spectral shape at around 1 TeV and suggests that the morphology and spectral shape are consistent with proton interaction, though a leptonic contribution cannot be ruled out <cit.>.

The electron population that inverse Compton scatters will necessarily produce synchrotron X-rays in the presence of a magnetic field.
The Klein Nishina effect suppresses inverse Compton scattering for optical and UV photons so that the TeV emission mostly traces the diffuse dust radiation field rather than the stellar radiation field.
An X-ray counterpart to the Cocoon will trace the extended TeV emission if some of the TeV emission is in fact produced by inverse Compton scattering.

<cit.> used Suzaku to constrain the X-ray 2–10 keV intensity from the Cocoon to less than $2.35\times 10^{-8}\rm~erg\,cm^{-2}\,s^{-1}\,sr^{-1}$ based on two on-source observations within the Cocoon and two off-source observations outside the Cocoon.
The Suzaku 2–10 keV intensity decreases monotonically with Galactic latitude, but the brightest source is actually labelled as one of their background targets (BG1).
The 2–10 keV upper limits are calculated by taking the two on-source observations (Source1 and Source2) and subtracting the fainter of the two backgrounds (BG2).
The Suzaku two off-source regions are within the TeV emitting region found by HAWC which was not known at the time, so they may include some signals.
The upper limit on X-ray flux from the entire Cocoon is found by extrapolating the upper limit intensity to the size of the Cocoon ($4.38\times10^{-3} \,\rm~sr$) finding an upper limit flux density of $6.41\times10^{-11}~\rm erg\,cm^{-2}\,s^{-1}$ at 1 keV assuming a spectral index equal to 2.

Multi-wavelength map of the Cygnus Cocoon. Red is the Planck CO map <cit.>, green is a WISE 12 micron mosaic <cit.>, and blue is the HAWC significance map above 1 TeV using 1,343 days of HAWC data <cit.>. The blue cutoff   approximately corresponds to $3\sigma$. White circles indicate the 11 fields observed by Swift-XRT and analyzed in this work.

We observed eleven fields within and around the Cygnus Cocoon with the Neil Gehrels Swift Observatory's X-Ray Telescope (Swift-XRT) to search for evidence of leptonic emission.
We establish 2–10 keV upper limits on the intrinsic (unabsorbed) Cocoon intensity in 7 fields assuming an $E^{-2.5}$ spectrum as suggested by the leptonic model by <cit.>.

Additionally, we derive upper limits on the solid angle integrated Cocoon emission assuming that the X-ray emission traces the TeV emission.

We describe our observations and data processing in Section <ref>. Our spectral fitting is described in Section <ref>. We discuss the implications for the leptonic TeV emission in Section <ref>.

Summary of Observations.

Field     Target ID     l [Deg]     b [Deg]     Exposure [s]
1     95932     81.418     +0.921     10007

2$^*$     95933     80.341     +1.597     9817

3$^*$     95934     79.714     +1.510     8256

4     95935     78.295     -0.975     8309

5     95936     78.327     +0.586     9282

6     95937     79.710     +2.684     9610

7$^{\dagger}$     95938     81.683     +3.110     9025

8$^{\dagger}$     95939     77.457     +3.945     7554

9$^{*\dagger}$     95940     76.973     -0.402     9012

10$^{\dagger}$     95941     81.297     -0.622     7994

11     15110     81.246     +1.521     10723

Dagger ($\dagger$) indicates that the field is outside of the Cocoon and used as background.
Star ($*$) indicates a field that is contaminated by stray light. The treatment of stray light contamination is discussed in Section <ref>.

§ OBSERVATIONS AND DATA PROCESSING

§.§ Observations

Swift observed ten fields in or near the Cygnus Cocoon from May 27, 2021 to December 19, 2021, and one additional field as a target of opportunity from May 4, 2022 to June 13, 2022.
Each field was observed for approximately 10 ksec.
Seven targets (labeled 1-6, 11) are within the TeV emitting region and four (labeled 7-10) surround the Cocoon to provide a reference for the X-ray background.
The fields are shown in Figure <ref> with infrared, radio, and TeV observations of the Cocoon and summarized in Table <ref>.

§.§ Data Processing

We retrieved event lists and exposure maps from the HEASARC <cit.>.
Each field was observed multiple times, so we produced an image by combining each of the event lists using XSelect.
We ran a point source detection algorithm in XImage to search for point sources with signal to noise greater than three.
For each point source found, we excluded a circle with 30$''$ radius around it from further processing.
We also excluded detector regions affected by stray light (see Section <ref>).
Fields 2, 3, and 9 were affected by stray light from nearby X-ray sources.
After stray light removal, Field 9 had too few counts to provide a useful constraint so we do not use it in the following analysis.
We extract a spectrum from the entire remaining detector after point source and stray light filtering.
We produced an exposure map for each field by summing the individual exposure maps.
We produced calibration files using the XRT tools with the extended source option enabled in the XRT xrtmkarf task.

Spectrum from Field 11. The data are shown with error bars overlaid with the overall model (black line). The model components are the Solar Wind Charge Exchange (SWCX; pink dot-dot-dot-dash), Local Hot Bubble (LHB; red solid), Cold Galactic Ridge X-Ray Background (GRXB; green dash), hot GRXB (cyan dot), and Cosmic X-Ray
Background (CXB; blue dot-dash). The non-X-ray background is plotted with error bars (black crosses) and dominates above 6 keV.

§.§ Stray Light Filtering

Swift-XRT observations can be contaminated by stray light from X-ray sources outside the telescope field of view <cit.>.

Because we are searching for a source which is extracted from the entire detector, stray light contamination can strongly affect the spectrum.
In unprocessed images of Fields 2, 3 and 9, the stray light appears in a characteristic ring-like pattern.
The centers of Fields 2 and 3 are approximately 45$'$ from Cyg X-3 which is likely the source of the stray light which contaminates up to half of the detector for those observations.
The stray light in those fields is temporally variable with a hard power law spectrum consistent with a high mass X-ray binary like Cyg X-3.
The source of stray light in Field 9 is unknown, but it features a similar hard power law spectrum suggesting another X-ray binary.
We excluded contaminated detector regions in fields 2 and 3, and produced spectra from the remaining detector area.
As noted above, field 9 was too contaminated to produce a useful spectrum.
There are ten observations with useful data (1-8, 10, 11).
Below, “all fields" refers to fields 1-8, 10, and 11.

§ SPECTRAL FITTING

§.§ Modelling

We include four components in the astrophysical background model: the Galactic Ridge X-Ray Background (GRXB), the Local Hot Bubble (LHB), Cosmic X-Ray Background (CXB), and the Solar Wind Charge Exchange (SWCX). The model for the Cocoon is an absorbed power law. These are summarized below:

* GRXB: The GRXB is modeled with two absorbed Astrophysical Plasma Emission Code (APEC) <cit.> plasmas following previous studies of the GRXB <cit.>.
We refer to these as the “hot” and “cold” GRXB.
Both components are absorbed by the tbabs or tbgrain component.
The column density of the cold GRXB is allowed to float in the fits while the hot GRXB absorption column density is frozen to the Galactic value (see Sec <ref>)

* LHB: The LHB is modeled by an unabsorbed APEC plasma which dominates the spectra below 0.5 keV.

* CXB: The CXB is modeled by an absorbed power law with normalization and spectral index frozen to the results of <cit.>.
The absorption column density is frozen to the Galactic column density.

* SWCX: We model the SWCX with a single line at 0.5 keV.
In reality, the SWCX features many lines, but the Swift-XRT energy resolution at low energy is broad enough that the single line is adequate.

* Cocoon: A hypothetical Cocoon component is modelled by an absorbed power law.
Over the fit range 0.3–6.0 keV, a power law adequately represents the expected synchrotron emission.
The spectral index was frozen to 2.5 motivated by the HAWC TeV result and the neutral hydrogen column density was frozen to the nearby column density (see Sec <ref>). We also tested a spectral index of 2.0.

We include a non-X-ray background spectrum using 43 ksec of data collected with the telescope sun shutter closed.
This spectrum and its response is included in XSpec as a background for all the spectra.
The non-X-ray background rate is approximately $2\times 10^{-8} \rm~counts\,s^{-1}\,keV^{-1}\,pixel^{-1}$.
Each of the observations include up to 38000 pixels (with 2.36” pixel scale), which equates to a background rate $0.007 \rm~counts\,s^{-1}\,keV^{-1}$.
The spectra are dominated by non-X-ray background above 6 keV.
A fit for one field is shown in Figure <ref>.

§.§ $N_\mathrm{H}$ Column Density

Atomic and molecular hydrogen  column densities derived from visual extinction of stars and kinematic separation of 21-cm and CO lines combined with $\gamma$-ray data. The gray dashed line would be a perfect correspondence. The molecular hydrogen column density is the column density of hydrogen atoms that are in molecules.

The choice of column density, particularly the column density of nearby gas absorbing emission from the Cygnus region, affects the upper limit of the intrinsic emission from the Cocoon.
To quantify this impact we estimated the column density two ways.
First, we estimate the absorption column densities from 3-dimensional dust maps by <cit.> using the  software package <cit.>.
We extract the reddening, E(g-r), at 1.4 kpc and 10 kpc for each of our fields and convert it to E(B-V).
We convert E(B-V) to the neutral hydrogen column density using the relation from <cit.>.
Second, we use the column density maps derived by <cit.>.
The latter uses the Doppler shift of the 21-cm and CO emission lines to kinematically separate HI and H$_2$ into two regions: the Local Arm (where the Cygnus region and Cocoon are embedded), and beyond.
Additional visual extinction suggests additional dark neutral gas which is included in the nearby H$_2$ map. The conversion between CO emission and visual extinction to column density is based on $\gamma$-ray data.
The column densities are compared in Figure <ref>.
The average nearby kinematic separation column density is 2.1 times higher than the dust map result, and the average total Galactic column density is 1.4 times larger suggesting larger attenuation within the Cocoon.
<cit.> suggest an H$_2$ fraction of 20% which is assumed in the  model;
however, the kinematic separation finds that the average fraction of hydrogen atoms in the molecular phase is approximately 50% for this region.

§.§ Hypothesis Testing

We performed a log-likelihood ratio test (LRT) for the models in the following sections where the null hypothesis is that there is no additional diffuse component and the test hypothesis is that there is additional emission as described in item 5 in Section <ref>.
Each model was fit by minimizing the C-statistic (equivalently maximizing the Poisson likelihood) with XSpec <cit.> between 0.3–6.0 keV.
The models are nested, but the fit parameter space is constrained to keep the Cocoon normalization positive or zero.
The test statistic (twice the log-likelihood ratio) thus follows a “$\chi^2/2$" distribution: $\rm~CDF(TS) =  {1}/{2}+\int_0^{TS} dx {\chi^2(x)}/{2} $ for $\rm~TS > 0$ <cit.>.
We verified that the test statistic does follow this distribution using Monte Carlo simulations where we generated synthetic data sets based on the best fit background model and performed a likelihood ratio test on each data set.
The resulting test statistic distribution had a $ {\chi^2}/{2}$ distribution with half the trials having $\rm~TS = 0$.
Intuitively, the unconstrained Cocoon normalization would be normally distributed around zero.
In the constrained case, the negative values are forced to be zero and thus have test statistic equal to zero.
For significance lower than $3\sigma$, upper limits for the parameters are defined by the 90% confidence upper confidence interval on the parameter fit.

§.§ Individual Fields

Intensity for each field between 2–10 keV after subtracting CXB compared with the position of the center of each observation. Unlike <cit.>, we do not find a monotonic relationship between Galactic latitude and hard X-ray emission.

We fit each field with a background and a Cocoon plus background model.
The plasma temperatures, normalizations, and absorption column densities are allowed to vary in both fits.
In the latter case, the normalization of the Cocoon component is allowed to vary.
The best fit parameters for the background-only model are tabulated in Table <ref>.
In general, the temperatures are similar except where the fit is poorly constrained by limited statistics at higher energy.
The background fields (7-10) in particular poorly constrain the hot GRXB.
The 2–10 keV intensity after subtracting CXB for each field is shown in Figure <ref>.
Unlike <cit.>, we find no monotonic relationship between the intensity and Galactic latitude.
Our Field 3 nearly overlaps the Suzaku Src 1 field and we find similar intensity in these two fields above 2 keV.
Field 3 and 6 are the nearest to Suzaku Src 2 but are 48$'$ and 65$'$ away.
The hard band intensity in Src 2 is similar to the intensity in Field 6.
We performed the hypothesis test described in <ref> for an additional Cocoon component in the on-source fields modeled by an absorbed power law.
No field has a Cocoon component with significance greater than $3\sigma$.
The 90% confidence upper limits on the Cocoon power law normalization are shown in Table <ref>.
The mean upper limit is $2.0\times 10^{-8}\rm~erg\,cm^{-2}\,s^{-1}\,sr^{-1}$ between 2–10 keV.

§.§ Integrated Limit on the Cygnus Cocoon

In addition to the individual fits, we jointly fit the spectra from all the fields assuming that a potential X-ray signal is proportional to the integrated 1-100 TeV surface brightness, $F_X\propto F_\gamma$.
For the fields outside the Cocoon, $F_X$ of the Cocoon component is set to 0.
For the fields that have a non-zero HAWC flux, we tie the normalization of the Cocoon component to the Cocoon normalization in Field 1 multiplied by the ratio of the HAWC flux in each field divided by the HAWC flux in Field 1.
The model adds the constraint that the GRXB and LHB temperatures are the same across all the fields within the Cocoon region.
The 90% confidence upper limit on the Cocoon normalization parameter is $6.7\times 10^{-4} \rm~cm^{-2}\,s^{-1}\,keV^{-1}$.
Multiplying this normalization parameter by the integrated TeV flux over the entire Cocoon (divided by the TeV flux in Field 1) and using the dust map based column density, we find the total flux

\begin{equation}
F = \int dE E \frac{dN}{dE\,dA\,dt} < 8.7\times 10^{-11}\,\rm erg\,cm^{-2}\,s^{-1}
\end{equation}

between 2–10 keV when assuming $dN/dE \propto E^{-2.5}$.
This limit is shown in Figure <ref>.
Using a spectral index of 2 increased the 2–10 keV flux to $1.0~\times~10^{-10}~\rm erg\,cm^{-2}\,s^{-1}$.
The statistical significance of the likelihood ratio tests for $E^{-2.5}$ and $E^{-2}$ were $2.3\sigma$ and $3.4\sigma$, respectively.

We repeated the test using the column densities based on kinematic separation and calibrated on $\gamma$-ray data.
We also replaced the  model with  which allows the user to set the fraction of hydrogen atoms in the molecular phase.
The significance of the log-likelihood ratio test reduces to $1.6\sigma$ and $1.7\sigma$ for $E^{-2.5}$ and $E^{-2.0}$ spectra respectively.
The corresponding 2–10 keV upper limits are $5.2\times 10^{-11}~\rm erg\,cm^{-2}\,s^{-1}$ and $4.4\times 10^{-11}~\rm erg\,cm^{-2}\,s^{-1}$
The results of both models are summarized in Table <ref>.

Multi-wavelength spectral energy distribution with data and flux limits from HAWC <cit.>, Fermi-LAT<cit.>, Suzaku-XIS <cit.>, Swift-XRT, and CGPS <cit.>. The red shaded region bounds the X-ray upper limits in this work for the two hydrogen column density models considered. The black curves represent the maximum synchrotron (thick black) and inverse Compton (thin black) emission allowed by our limit obtained using the kinematic separation column density. The injection electron spectrum is assumed to be $dN_e/dE_e\propto E_e^{-2}$ and magnetic field is set to 20 $\mu$G. The grey dash-dotted curve show the Bremsstrahlung emission by these electrons. The grey dashed curve indicates the pion decay emission by protons. The grey solid curve shows the sum of the emission components. Details of the background radiation model may be found in <cit.>.

§ DISCUSSION

The X-ray upper limits derived from the Swift-XRT observations constrain the electron contribution to gamma rays above 1 TeV.
We adopt the leptonic model presented in <cit.> which assumes that a population of relativistic electrons are continuously injected by the star clusters, diffuse in the Cocoon magnetic field, and inverse Compton scatter the stellar light fields and dust radiation field.
When fixing the magnetic field strength to 20 $\mu$G, which is estimated using the gas pressure based on infrared imaging <cit.>, the integrated Swift-XRT limit constrains the leptonic $\gamma$-ray emission at 1 TeV to $25\%$ of observed emission.

Alternatively, when freeing the magnetic field and allowing a dominant leptonic contribution above 10 TeV, the X-ray observation constrains the synchrotron radiation by the electrons and limits the field strength to $B\lesssim 10\,\mu$G. The maximum electron energy that can be accelerated by the stellar winds in this field is $E_{e,\rm max} \approx 34\,\eta^{1/2} (B/10\,\mu{\rm G})^{-1/2}(v_w/10^3\,\rm km\,s^{-1})\,\rm TeV$ <cit.>, where $\eta = \delta B^2/B^2$ denotes the level of magnetic-field fluctuations and $v_w\sim10^3\,\rm km\,s^{-1}$ is the velocity of stellar winds. These electrons cannot explain the photons with energies above  $E_{e,\rm max}$.

Although the spectrum in Field 1 cannot rule out the purely background hypothesis, it has significantly brighter hard emission than the rest of the field.
<cit.> searched for a X-ray counterparts to  GeV sources with ASCA including GeV J2035+4214 which lies within the Cygnus Cocoon.
Their observation of this source overlaps with Field 1 where they found three X-ray sources.
One of these sources (Src1) is clearly extragalactic because it is strongly absorbed and has a spectral index of 1.41.
This source was detected by Swift-XRT and excluded during data processing.
The other two sources (Src2 and Src3) have absorption column densities equal to 0.98 and 0.84 which suggest distances consistent with the Cygnus Cocoon.
These two sources have spectral indices equal to 2.02 and 2.44 and the latter is identified as an extended source.
In Chandra imaging of the same region, the extended source is in fact several confused point sources <cit.>.
We did not detect Src2 in Swift-XRT images.
A comparison of spectra including and excluding Src3 made no apparent change in the shape of the Swift-XRT spectra.

§ SUMMARY

We observed 11 fields in and near the Cygnus Cocoon with Swift-XRT to search for non-thermal X-ray emission.
The X-ray emission for each field is well described by a background model consisting of the Local Hot Bubble, Galactic Ridge X-Ray Background, and Cosmic X-Ray Background.
We establish upper limits for an additional power law component for each of the fields within the Cygnus Cocoon.
We also tested the hypothesis that non-thermal X-ray emission traces the TeV emission measured by HAWC and established an upper limit to this component. The limits thus obtained are sensitive to the hypotheses made to model X-ray absorption by interstellar matter, the impact of which is quantified.
The X-ray upper limits constrain the contribution of inverse Compton scattering to the Cocoon's TeV emission.
The X-ray limits can be explained by a magnetic field below equipartition within the Cocoon or by reducing the relativistic electron density.
In the latter case, the inverse Compton emission can explain up to one quarter of the TeV emission.

We thank Binita Hona for discussion on the HAWC significance map. The work of D.G. and K.F. is supported by the Office of the Vice Chancellor for Research and Graduate Education at the University of Wisconsin-Madison with funding from the Wisconsin Alumni Research Foundation. D.G. and K.F. acknowledge support from NASA (NNH20ZDA001N-Swift) and National Science Foundation (PHY-2110821). A.B. and K.L.P. acknowledge funding from the UK Space Agency.

The work of K.F and F.H is supported by the Office of the Vice Chancellor for Research and Graduate Education at the University of Wisconsin-Madison with funding from the Wisconsin Alumni Research Foundation. K.F. acknowledges support from NASA (NNH19ZDA001N-Fermi, NNH20ZDA001N-Fermi, NNH20ZDA001N-Swift) and National Science Foundation (PHY-2110821). J.S.G. thanks the University of Wisconsin College of Letters and Science for partial support of his IceCube-related research. The research of F.H was also supported in part by the U.S. National Science Foundation under grants PLR-1600823 and PHY-1913607.


# Spin-Orbit Coupling in Transition Metal Dichalcogenide Heterobilayer Flat
Bands

Louk Rademaker  Department of Theoretical Physics, University of Geneva,
CH-1211 Geneva, Switzerland

###### Abstract

The valence flat bands in transition metal dichalcogenide (TMD) heterobilayers
are shown to exhibit strong intralayer spin-orbit coupling. This is reflected
in a simple tight-binding model with spin-dependent complex hoppings based on
the continuum model. A perpendicular electric field causes interlayer
hybridization, such that the effective model is equivalent to the Kane-Mele
model of topological insulators. The proposed model can be used as a starting
point to understand interactions and the experimentally observed topological
phases.

Heterostructures with Moiré patterns, emerging from the lattice mismatch or
misalignment of two-dimensional materials, have proven to be a fertile ground
for novel electronic states.Andrei:2020asg ; Carr:2020cm ; Balents:2020.725B
Amongst possible materials, ‘flat bands’ are expected to occur in bilayers of
transition metal dichalcogenide (TMD) monolayer semiconductors. TMD
heterobilayers, when the two layers are of a different material, have been
proposed as an ideal realization of the triangular lattice Hubbard model with
real short-range hoppings.Wu:2018ic

Indeed, experimental results in WS2/WSe2 and MoTe2/WSe2 bilayers include the
observation of a Mott insulating state at half-filling of the flat
bandsTang:2020bb ; Li:2021cd  and generalized Wigner crystals at other
fillingsRegan:2020fk ; Xu:2020dx ; Huang:2021io ; Jin:2021es . Recently, it
came as quite a surprise that upon tuning a perpendicular electric field, both
the quantum spin Hall effect (QSH) and quantum anomalous Hall (QAH) effect
were observed in bilayer MoTe2/WSe2.Li:2021vy  After all, the original
heterobilayer modelWu:2018ic  did not include any spin-orbit coupling, in
contrast to the homobilayer model that predicted topological
bandstructures.Wu:2019gm

There is a growing literature aiming to understand and characterize TMD
heterobilayers,Shih:2020ia ; Weston:2020fz ; Zhang:2020hk ; Shabani:2021fh ;
Zhang:2019wk ; Zhang:2021hx ; Padhi:2021bs ; Vitale:2021vg ; Li:2021up ;
Li:2021cc ; Naik:2020bq ; Naik:2018eha ; Pan:2021uh ; Pan:2020cd ; Pan:2020kga
where the observed topological transitions are proposed to originate in
lattice relaxationXie:2021ay  or interlayer coupling.Zhang:2021mote2 ;
2021arXiv211101152P  However, as of yet there is no simple tight-binding model
describing topological phases.

In this Letter, I reinterpret the continuum model of Ref. Wu:2018ic , using
symmetry arguments related to the backfolding of the monolayer momenta to the
Moiré mini-Brillouin zone. This directly yields strong spin-orbit coupling in
the form of a complex spin-dependent hopping. The corresponding tight-binding
model is equivalent to the next-nearest neighbor hopping in the Kane-Mele (KM)
model. By tuning the band offset through a perpendicular electric field, the
second layer comes into play and we find a full realization of the KM model
with a corresponding topological transition and band inversion to a
topological insulator phase. The tight-binding model derived in this Letter
can serve as a starting point for studying the interplay between topology and
interaction effects in TMD heterobilayers.

Continuum model –  Let me first recap the essence of the continuum model of
flat bands in TMD heterobilayers.Wu:2018ic  Monolayer transition metal
dichalcogenide (MX2 with M=W, Mo and X = S, Se, Te) have a hexagonal lattice
with $C_{3}$ rotational symmetry. They are semiconductors with a relatively
large direct bandgap at the ${\bf K}$ and ${\bf K^{\prime}}$ point. The
valence band states have a large Ising spin splitting, such that states at
${\bf K}$ are spin-$\uparrow$ and the states at ${\bf K^{\prime}}$ are
spin-$\downarrow$. Spin ($\uparrow$ or $\downarrow$) and valley (${\bf K}$ or
${\bf K^{\prime}}$) are therefore inextricably coupled.

Whenever two different TMD monolayers are combined in a heterobilayer with
twist angle $\theta$, a Moiré pattern is created with lengthscale
$a_{M}=1/\sqrt{\frac{1}{a_{1}^{2}}+\frac{1}{a_{2}^{2}}-\frac{2\cos\theta}{a_{1}a_{2}}}$
where $a_{1,2}$ are the lattice constants of the monolayers. Because different
TMD monolayers have different band gaps, the top of the valence band of one
layer lies in the gap of the other (a type-I or type-II band alignment).
Consequently, the valence band in a heterobilayer exclusively consists of
states localized in a single layer.

These heterobilayer valence band states can be described, per spin/valley, by
imposing a Moiré potential $V({\bf r})$ on the monolayer valence band states,

$\displaystyle H$ $\displaystyle=$ $\displaystyle-\frac{\hbar^{2}{\bf
Q}^{2}}{2m^{*}}+V({\bf r})$ $\displaystyle V({\bf r})$ $\displaystyle=$
$\displaystyle\sum_{{\bf g}_{j}}V_{j}\exp\left[i{\bf g}_{j}{\bf r}\right]$ (1)

where ${\bf
g}_{j}=\frac{4\pi}{\sqrt{3}}(-\sin\frac{2\pi(j-1)}{6},\cos\frac{2\pi(j-1)}{6})$
are the six reciprocal Moiré vectors and ${\bf Q}$ is the momentum relative to
the ${\bf K}/{\bf K^{\prime}}$ point of the monolayer. Because $\Delta({\bf
r})$ must be real and the $C_{3}$ symmetry, we have $V_{1}=V_{3}=V_{5}$ and
$V_{2}=V_{4}=V_{6}=V_{1}^{*}$,Wu:2018ic ; 2018PhRvB..97c5306W ;
2017PhRvL.118n7401W  which means one can parametrize the Moiré potential using
only two parameters $(V,\psi)$ such that $V_{1}=Ve^{i\psi}$. For MoSe2/WSe2
bilayers, Ref. Wu:2018ic  calculated $(V,\psi)=($6.6 meV$,-94^{\circ})$. This
results in a topmost valence flat band with almost perfect Gaussian Wannier
orbitals centered on a Moiré triangular lattice.

Figure 1: a) When two layers with inequivalent lattice constant are combined,
in general the ${\bf K}$-points of the monolayers map onto the
${\bm{\kappa}}$-points of the mini-Brillouin zone. This is specifically the
case for aligned WS2/WSe2 and MoTe2/WSe2. b) In case of an incommensurate
Moiré pattern, for example due to a twist angle $\theta$, the natural mapping
is similar. c) Only when the Moiré length satisfies $a_{M}=pa_{1}$ with $p$
mod $3=0$, the ${\bf K}$-points of the layer 1 map onto the center of the
mini-BZ ${\bm{\gamma}}$, as was done in Ref. Wu:2018ic . Figure 2: a) Flat
valence band structure of aligned AA-stacked WS2/WSe2. The colors red/blue
indicate the spin $\downarrow$/$\uparrow$ of the bands, showing significant
spin-orbit coupling. The thin lines indicate the monolayer valence band states
in the absence of the Moiré potential. Upon wannierization, the nearest
neighbor hopping equals $|t_{1}|=1.8$ meV. b) Flat valence band structure of
aligned AB-stacked MoTe2/WSe2, with a similar spin-orbit splitting. The
nearest neighbor hopping equals $|t_{1}|=4.0$ meV.

Momentum backfolding –  The continuum model is expressed using the momenta of
the monolayer, here denoted in capital letters ${\bf K}_{\ell}$, ${\bf
K^{\prime}}_{\ell}$ with layer index $\ell$. The Moiré potential folds the
momenta back to a small mini-Brillouin zone (BZ) whose high-symmetry points
are ${\bm{\gamma}}=(0,0)$, ${\bm{\kappa}}=(\frac{4\pi}{3a_{M}},0)$,
${\bm{\kappa}^{\prime}}=(\frac{2\pi}{3a_{M}},\frac{2\pi}{\sqrt{3}a_{M}})$ and
${\bm{\mu}}=\frac{1}{2}{\bf g}_{1}$. The question is which monolayer momenta
correspond to which mini-BZ momenta.

In Ref. Wu:2018ic , the authors map the monolayer ${\bf K}$ to the
${\bm{\gamma}}$ point of the mini-BZ. Such backfolding only occurs when the
Moiré length $a_{M}$ is a multiple of three times the monolayer lattice
constant, see Fig. 1c. In general, however, this is not true: the Moiré length
in aligned commensurate cases can be expressed as $a_{M}=pa_{1}=qa_{2}$ for
coprime integers $p,q$. When both $p,q$ are not multiples of three, the ${\bf
K}$-points of the single layers map onto the ${\bm{\kappa}}$-points of the
mini-BZ. This is the case for aligned WS2/WSe2 with
$a_{M}=25a_{\mathrm{WSe}_{2}}=26a_{\mathrm{WS}_{2}}$Li:2021cc  and aligned
MoTe2/WSe2 with
$a_{M}=13a_{\mathrm{MoTe}_{2}}=14a_{\mathrm{WSe}_{2}}$.Zhang:2021mote2
Therefore the relevant momentum backfolding to describe the experiments is
${\bf K}_{1},{\bf K^{\prime}}_{2}\rightarrow{\bm{\kappa}}$ and ${\bf
K}_{2},{\bf K^{\prime}}_{1}\rightarrow{\bm{\kappa}^{\prime}}$, as is shown in
Fig. 1a. For general incommensurate and possibly twisted Moiré patterns we
retain this mapping of momenta as shown in Fig. 1b. Note that in the continuum
theory of TMD homobilayers the same mapping of momenta is used.Wu:2019gm

Figure 3: a) The effective model of a plain heterobilayer contains orbitals at
triangular lattice sites with nearest neighbor hopping
$t_{1}=|t_{1}|e^{i\phi}$ where the phase is spin-dependent,
$\phi=\frac{2\pi}{3}\sigma^{z}$. The states live completely within one
monolayer. In the case of MoTe2/WSe2, the states are centered at MM’ sites in
the MoTe2 layer. b) By tuning the vertical displacement field the second layer
becomes relevant. The second layer, in the case of AB-stacked MoTe2/WSe2, is
described by the same triangular lattice model with opposite spin-orbit
coupling phases and a different Moiré unit cell center XX’. As a result, the
combination of in-plane spin-orbit coupling Eq. (3) and interlayer hopping Eq.
(4) leads to an ideal realization of the Kane-Mele model.

With the correct momentum backfolding, the valence band structure of aligned
AA-stacked WS2/WSe2 and AB-stacked MoTe2/WSe2 are shown in Fig. 2. For
WS2/WSe2, I use $a_{M}=7.98$ nm, $m^{*}=0.36m_{e}$, and the Moiré potential
$(V,\psi)=(7.7$ meV, $-106^{\circ})$ is calculated using density functional
theory with Quantum ESPRESSO2017JPCM…29T5901G ; 2009JPCM…21M5502G  with a
Coulomb cut-offSohier:2017fl  and the method of Ref. Wu:2018ic . For
MoTe2/WSe2, I use the parameters $a_{M}=4.55$ nm, $m^{*}=0.65m_{e}$, and
$(V,\psi)=(7$ meV, $-14^{\circ})$ based on Ref. 2021arXiv211101152P . The
resulting flat bands have a very clear spin-orbit splitting at ${\bm{\kappa}}$
and ${\bm{\kappa}^{\prime}}$. The general structure of spin/valley-split bands
except along the line $\gamma-\mu$ was also observed in large-scale DFT
calculations of WS2/WSe2Li:2021cc  and MoTe2/WSe2Zhang:2021mote2 .

Tight-binding model – The tight-binding model for the valence flat bands can
be derived using symmetry arguments. I will focus only on nearest neighbor
hopping $t_{1}$ since it is much larger than the longer-ranged hoppings
$t_{2},t_{3}$. The $C_{3}$ symmetry of the heterobilayer implies that $t_{1}$
can be complex, $t_{1}=|t_{1}|e^{i\phi}$, as shown in Fig. 3a. Within the
continuum model of Eq. (1), the band-structure is six-fold symmetric around
the top of the flat band. Given that the top of the flat band is positioned at
${\bm{\kappa}}$ or ${\bm{\kappa}^{\prime}}$, it directly follows that the
complex phase equals $\phi=\pm\frac{2\pi}{3}$, respectively. The resulting
tight-binding model is thus a triangular lattice with spin-orbit coupled
hopping

$H=t_{1}\sum_{\langle ij\rangle\sigma}e^{i\phi\sigma^{z}\nu_{\langle
ij\rangle}}c^{\dagger}_{i\sigma}c^{\phantom{\dagger}}_{j\sigma}$ (2)

where $\nu_{\langle ij\rangle}=\pm 1$ depending on the direction of the bond,
$\sigma^{z}=\pm 1$ indicates the $z$-component of the electron spin, and
$\phi=\frac{2\pi}{3}$. What is left is to derive the magnitude of the hopping
$|t_{1}|$, and the positioning of the triangular lattice sites, both of which
depend on the specific choice of heterobilayer compounds.

Because a honeycomb lattice is equivalent to two intertwined triangular
lattices, we see that Eq. (2) corresponds for a single spin to the hopping
model introduced by HaldaneHaldane:1988gh  for the Chern insulator. Including
both spins, one find that the Kane-Mele model of topological
insulators,Kane:2005hl ; Kane:2005gb  restricted to one honeycomb sublattice,
is equal to Eq. (2) with $\phi=\pi/2$.

Figure 4: a. The flat valence bands in the MoTe2/WSe2 heterobilayer for the
${\bf K}$-valley (the opposite valley is shown in thin lines). The top flat
band contains exclusively states in the MoTe2 layer, whereas the bottom band
contains exclusively states in the WSe2 layer. b. By tuning the perpendicular
electric field $V$ through a critical value, a band inversion occurs around
the ${\bm{\kappa}}^{\prime}$ point. The colorscale indicates whether states
live in the MoTe2 layer (red) or WSe2 layer (blue). The resulting bands carry
a nonzero Chern number with opposite Chern numbers for opposite spin. This
transition to a topological insulator bandstructure is described by the Kane-
Mele model of Eqs. (3)-(5).

Topological transition in MoTe2/WSe2 –  It has been argued that the observed
topological transition in MoTe2/WSe2Li:2021vy  is due to a coupling to the
states in the second layer,Zhang:2021mote2 ; 2021arXiv211101152P  induced by a
perpendicular electric field. There is, however, not yet a simple tight-
binding model explaining this topological transition. Building on the model of
Eq. (2), we can show how topologically nontrivial bands arise in a hopping
model.

In the absence of an electric field, the valence flat band of MoTe2/WSe2 is
localized in the MM’ region of the MoTe2 layer. In the other layer, the Moiré
potential localizes the top of the valence band at the XX’ region.
Furthermore, because the valence band top in the ${\bf K}$-valley of WSe2 maps
onto ${\bm{\kappa}^{\prime}}$, we find that the effective tight-binding model
for the WSe2 states has opposite phase $\phi=\mp\frac{2\pi}{3}$. The resulting
two flat bands can therefore be represented by a honeycomb lattice as is shown
in Fig. 3b. The intralayer hopping in the MoTe2/WSe2 bilayer is represented by
the spin-orbit coupling term

$H=\sum_{\langle\langle
ij\rangle\rangle\ell\sigma}t_{\ell}e^{i\ell\sigma^{z}\nu_{\langle\langle
ij\rangle\rangle}\phi}c^{\dagger}_{i\ell\sigma}c^{\phantom{\dagger}}_{j\ell\sigma}$
(3)

where $\langle\langle ij\rangle\rangle$ represents the next-nearest neighbor
on the honeycomb lattice where $\nu_{\langle\langle ij\rangle\rangle}=\pm 1$
depends on the direction, $\ell=\pm 1$ the layer index, and $\sigma^{z}=\pm 1$
the spin. This is exactly the Kane-Mele spin-orbit coupling term.Kane:2005hl ;
Kane:2005gb

The interlayer hopping is now given by nearest neighbor hopping on the
honeycomb lattice. The $C_{3}$ symmetry once again constrains the possible
complex phases of the interlayer hopping, as described in Ref. Zhang:2021mote2
. For the tight-binding model this yields

$H_{\perp}=t_{\perp}\sum_{\langle
ij\rangle\sigma}e^{i\frac{2\pi}{3}\nu_{\langle
ij\rangle}}c^{\dagger}_{i\ell\sigma}c^{\phantom{\dagger}}_{j\overline{\ell}\sigma}$
(4)

where $\langle ij\rangle$ now couples the nearest neighbors on the honeycomb
lattice. The parameter $\nu_{\langle ij\rangle}=0,1$ or $2$ depends on the
direction and increases counterclockwise when going around the MoTe2 lattice
sites. Eqs. (3)-(4) form an effective tight-binding model constrained by
symmetry and is visualized in Fig. 3b.

In the absence of an electric field the flat bands from the WSe2 layer are at
a much lower energy than the MoTe2 states. A perpendicular electric field $V$
can shift the WSe2 states upward,

$H_{V}=(V+\Delta)\sum_{i}n_{i,\ell=2}$ (5)

where $\Delta<0$ is the band offset. By increasing $V$ the bands from the two
layers will overlap, at which point the interlayer coupling Eq. (4) becomes
relevant. With $t_{2}=3.4$ meV extracted from the effective mass of the WSe2
valence band and $t_{\perp}=4$ meV, we calculate the bandstructure as a
function of perpendicular field in Fig. 4. Above the critical field value
$V_{c}$, a band inversion happens at the ${\bm{\kappa}^{\prime}}$ point for
the states in the ${\bf K}$-valley. For $V>V_{c}$, the top valence flat bands
obtain a nonzero Chern number which is opposite for the two spin
species.Fukui:2005cn  As such, the system has become a topological insulator
fully described by the honeycomb lattice model of Kane-Mele.

The effective Kane-Mele model applies to all TMD heterobilayers provided the
applied perpendicular electric field is sufficient to overcome the valence
band offset.

Interactions –  The presence of spin-orbit coupling has some implications for
the possible interacting states. In particular, at half-filling of the topmost
flat band the complex hopping phases lead to an effective spin model with
Dzyaloshinskii-Moriya (DM) interactions.Pan:2020kga  Its strength,
characterized by the hopping phase $\phi=\frac{2\pi}{3}$, stabilizes an in-
plane 120∘ antiferromagnetic order. At the same time, Chern bands at half-
filling are known to be susceptible to full spin polarization leading to a
Quantum Anomalous Hall (QAH) effect.Bultinck:2019wp ; Rademaker:2020fr  The
Mott-to-QAH transition observed in Ref. Li:2021vy  can possibly be understood
as a metamagnetic transition from in-plane Néel to Ising ferromagnetic order.

Note that within the current model, there is no band gap between the
topologically nontrivial bands. However, interaction-driven renormalization of
the bands can open up the gap at full filling of the top valence flat
bands.2021arXiv211101152P

Conclusion – I showed that by correctly mapping the monolayer momenta onto the
mini-BZ, the effective model of TMD heterobilayers obtains a strong spin-orbit
coupling as given by Eq. (2). The interlayer coupling can be described by a
Kane-Mele model following Eqs. (3)-(5), yielding a topological transition as a
function of perpendicular electric field.

Having an effective tight-binding model is an important step towards a full
understanding of possible strongly correlated phases. This Letter shows that
spin-orbit coupling cannot be ignored in further studies of TMD
heterobilayers.

I thank Johannes Motruk, Kin Fai Mak, Fengcheng Wu, Allan MacDonald, and
Vladimir Dobrosavljević for useful discussion. I acknowledge support by Swiss
National Science Foundation via an Ambizione grant PZ00P2_174208.

## References

  * (1) E. Y. Andrei and A. H. MacDonald,  _Graphene bilayers with a twist_, Nat. Mater. 19, 1265 (2020).
  * (2) S. Carr, S. Fang, and E. Kaxiras,  _Electronic-structure methods for twisted moiré layers_, Nature Reviews Materials 556, 1 (2020).
  * (3) L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young,  _Superconductivity and strong correlations in moiré flat bands_, Nat. Phys. 16, 725 (2020).
  * (4) F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald,  _Hubbard Model Physics in Transition Metal Dichalcogenide Moiré Bands_, Phys. Rev. Lett. 121, 026402 (2018).
  * (5) Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, and K. F. Mak,  _Simulation of Hubbard model physics in WSe2/WS2 moiré superlattices_, Nature 579, 353 (2020).
  * (6) T. Li, S. Jiang, L. Li, Y. Zhang, K. Kang, J. Zhu, K. Watanabe, T. Taniguchi, D. Chowdhury, L. Fu, J. Shan, and K. F. Mak,  _Continuous Mott transition in semiconductor moiré superlattices_, Nature 597, 350 (2021).
  * (7) E. C. Regan, D. Wang, C. Jin, M. I. Bakti Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang, K. Yumigeta, M. Blei, J. D. Carlström, K. Watanabe, T. Taniguchi, S. Tongay, M. Crommie, A. Zettl, and F. Wang,  _Mott and generalized Wigner crystal states in WSe2/WS2 moiré superlattices_, Nature 579, 359 (2020).
  * (8) Y. Xu, S. Liu, D. A. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, V. Elser, K. F. Mak, and J. Shan,  _Correlated insulating states at fractional fillings of moiré superlattices_, Nature 587, 214 (2020).
  * (9) X. Huang, T. Wang, S. Miao, C. Wang, Z. Li, Z. Lian, T. Taniguchi, K. Watanabe, S. Okamoto, D. Xiao, S.-F. Shi, and Y.-T. Cui,  _Correlated insulating states at fractional fillings of the WS2/WSe2 moiré lattice_, Nat. Phys. 17, 715 (2021).
  * (10) C. Jin, Z. Tao, T. Li, Y. Xu, Y. Tang, J. Zhu, S. Liu, K. Watanabe, T. Taniguchi, J. C. Hone, L. Fu, J. Shan, and K. F. Mak,  _Stripe phases in WSe2/WS2 moiré superlattices_, Nat. Mater. 20, 940 (2021).
  * (11) T. Li, S. Jiang, B. Shen, Y. Zhang, L. Li, T. Devakul, K. Watanabe, T. Taniguchi, L. Fu, J. Shan, and K. F. Mak,  _Quantum anomalous Hall effect from intertwined moiré bands_, arXiv (2021).
  * (12) F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. H. MacDonald,  _Topological Insulators in Twisted Transition Metal Dichalcogenide Homobilayers_, Phys. Rev. Lett. 122, 086402 (2019).
  * (13) E.-M. Shih, Y. Bai, K. Watanabe, T. Taniguchi, C. Dean, and C. E. Dreyer,  _Visualization of moiré superlattices_, Nature Nanotechnology 15, 580 (2020).
  * (14) A. Weston, Y. Zou, V. Enaldiev, A. Summerfield, N. Clark, V. Zólyomi, A. Graham, C. Yelgel, S. Magorrian, M. Zhou, J. Zultak, D. Hopkinson, A. Barinov, T. H. Bointon, A. Kretinin, N. R. Wilson, P. H. Beton, V. I. Fal’ko, S. J. Haigh, and R. Gorbachev,  _Atomic reconstruction in twisted bilayers of transition metal dichalcogenides_, Nature Nanotechnology 15, 592 (2020).
  * (15) Z. Zhang, Y. Wang, K. Watanabe, T. Taniguchi, K. Ueno, E. Tutuc, and B. J. LeRoy,  _Flat bands in twisted bilayer transition metal dichalcogenides_, Nat. Phys. 16, 1093 (2020).
  * (16) S. Shabani, D. Halbertal, W. Wu, M. Chen, S. Liu, J. Hone, W. Yao, D. N. Basov, X. Zhu, and A. N. Pasupathy,  _Deep moiré potentials in twisted transition metal dichalcogenide bilayers_, Nat. Phys. 17, 720 (2021).
  * (17) Y. Zhang, N. F. Q. Yuan, and L. Fu,  _Moiré quantum chemistry: Charge transfer in transition metal dichalcogenide superlattices_, Phys. Rev. B 102, 201115 (2020).
  * (18) Y. Zhang, T. Liu, and L. Fu,  _Electronic structures, charge transfer, and charge order in twisted transition metal dichalcogenide bilayers_, Phys. Rev. B 103, 155142 (2021).
  * (19) B. Padhi, R. Chitra, and P. W. Phillips,  _Generalized Wigner crystallization in moiré materials_, Phys. Rev. B 103, 125146 (2021).
  * (20) V. Vitale, K. Atalar, A. A. Mostofi, and J. Lischner,  _Flat band properties of twisted transition metal dichalcogenide homo- and heterobilayers of MoS 2, MoSe2, WS2 and WSe2_, arXiv (2021).
  * (21) H. Li, S. Li, M. H. Naik, J. Xie, X. Li, E. Regan, D. Wang, W. Zhao, K. Yumigeta, M. Blei, T. Taniguchi, K. Watanabe, S. Tongay, A. Zettl, S. G. Louie, M. F. Crommie, and F. Wang,  _Imaging local discharge cascades for correlated electrons in WS2/WSe2 moiré superlattices_, arXiv (2021).
  * (22) H. Li, S. Li, M. H. Naik, J. Xie, X. Li, J. Wang, E. Regan, D. Wang, W. Zhao, S. Zhao, S. Kahn, K. Yumigeta, M. Blei, T. Taniguchi, K. Watanabe, S. Tongay, A. Zettl, S. G. Louie, F. Wang, and M. F. Crommie,  _Imaging moiré flat bands in three-dimensional reconstructed WSe 2 /WS 2 superlattices_, Nat. Mater. , 1 (2021).
  * (23) M. H. Naik, S. Kundu, I. Maity, and M. Jain,  _Origin and evolution of ultraflat bands in twisted bilayer transition metal dichalcogenides: Realization of triangular quantum dots_, Phys. Rev. B 102, 075413 (2020).
  * (24) M. H. Naik and M. Jain,  _Ultraflatbands and Shear Solitons in Moiré Patterns of Twisted Bilayer Transition Metal Dichalcogenides_, Phys. Rev. Lett. 121, 266401 (2018).
  * (25) H. Pan and S. Das Sarma,  _Interaction-driven filling-induced metal-insulator transitions in 2D moiré lattices_, Phys. Rev. Lett. 127, 096802 (2021).
  * (26) H. Pan, F. Wu, and S. Das Sarma,  _Quantum phase diagram of a Moiré-Hubbard model_, Phys. Rev. B 102, 201104 (2020).
  * (27) H. Pan, F. Wu, and S. Das Sarma,  _Band topology, Hubbard model, Heisenberg model, and Dzyaloshinskii-Moriya interaction in twisted bilayer WSe2_, Phys. Rev. Research 2, 033087 (2020).
  * (28) Y.-M. Xie, C.-P. Zhang, J.-X. Hu, K. F. Mak, and K. T. Law,  _Theory of Valley Polarized Quantum Anomalous Hall State in Moire MoTe2/WSe2Heterobilayers_, arXiv , arXiv:2106.13991 (2021).
  * (29) Y. Zhang, T. Devakul, and L. Fu,  _Spin-textured Chern bands in AB-stacked transition metal dichalcogenide bilayers_, Proc. Natl. Acad. Sci. USA 118, e2112673118 (2021).
  * (30) H. Pan, M. Xie, F. Wu, and S. Das Sarma,  _Topological Phases in AB-stacked MoTe 2/WSe2: $\mathbb{Z}_{2}$ Topological Insulators, Chern Insulators, and Topological Charge Density Waves_, arXiv , arXiv:2111.01152 (2021).
  * (31) F. Wu, T. Lovorn, and A. H. MacDonald,  _Theory of optical absorption by interlayer excitons in transition metal dichalcogenide heterobilayers_, Phys. Rev. B 97, 035306 (2018).
  * (32) F. Wu, T. Lovorn, and A. H. MacDonald,  _Topological Exciton Bands in Moiré Heterojunctions_, Phys. Rev. Lett. 118, 147401 (2017).
  * (33) P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. de Gironcoli, P. Delugas, R. A. J. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Kucukbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni,  _Advanced capabilities for materials modelling with Quantum ESPRESSO_, J. Phys.: Condens. Matter 29, 465901 (2017).
  * (34) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch,  _QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials_, J. Phys.: Condens. Matter 21, 395502 (2009).
  * (35) T. Sohier, M. Calandra, and F. Mauri,  _Density functional perturbation theory for gated two-dimensional heterostructures: Theoretical developments and application to flexural phonons in graphene_, Phys. Rev. B 96, 075448 (2017).
  * (36) F. D. M. Haldane,  _Model for a Quantum Hall Effect without Landau Levels: Condensed-Matter Realization of the ”Parity Anomaly”_, Phys. Rev. Lett. 61, 2015 (1988).
  * (37) C. L. Kane and E. J. Mele,  _Quantum Spin Hall Effect in Graphene_, Phys. Rev. Lett. 95, 226801 (2005).
  * (38) C. L. Kane and E. J. Mele,  _Z2Topological Order and the Quantum Spin Hall Effect_, Phys. Rev. Lett. 95, 343 (2005).
  * (39) T. Fukui, Y. Hatsugai, and H. Suzuki,  _Chern Numbers in Discretized Brillouin Zone: Efficient Method of Computing (Spin) Hall Conductances_, J. Phys. Soc. Jpn. 74, 1674 (2005).
  * (40) N. Bultinck, S. Chatterjee, and M. P. Zaletel,  _Mechanism for Anomalous Hall Ferromagnetism in Twisted Bilayer Graphene_, Phys. Rev. Lett. 124, 166601 (2020).
  * (41) L. Rademaker, I. V. Protopopov, and D. A. Abanin,  _Topological flat bands and correlated states in twisted monolayer-bilayer graphene_, Phys. Rev. Research 2, 033150 (2020).