# A Point in the Right Direction: Vector Prediction for Spatially-aware Self-
supervised Volumetric Representation Learning

###### Abstract

High annotation costs and limited labels for dense 3D medical imaging tasks
have recently motivated an assortment of 3D self-supervised pretraining
methods that improve transfer learning performance. However, these methods
commonly lack spatial awareness despite its centrality in enabling effective
3D image analysis. More specifically, position, scale, and orientation are not
only informative but also automatically available when generating image crops
for training. Yet, to date, no work has proposed a pretext task that distills
all key spatial features. To fulfill this need, we develop a new self-
supervised method, VectorPOSE, which promotes better spatial understanding
with two novel pretext tasks: Vector Prediction (VP) and Boundary-Focused
Reconstruction (BFR). VP focuses on global spatial concepts (i.e., properties
of 3D patches) while BFR addresses weaknesses of recent reconstruction methods
to learn more effective local representations. We evaluate VectorPOSE on three
3D medical image segmentation tasks, showing that it often outperforms state-
of-the-art methods, especially in limited annotation settings.

Index Terms—  3D Self-Supervised Learning, Representation Learning, Volumetric
Medical Image Segmentation.

##  1 Introduction

Modern 3D medical image analysis techniques have made great performance
strides by extracting appearance-based features from local image patches.
However, they fail to explicitly attend to a key capability that facilitates
volumetric analysis: _spatial awareness_. Given an arbitrary image crop (e.g.,
see Fig. 1), we can often leverage anatomical priors and infer a large amount
of information from its position, scale, and orientation. Crop position such
as its axial window (e.g., the dashed red and blue lines in Fig. 1 that span
the CT volume depth) reveals which anatomical structures to expect. Image zoom
is commonly adjusted by practitioners, and the knowledge of a crop’s scale is
essential for accurate assessments of the structure size, shape, and extent.
Finally, crop orientation informs object poses and valid adjacent structures.
From these insights, we hypothesize that _networks which better predict 3D
crop position, scale, and orientation can more accurately segment images with
superior anatomical knowledge_.

Fig. 1: Spatial properties for effective 3D crop analysis. On the right, one
may see the insights brought by crop position, scale, and orientation using
spatial anatomical priors. On the left, one may see the importance of boundary
delineation.

Fig. 2: An overview of our proposed VectorPOSE approach.

To distill relevant spatial information without additional annotation burden,
we leverage self-supervised learning (SSL). Recent 3D SSL methods have
improved downstream task performances. They generally fall into three types.
1) Contrastive learning [1, 2, 3, 4, 5, 6] stems from metric learning where
representations of similar images (e.g., two augmented versions of an image)
are encouraged to be close while features of distinct images are contrasted.
However, many 3D medical tasks have prohibitively small datasets while
requiring cubically more memory, which is incongruous with contrastive
methods’ reliance on large datasets, batch sizes, and models [7]. 2)
Reconstruction [8, 9, 10, 2, 11] has been shown to outperform supervised
transfer learning on 3D data by restoring an original image from a noised or
masked version. Still, when analyzing reconstruction results in [9] (see the
left side of Fig. 1), one can see significant losses of boundary acuity as
well as tissue details. This may suggest an over-focus on local intensity
statistics over more robust semantic or spatial understandings of the whole
structures. 3) Prediction-based pretext tasks entail solving spatial “puzzles”
such as relative position prediction [1], jigsaw unscrambling [1, 12, 13], or
scale classification [10]. However, these methods are prone to short-cut
learning and over-emphasize discriminative elements (e.g., patch boundaries).
Moreover, some approaches propose to learn spatial properties incompletely
[10] or indirectly through disparate pretext tasks [9, 14], but implicit
learning is suboptimal and uncontrolled. To our best knowledge, no approach
adequately covers all key spatial properties important to 3D medical image
analysis in practice.

To remedy this gap, we propose a new 3D self-supervised pretraining approach,
VectorPOSE (POSE for Position, Orientation, Scale, and Extent), with two novel
pretext tasks for learning improved spatial and appearance features. Our first
pretext task, Vector Prediction (VP), infers vectors that originate from
predefined points within a sampled input crop (e.g., crop center, corners) and
terminate at the center of the volume. Vectors encapsulate both magnitude and
direction, so predicting multiple vectors with consistent origin points
effectively encodes position (with respect to the terminal point), orientation
(via vector angles to the terminal point), and scale (where differentials
between predicted vector origins represent crop extent). To learn effective
local features and address boundary degredations in reconstruction, our second
pretext task, Boundary-Focused Reconstruction (BFR), conducts: 1) boundary
reconstruction where extracted edges of the patch are explicitly predicted,
and 2) voxel reconstruction where original intensity values are predicted with
a loss that better preserves boundaries. To demonstrate efficacy, we pretrain
and evaluate on three 3D semantic medical image segmentation tasks (i.e., CT
cardiac, CT abdominal, MR prostate). VectorPOSE generally outperforms state-
of-the-art single- and multi-task methods while exhibiting superior data
efficiency compared to known contrastive learning methods.

Our contributions are summarized below.

  1. 1.

We propose a self-supervised approach, VectorPOSE, that places a holistic
emphasis on three key spatial properties (position, orientation, scale) while
synergizing learned global information with local, boundary-aware cues.

  2. 2.

As a cornerstone to spatial understanding, we introduce a novel Vector
Prediction (VP) pretext task that succinctly encapsulates key spatial
properties and is both effective and extensible. Further, we propose the
Boundary-Focused Reconstruction (BFR) pretext task to address pitfalls of
recent works and learn improved local features.

  3. 3.

VectorPOSE is comprehensively evaluated on two 3D CT segmentation datasets
(MMWHS [15, 16] & BCV [17]) and one 3D MR prostate dataset [18]. Not only does
it outperform recent self-supervised methods in the majority of settings, it
also exhibits improved data efficiency without additional training time.

##  2 Methodology

Before detailing the proposed pretext tasks, we first overview our general
pipeline and its major components. Following established self-supervised
methods, we pretrain a randomly initialized neural network consisting of an
encoder $E$, a decoder $D$, the prediction head $D_{\textit{vp}}$ for Vector
Prediction (VP), and the prediction head $D_{\textit{bfr}}$ for Boundary-
Focused Reconstruction (BFR). Note that $D_{\textit{vp}}$ employs global
information and processes pooled & flattened output features from $E$, while
$D_{\textit{bfr}}$ utilizes $D$’s output that has the same resolution as $E$’s
input. To pretrain $\\{E,D,D_{\textit{vp}},D_{\textit{bfr}}\\}$ with
VectorPOSE (see Fig. 2), we input patches $x^{i(j)},x^{i(k)},\ldots$ from the
$i$th image $X^{i}$ of an unlabeled dataset $\mathcal{D}_{u}$. Stochastic
spatial augmentations $\mathcal{T}_{S}$ (e.g., flipping or rotating along each
axis) and intensity noising $\mathcal{T}_{I}$ (e.g., local pixel shuffle,
masking, intensity shifting) are applied to all input patches independently
and in order. The targets for VP & BFR are computed after $\mathcal{T}_{S}$ is
applied. After pretraining, $D_{\textit{vp}},D_{\textit{bfr}}$ are discarded
while $E,D$ are transferred and fine-tuned on a downstream target task.

###  2.1 Vector Prediction (VP)

The objective of VP is to distill position, scale, and orientation information
from crops by predicting multiple vectors which all point to one standardized
location in all dataset volumes (e.g., volume center or center of a chosen
object). We refer to this point as the reference landmark,
$(x^{i(j)}_{\dashv},y^{i(j)}_{\dashv},z^{i(j)}_{\dashv})$, which denotes the
landmark point for crop $x^{i(j)}$ from image $X^{i}$. Reference landmarks
across images have anatomical correspondences from which we learn useful
priors and structural consistencies. In many 3D tasks, reference landmarks can
simply be taken as the center of volumes given rough image alignments. In
others, positions of landmarks can be extracted in an unsupervised manner
(e.g., center of the left lung as a landmark after lung masks are
automatically extracted [19]). To regularize the reference landmark, we
randomly jitter coordinates by $\eta$% of each image dimension’s total length.

After image landmarks are obtained, a vector’s origin point concretely defines
distance (via its magnitude) and position (via angles) to the reference. For a
sampled crop $x^{i(j)}$, we define $n$ vector origin points as
$(x^{i(j)}_{0\vdash},y^{i(j)}_{0\vdash},z^{i(j)}_{0\vdash})$, $\ldots$,
$(x^{i(j)}_{n-1\vdash},y^{i(j)}_{n-1\vdash},z^{i(j)}_{n-1\vdash})$. In
practice, we set these origin points in consistent positions for all crops
(e.g., $m$=$0$ at patch center, $m$=$1$ at patch top-right corner, etc.).
Mathematically, predicting a single vector is sufficient to describe position
and orientation, while two are enough to also cover scale. However, we
empirically find that $n$=$9$ vectors (vectors originating from the patch
center and 8 corners) improve feature learning (see Section § 3, Tab. 2).

Thus, all vectors belonging to crop $x^{i(j)}$ are defined as $v_{m}^{i(j)}$
($m=0,\ldots,n-1$) with origin point
$(x^{i(j)}_{m\vdash},y^{i(j)}_{m\vdash},z^{i(j)}_{m\vdash})$ and terminal
point $(x^{i(j)}_{\dashv},y^{i(j)}_{\dashv},z^{i(j)}_{\dashv})$. The terminal
points are identical for all crops from the same image. Below we omit “$i(j)$”
for readability. For prediction, we reparameterize vectors into spherical
coordinates: let
$(x_{m},y_{m},z_{m})=(x_{\dashv}-x_{m\vdash},y_{\dashv}-y_{m\vdash},z_{\dashv}-z_{m\vdash})$;
so, $v_{m}=$ ($r_{m}$, $\theta_{m}$, $\rho_{m}$), where $r_{m}=$
$\sqrt{x_{m}^{2}+y_{m}^{2}+z_{m}^{2}}$, $\theta_{m}=\arccos(z_{m}/r_{m})$, and
$\phi_{m}=\arctan(y_{m}/x_{m})$ with constraints $r\in[0,\inf)$,
$\theta\in[0,\pi]$, and $\phi\in[-\pi,\pi]$. The model regresses the $n$
normalized spherical coordinates of a crop via the VP loss:

$\displaystyle\mathcal{L}_{\textit{vp}}=$
$\displaystyle\frac{1}{n}\sum_{m=0}^{n-1}\left(||\frac{r_{m}}{R}-\sigma(\hat{r}_{m})||_{1}+||\frac{\theta_{m}}{\pi}-\sigma(\hat{\theta}_{m})||_{1}+\right.$
(1) $\displaystyle\min($
$\displaystyle||\frac{\phi_{m}}{\pi}-\tanh(\hat{\phi}_{m})||_{1},||\frac{\phi_{m}+2\pi}{\pi}-\tanh(\hat{\phi}_{m})||_{1},$
$\displaystyle||\left.\frac{\phi_{m}-2\pi}{\pi}-\tanh(\hat{\phi}_{m})||_{1})\right)$

where $\sigma$ is the sigmoid function, $R$ is the radius of a sphere
circumscribing the image volume, and $\hat{r}_{m}$, $\hat{\theta}_{m}$,
$\hat{\phi}_{m}$ are logits from $E\circ D\circ D_{\textit{vp}}$ (see the
bottom right of Fig. 2). To address the $\phi$ angles that are physically
close but distant in the parameter space (e.g., $-179\degree$ &
+179$\degree$), we take the minimum among targets
$\phi_{m},\phi_{m}-2\pi,\phi_{m}+2\pi$ (see Eq. (1)). Also, note that when
crops undergo spatial augmentations (e.g., flipping, rotations), the index of
each vector’s post-transform origin corresponds with the origin before
transforming.

By encapsulating the three key spatial properties into an interdependent task,
we prevent training instability and feature interference from independent
predictions like when classifications were used [10]. Additionally, vectors
represent spatial properties in a continuous space which is a generalized
extension of discrete rotation or scale prediction [13].

###  2.2 Boundary-Focused Reconstruction (BFR)

Complementing the global features from VP, two appearance-focused losses are
employed in BFR. To address previous issues of boundary degradation and
reduced texture acuity, we first switch the reconstruction criterion from
$L_{2}$ loss (as used in [9]) to $L_{1}$ loss, which has been shown to improve
both reconstuction quality and boundary accuracy [11]. To further emphasize
anatomical delineation, we add an additional boundary reconstruction task
where boundaries are extracted via a 3D Scharr edge detector. We select the
Scharr transform since it exhibits superior rotational invariance compared to
other filters like Sobel. In conjunction with strong augmentations such as
texture noising (e.g., local pixel shuffle) and masking (e.g., in-painting or
out-painting), explicit edge reconstruction facilitates the learning of
general anatomical shapes and boundary extents. The overall BFR loss on a crop
is defined as:

$\mathcal{L}_{\textit{bfr}}=\frac{1}{n}\sum_{m=0}^{n-1}(||{y}_{m}^{v}-\sigma(\hat{y}^{v}_{m})||_{1}+\alpha||y_{m}^{b}-\sigma(\hat{y}^{b}_{m})||_{1}),\vspace{-2.2mm}$
(2)

where ${y}_{m}^{v}$ and $y_{m}^{b}$ are the targets of the voxel and boundary
reconstruction tasks, $\hat{y}^{v}_{m}$ and $\hat{y}^{b}_{m}$ (from $E\circ
D\circ D_{\textit{bfr}}$) are the voxel and boundary logits, respectively.
$\alpha$ is a scaling term to address the small proportion of boundaries in
crops.

###  2.3 Overall Loss and Implementation

The overall VectorPOSE loss for an image crop is:

$\mathcal{L}=\lambda\mathcal{L}_{\textit{bfr}}+(1-\lambda)\mathcal{L}_{\textit{vp}}.\vspace{-2mm}$
(3)

For our UNet-like model, we employ a 3D version of ResNet-50 [20] as the
encoder and attach it to a light decoder with bilinear upsampling layers &
additions for feature fusion. For $D_{\textit{vp}}$, the lowest encoder output
is pooled, flattened, and processed through a 2-layer MLP with 256 hidden
dimensions, while $D_{\textit{bfr}}$ uses a single 1x1 convolution. We use
$m$=$9$ vectors per patch with a jitter $\eta$=$5$%, and train with
$\alpha$=$5$ & $\lambda$=$0.5$.

Method | CT MMWHS |  | CT BCV |  | MRI MSD-Prostate
---|---|---|---|---|---
10% | 25% | 50% | 100% |  | 10% | 25% | 50% | 100% |  | 10% | 25% | 50% | 100%
Random Init. | 79.31 | 87.29 | 89.15 | 90.70 |  | 30.06 | 58.48 | 66.23 | 75.81 |  | 40.05 | 58.78 | 69.25 | 74.24
Models Genesis [9] | 80.05 | 88.14 | 90.54 | 91.05 |  | 32.03 | 58.80 | 66.87 | 75.79 |  | 44.00 | 60.89 | 69.55 | 75.86
SAR [10] | 81.48 | 88.57 | 90.81 | 91.14 |  | 35.63 | 60.29 | 67.52 | 76.77 |  | 42.53 | 61.02 | 70.16 | 74.62
Rubik++ [11] | 81.08 | 88.40 | 90.94 | 91.23 |  | 33.28 | 59.83 | 67.72 | 76.65 |  | 43.70 | 61.97 | 71.35 | 75.81
PGL [21] | 79.76 | 87.50 | 89.48 | 90.57 |  | 30.99 | 58.56 | 66.39 | 75.98 |  | 41.36 | 59.08 | 69.97 | 74.08
MoCo [4] | 80.22 | 87.94 | 89.60 | 91.05 |  | 31.89 | 58.95 | 66.91 | 76.01 |  | 42.07 | 60.69 | 70.17 | 74.31
VectorPOSE (Ours) | 83.30 | 89.43 | 91.27 | 91.60 |  | 37.88 | 61.61 | 68.01 | 76.53 |  | 46.35 | 63.14 | 71.22 | 75.97
p-value | .0335 | .0286 | .0437 | .0564 |  | .0179 | .0156 | .0388 | .0593 |  | .0082 | .0098 | .0571 | .0681
Table 1: Class-averaged Dice score (%) comparisons vs. recent self-supervised 3D medical image segmentation methods. Entries that are bolded & underlined are the best & second-best scores, respectively. p-values are obtained from independent t-tests of the best and second-best scores. Voxel Rec. | $\mathcal{L}_{1}$ | $\mathcal{L}_{1}$ | $\mathcal{L}_{1}$ | $\mathcal{L}_{1}$ | $\mathcal{L}_{1}$ | $\mathcal{L}_{1}$
---|---|---|---|---|---|---
Bound. Rec. |  | $\mathcal{L}_{1}$ | $\mathcal{L}_{1}$ | $\mathcal{L}_{1}$ | $\mathcal{L}_{1}$ | $\mathcal{L}_{1}$
Center Vector |  |  | ✓ | ✓ | ✓ | ✓
Corner Vector(s) |  |  |  | 1 | 4 | 8
Dice (%) | 80.16 | 80.80 | 81.37 | 82.29 | 82.75 | 83.30
Table 2: Ablations on the proposed components on 10% MMWHS. The top two rows
contain appearance-based components while the middle two rows have VP
settings.

##  3 Experiments & Results

We implemented all experiments in PyTorch and trained on NVIDIA Tesla P100s
(16GB VRAM). For pretraining, volume intensities were normalized between 0 and
1, uninformative regions were excluded, (96, 96, 96) crops were sampled, and
augmentations in [9] were applied. We trained using AdamW with a $0.0002$
learning rate & $0.0001$ weight decay for 300 epochs ($\sim$21 hours) with a
batch size of 12. For fine-tuning, we normalized volume intensities between 0
and 1, sampled (64, 128, 128) crops, and applied random flipping, brightness,
gamma, and blurring. We optimized using AdamW with a learning rate $0.001$ and
a weight decay $0.001$ for 200 epochs ($\sim$14 hours) with a batch size of 4
over 8 runs. For evaluation, we used the average Dice-Sørensen Coefficient
across foreground classes.

###  3.1 Datasets

All the datasets are split using a 6:2:2 training:validation:test ratio.
Models are pretrained using the full training sets and fine-tuned using 10%,
25%, 50%, or 100% of the training data (if the number of training samples is
not divisible, the floor function is applied).

CT MMWHS [16, 15] provides 20 annotated CT volumes segmenting seven cardiac
structures (left ventricle, left atrium, right ventricle, right atrium,
myocardium, ascending aorta, and pulmonary artery). CT BCV [17] is an
abdominal organ segmentation challenge with 30 annotated CT volumes and 13
classes. MRI MSD-Prostate [18] contains 32 labeled MRI (T2, ADC) volumes with
two foreground classes.

###  3.2 Results and Discussion

For fair comparisons, all the methods were pretrained and fine-tuned on the
same data splits with the same number of iterations, and allowed an equal
budget for hyperparameter tuning. The performances on the three downstream
datasets are summarized in Table 1.

From our experimental results, we make the following observations. 1) Our
approach generally outperforms recent state-of-the-art multi-task methods and
handily beats training from scratch, which supports our hypothesis that more
spatially aware models yield more effective representations of 3D crops. More
details on the effects of spatially-aware components are discussed in Section
§ 3.3. 2) The fewer the annotations, the bigger the performance improvements.
This is a desirable property for tasks with many biomedical objects with
extremely limited labels. Our outperformance over other methods that also use
reconstruction [9, 10, 11] shows both the efficacy of our spatial pretext task
and superior data efficiency. 3) Both our encoder and decoder are initialized
with multi-scale and semantically diverse tasks. Many methods only learn
encoder features while suboptimally neglecting decoder parameters (e.g.,
[21]). Additionally, the decoder receives heterogeneous learning signals from
both local-focused BFR and global, spatially-aware VP. This allows for each
segment of the model to learn relevant and transferable features and serves as
a natural regularizer that prevents overfitting to any single task. In fact,
we notice that joint training leads to synergistic improvements. 4) Lesser
performance from PGL & MoCo may indicate innate disadvantages in using metric
learning techniques with 3D radiology tasks. We hypothesize that this is in
part due to the fact that PGL only uses positive samples which inhibits
feature learning. The small amount of pretraining images also limits
effectiveness. Given that many 3D medical images are highly structured and
prior-dependent, we can more effectively learn these priors through spatial
prediction tasks that holistically cover multiple key properties.

###  3.3 Ablations

Here, we explore the contributions of each individual component in our method.
Notably, we intend to test our hypothesis about the importance of explicitly
learning spatial parameters and compare its efficacy against other proposed
methods (e.g., reconstruction). In Table 2, one can see $\sim$0.6% improvement
after adding the boundary awareness task, which supports the benefit of
explicitly predicting boundaries over voxel reconstruction alone. Next, we see
that with increasing $m$ (the number of vectors per crop), performance notably
rises. This supports our hypothesis of using spatial awareness as an effective
feature-learning principle.

##  4 Conclusions

We proposed VectorPOSE, a new self-supervised approach for 3D medical image
segmentation with two novel pretext tasks: 1) Vector Prediction (VP) to
address the lack of spatial understanding (i.e., position, scale, orientation)
in existing works; 2) Boundary-Focused Reconstruction (BFR) for improving the
understanding of anatomical structure extents through boundary regression.
After evaluation on three 3D semantic segmentation datasets, VectorPOSE’s
general outperformance over state-of-the-art methods shows that spatial
understanding is pertinent for feature learning and may be a promising avenue
for future research progress.

##  5 Compliance with ethical standards

This research was conducted retrospectively using open access human subject
data by three publicly available datasets [16, 17, 18]. Ethical approval was
not required as confirmed by the licenses attached with the open access
datasets.

## References

  * [1] A. Taleb, W. Loetzsch, N. Danz, J. Severin, T. Gaertner, B. Bergner, and C. Lippert,  “3D self-supervised methods for medical imaging,”  NIPS, vol. 33, pp. 18158–18172, 2020.
  * [2] H. Zheng, J. Han, H. Wang, L. Yang, Z. Zhao, C. Wang, and D. Chen,  “Hierarchical self-supervised learning for medical image segmentation based on multi-domain data aggregation,”  in MICCAI, 2021.
  * [3] D. Zeng, Y. Wu, X. Hu, X. Xu, H. Yuan, M. Huang, J. Zhuang, J. Hu, and Y. Shi,  “Positional contrastive learning for volumetric medical image segmentation,”  in MICCAI, 2021.
  * [4] X. Chen, L. Yao, T. Zhou, J. Dong, and Y. Zhang,  “Momentum contrastive learning for few-shot COVID-19 diagnosis from chest CT images,”  Pattern Recognition, vol. 113, pp. 107826, 2021.
  * [5] K. Chaitanya, E. Erdil, N. Karani, and E. Konukoglu,  “Contrastive learning of global and local features for medical image segmentation with limited annotations,”  NIPS, vol. 33, 2020.
  * [6] H. Zhou, C. Lu, S. Yang, X. Han, and Y. Yu,  “Preservational learning improves self-supervised medical image models by reconstructing diverse contexts,”  ICCV, pp. 3479–3489, 2021.
  * [7] T. Chen, S. Kornblith, M. Norouzi, and G. Hinton,  “A simple framework for contrastive learning of visual representations,”  in ICML. PMLR, 2020, pp. 1597–1607.
  * [8] F. Matzkin, V. F. J. Newcombe, S. Stevenson, A. Khetani, T. Newman, R. Digby, A. Stevens, B. Glocker, and E. Ferrante,  “Self-supervised skull reconstruction in brain CT images with decompressive craniectomy,”  ArXiv, vol. abs/2007.03817, 2020.
  * [9] Z. Zhou, V. Sodha, M. Siddiquee, R. Feng, N. Tajbakhsh, M. Gotway, and J. Liang,  “Models genesis: Generic autodidactic models for 3D medical image analysis,”  MICCAI, vol. 11767, pp. 384–393, 2019.
  * [10] X. Zhang, S. Feng, Y. Zhou, Y. Zhang, and Y. Wang,  “SAR: Scale-aware restoration learning for 3D tumor segmentation,”  in MICCAI, 2021.
  * [11] X. Tao, Y. Li, W. Zhou, K. Ma, and Y. Zheng,  “Revisiting Rubik’s cube: Self-supervised learning with volume-wise transformation for 3D medical image segmentation,”  in MICCAI. Springer, 2020, pp. 238–248.
  * [12] X. Zhuang, Y. Li, Y. Hu, K. Ma, Y. Yang, and Y. Zheng,  “Self-supervised feature learning for 3D medical images by playing a Rubik’s cube,”  in MICCAI. Springer, 2019, pp. 420–428.
  * [13] J. Zhu, Y. Li, Y. Hu, K. Ma, S. Zhou, and Y. Zheng,  “Rubik’s Cube+: A self-supervised feature learning framework for 3D medical image analysis,”  Medical Image Analysis, vol. 64, pp. 101746, 2020.
  * [14] F. Haghighi, M. Taher, Z. Zhou, M. Gotway, and J. Liang,  “Transferable visual words: Exploiting the semantics of anatomical patterns for self-supervised learning,”  IEEE TMI, vol. 40, pp. 2857–2868, 2021.
  * [15] X. Zhuang and J. Shen,  “Multi-scale patch and multi-modality atlases for whole heart segmentation of MRI,”  Medical Image Analysis, vol. 31, pp. 77–87, 2016.
  * [16] X. Zhuang,  “Challenges and methodologies of fully automatic whole heart segmentation: A review,”  Journal of Healthcare Engineering, vol. 4, pp. 371–407, 2013.
  * [17] “Multi-atlas labeling beyond the cranial vault - workshop and challenge,” 2015,  Accessed January 2021 at https://www.synapse.org/#!Synapse:syn3193805/wiki/89480.
  * [18] A. L. Simpson, M. Antonelli, S. Bakas, M. Bilello, K. Farahani, B. Ginneken, A. Kopp-Schneider, B. A. Landman, G. J. S. Litjens, B. H. Menze, O. Ronneberger, R. M. Summers, P. Bilic, P. F. Christ, R. K. G. Do, M. J. Gollub, J. Golia-Pernicka, S. Heckers, W. R. Jarnagin, M. McHugo, S. Napel, E. Vorontsov, L. Maier-Hein, and M. J. Cardoso,  “A large annotated medical image dataset for the development and evaluation of segmentation algorithms,”  ArXiv, vol. abs/1902.09063, 2019.
  * [19] B. Zhao, G. Gamsu, M. S. Ginsberg, L. Jiang, and L. H. Schwartz,  “Automatic detection of small lung nodules on CT utilizing a local density maximum algorithm,”  Journal of Applied Clinical Medical Physics, vol. 4, no. 3, pp. 248–260, 2003.
  * [20] K. He, X. Zhang, S. Ren, and J. Sun,  “Deep residual learning for image recognition,”  CVPR, pp. 770–778, 2016.
  * [21] Y. Xie, J. Zhang, Z. Liao, Y. Xia, and C. Shen,  “PGL: Prior-guided local self-supervised learning for 3D medical image segmentation,”  ArXiv abs/2011.12640, 2020.


# Proximity-driven ferromagnetism and superconductivity in the triangular
Rashba-Hubbard model

Mehdi Biderang<EMAIL_ADDRESS>Department of Physics and
Astronomy, University of Manitoba, Winnipeg R3T 2N2, Canada  Manitoba Quantum
Institute, University of Manitoba, Winnipeg R3T 2N2, Canada  Mohammad-Hossein
Zare<EMAIL_ADDRESS>Department of Physics, Qom University of Technology, Qom
37181-46645, Iran  Jesko Sirker<EMAIL_ADDRESS>Department of
Physics and Astronomy, University of Manitoba, Winnipeg R3T 2N2, Canada
Manitoba Quantum Institute, University of Manitoba, Winnipeg R3T 2N2, Canada

###### Abstract

Bilayer Moiré structures are a highly tunable laboratory to investigate the
physics of strongly correlated electron systems. Moiré transition metal
dichalcogenides at low-energies, in particular, are believed to be described
by a single narrow band Hubbard model on a triangular lattice with spin-orbit
coupling. Motivated by recent experimental evidence for superconductivity in
twisted bilayer materials, we investigate the possible superconducting
pairings in a two-dimensional single band Rashba-Hubbard model. Using a
random-phase approximation in the presence of nearest and next-nearest
neighbor hopping, we analyze the structure of spin fluctuations and the
symmetry of the superconducting gap function. We show that Rashba spin-orbit
coupling favors ferromagnetic fluctuations which strengthen triplet
superconductivity. If parity is violated due to the absence of spatial
inversion symmetry, singlet ($d$-wave) and triplet ($p$-wave) channels of
superconductivity will be mixed. Moreover, we show that time-reversal symmetry
can be spontaneously broken leading to a chiral superconducting state.
Finally, we consider quasiparticle interference as a possible experimental
technique to observe the superconducting gap symmetry.

##  I Introduction

There has been a growing interest in exotic quantum phasesof twisted van der
Waals bilayer materials in recent years [1, 2, 3, 4, 5, 6, 7, 8, 9]. In such
systems, the band structure is highly tunable by twisting one layer with
respect to the other [10, 11, 12]. In particular, the bandwidth can be tuned
to be smaller than the magnitude of the electron interactions and, as a
result, a strongly correlated electron system is formed [13, 14, 15, 16, 17].
By electron or hole-doping the Mott insulating state, unconventional
superconducting phases can potentially be realized [18, 19, 20, 21]. There
have been many recent theoretical and experimental reports on doped twisted
bilayer graphene (TBG) that are consistent with the formation of a chiral
$d+{\rm i}d$-wave topological superconductivity (TSC) near half-filling [22,
23, 24, 25]. A related Moiré system, which has also been studied both
theoretically and experimentally in recent years, are bilayers of group-VI
semiconducting transition metal dichalcogenides (TMDs) such as WSe2 [26, 27,
27], which realize an effective triangular lattice [18, 10, 27, 28, 29, 30].
These materials can host exotic quantum phases such as correlated insulating
states, unconventional superconductors, fractional quantum Hall states, and
quantum spin liquids [18, 10, 31, 32, 27, 33, 34, 35, 30, 36]. While graphene
behaves like a Dirac semimetal with SU(2) spin rotational symmetry, TMDs are
semiconductors with a large band gap and a large spin-orbit coupling (SOC)
[26, 10, 37]. Under the assumption that disorder can be ignored, TMDs
therefore appear to be easier subjects for theoretical studies than TBG
because of the smaller number of low-energy degrees of freedom and the
presence of narrow bands for a wide range of twisting angles [26, 10, 37]. In
particular, it appears to be possible to understand many of their physical
properties based on an appropriate single-band Hubbard model.

One of the main open questions in condensed matter physics is the pairing
mechanism and the symmetry of the superconducting phase in unconventional
high-temperature superconductors such as the cuprates or the pnictides [38,
39, 40, 41]. There appears to be a consensus that the electron-phonon coupling
is not strong enough to overcome the Coulomb interaction between the electrons
in order to form Cooper pairs. A common denominator of these systems is that
superconducting phases appear in the vicinity of magnetic instabilities and
magnetic ordered phases. A possible candidate for a pairing mechanism in
unconventional superconductors is therefore the exchange of magnons or
paramagnons to form Cooper pairs [42]. Another interesting aspect of non-BCS
type superconductors is that they are capable of hosting TSC such as chiral
$p$-wave and chiral $d$-wave, or helical/chiral parity-mixed states [43, 44,
45]. Specifically, the twisted bilayer of TMDs with an effective triangular or
honeycomb lattice can potentially produce chiral superconductivity if time-
reversal symmetry (TRS) becomes spontaneously broken in their superconducting
ground state.

Figure 1:  (a) Geometry of the considered system: The blue spheres represent
${\cal R}^{M}_{M}$ high-symmetry positions of a twisted bilayer of TMD,
forming an effective triangular lattice on top of a substrate. $t_{1}$ and
$t_{2}$ are the nearest- and the next-nearest- hopping, respectively. Broken
inversion symmetry induces and effective electric field. Also shown is the
band structure of the non-interacting model for $t_{2}=0$, for (b) $V_{\rm
SOC}=0$ and (d) $V_{\rm SOC}=0.3$ at $T=0.03$. (c,e) The corresponding
densities of state. The green line indicates half filling and the red and blue
dashed lines the positions of the van Hove singularities. (f,g) Fermi surface
for $\delta=0$ and $V_{\rm SOC}=0$ and $V_{\rm SOC}=0.3$. The central hexagons
(black dashed lines) indicate the Brillouin Zone. (h,i) In the hole-doped
region, the Fermi surface consists of two hole pockets centered around the
$\bm{\Gamma}$ point. (j,k) Close to the van Hove singularities, the Fermi
surface consists of hole- and an electron-like pockets. (l) For large electron
doping, both helical Fermi surfaces become electron pockets.

Parallel with the discovery of topological superconductors, progress in the
fabrication of heterostructures has been made based on a wide variety of
crystal growth methods, such as molecular beam epitaxy (MBE) [46, 47, 48, 49,
50]. A very important aspect of these structures is that inversion symmetry
is, in general, broken leading to the emergence of antisymmetric Rashba or
Ising SOCs [51, 52, 53, 54, 53, 53, 54, 55, 56]. As a consequence, spin-flip
scattering processes are present and can strongly affect the magnetic
structure and the symmetry of the superconducting state [57, 58]. In
particular, spin is no longer a good quantum number and singlet and triplet
pairings in the superconducting state will be mixed [59, 60, 61].

The locations of the high-symmetry points in a TMD Moiré superlattice are
denoted by ${\cal R}^{\alpha}_{\beta}$, in which $\alpha$ and $\beta$ are
running over metal (M) and chalcogene (X) atoms of the top and bottom layers,
respectively. In general, for a twisted bilayer of TMD these high-symmetry
points are: ${\cal R}^{M}_{M}$, ${\cal R}^{X}_{M}$, and ${\cal R}^{M}_{X}$
[10, 37]. The structure has a D3 point-group symmetry for rotations around the
$\hat{z}$ axis along with a twofold C2v symmetry around the in-plane $\hat{y}$
axis [10, 37]. The effective Moiré lattice of the twisted bilayer of TMDs can
be easily changed by tuning the local density of states (LDOS) at high-
symmetry positions using a gate voltage between the layers [19, 13, 37]. There
are two possible structures for the effective lattice of a twisted bilayer of
TMDs: (i) honeycomb, or (ii) triangular lattices.

In this paper, we will focus on the effective triangular lattice. This case is
realized if the ${\cal R}^{M}_{M}$ points have a much larger LDOS than the
other high-symmetry points. We aim to study the magnetic and superconducting
instabilities in a twisted bilayer of TMD which has been deposited on top of a
bulk substrate. This setup, sketched in Fig. 1(a), will combine the physics of
strongly correlated electrons with Rashba SOC.

We note that in real twisted TMD heterostructures there is also a substrate
top layer and a top gate. Here we show schematically the structure without the
top layer. In the limit of a single narrow band being relevant for the low-
energy physics, we are thus trying to understand the different phases and
potential pairings of a Rashba-Hubbard model on an effective triangular
lattice. For large twist angles, which can be realized in TMDs, the Moiré unit
cell becomes relatively small and hopping integrals between next-nearest
neighbors might no longer be negligible. We therefore will include such
processes as well. We will consider both the magnetic susceptibility and a
magnon-mediated pairing using the random phase approximation (RPA). To find
the symmetry of the dominant superconducting pairing channel, we will then
solve the self-consistent gap equation within BCS theory.

Our paper is organized as follows: In Sec. II we will define the considered
model and calculate its non-interacting band structure and DOS for various
doping levels. In Sec. III we consider the bare and RPA spin susceptibilities
as a function of hopping amplitudes, SOC, and doping level. Also within the
RPA approximation, we then consider in Sec. IV the pairing of electrons,
separating the effective vertex into singlet and triplet channels and
neglecting its frequency dependence. In Sec. V we then discuss the response
which would be seen in a quasi-particle interference experiment. We note that
this scanning tunneling microscopy (STM) technique is only possible on a
heterostructure without a substrate top layer and without a top gate as shown
schematically in Fig. 1(a) and is thus, for now, a purely theoretical
consideration. In the final section we then provide a brief summary and
discussion of our main results.

##  II Model and non-interacting band structure

We start by introducing the two-dimensional Hubbard model which we will
consider in the following and discuss the band structure and DOS in the non-
interacting limit.

###  II.1 Rashba-Hubbard Model on the Triangular Lattice

Motivated by recent experiments on twisted heterostructures which have an
effective triangular Moiré lattice [25, 19, 13, 31, 29, 30, 62, 63], we
consider a two-dimensional lattice of atoms deposited on a thick substrate,
see Fig. 1(a). At small twist angles, which cause a Moiré pattern with an
enlarged unit cell [26], there are hundreds of energy bands due to the high
number of electrons inside the unit cell [22]. However, nearly flat low-energy
bands with a total bandwidth of only several meV [22] can arise which are well
isolated from the other high energy bands. To understand the low-energy
physics in this situation, it thus suffices to consider only the low-energy
bands. If these bands are partially filled, magnetic and superconducting
instabilities can occur [20, 16, 19]. In addition, the broken inversion
symmetry at the intersection of two subsystems generates a perpendicular
electric field, which induces the well-known Rashba SOC [51]. The single-
orbital Rashba-Hubbard model on a two-dimensional triangular lattice is then
expressed as

${\cal H}={\cal H}_{0}+{\cal H}_{\rm int},$ (1)

in which the first term denotes the non-interacting part of the Hamiltonian,
including kinetic and SOC terms,

${\cal H}_{0}={\cal H}_{\rm kin}+{\cal H}_{\rm SOC}.$ (2)

In real space, the kinetic part is given by

${\cal H}_{\rm kin}=-\mu-t_{1}\\!\\!\sum_{\langle
ij\rangle,\sigma}c^{\dagger}_{i\sigma}c_{j\sigma}-t_{2}\\!\\!\sum_{\langle\\!\langle
ij\rangle\\!\rangle,\sigma}c^{\dagger}_{i\sigma}c_{j\sigma}+{\rm h.c.},$ (3)

with $\mu$, $t_{1}$, and $t_{2}$ being the chemical potential, and the nearest
and next-nearest neighbor hopping amplitudes, respectively. The operator
$c^{\dagger}_{i\sigma}$ ($c_{i\sigma}$) creates (annihilates) an electron with
spin $\sigma=\uparrow,\downarrow$ at lattice site $i$. Hereafter, we set
$t_{1}=1$, unless specifically mentioned otherwise. As already eluded to
earlier, we include a hopping $t_{2}$ because for large twist angles such
hopping processes might no longer be negligible and they can be tuned using an
out-of-plane displacement field [37]. However, $t_{2}$ is in general expected
to be quite small. In momentum space, the kinetic part of the Hamiltonian is
given by ${\cal H}_{\rm kin}\\!=\\!\sum_{{\bf k},\sigma}\varepsilon_{{\bf
k}}c^{\dagger}_{{\bf k}\sigma}c_{{\bf k}\sigma}$ with dispersion

$\varepsilon_{{\bf k}}=-\mu-2t_{1}\sum_{\gamma_{1}=1}^{6}\cos({\bf
k}\cdot\bm{{\bf r}}_{\gamma_{1}})-2t_{2}\sum_{\gamma_{2}=1}^{6}\cos({\bf
k}\cdot\bm{{\bf R}}_{\gamma_{2}}).$ (4)

Here, ${\bf
r}\\!\\!\in\\!\\!\\{\pm\hat{x},\pm\frac{1}{2}\hat{x}\pm\frac{\sqrt{3}}{2}\hat{y}\\}$,
and ${\bf
R}\\!\\!\in\\!\\{\pm\frac{3}{2}\hat{x}\pm\frac{\sqrt{3}}{2}\hat{y},\pm\sqrt{3}\hat{y}\\}$
are the real space vectors connecting the nearest- and next-nearest sites,
respectively. Spin and orbital degrees of freedom are coupled by the SOC
Hamiltonian

${\cal H}_{\rm SOC}=V_{\rm SOC}\sum_{{\bf k},\sigma\sigma^{\prime}}({\bf
g}_{{\bf k}}\cdot\hat{\bm{\sigma}})_{\sigma\sigma^{\prime}}c^{\dagger}_{{\bf
k}\sigma}c_{{\bf k}\sigma},$ (5)

in which the antisymmetric Rashba SOC (ASOC) ${\bf g}$-vector

${\bf g}_{{\bf k}}=\Big{(}\frac{\partial\varepsilon_{{\bf k}}}{\partial
k_{y}},-\frac{\partial\varepsilon_{{\bf k}}}{\partial k_{x}},0\Big{)},$ (6)

is an odd function with respect to momentum and $V_{\rm SOC}$ is the coupling
constant. This type of ASOC plays the role of a momentum-dependent Zeeman
field that splits the two-fold spin degenerate bands into positive and
negative helical bands with energy dispersion

$\xi_{{\bf k},s}=\varepsilon_{{\bf k}}+s|{\bf g}_{{\bf k}}|,$ (7)

and helicity $s=\pm$, which is taken as a band index. We assume, furthermore,
that the Coulomb repulsion can be truncated to an on-site Hubbard interaction

${\cal H}_{\rm int}=U\sum_{i}n_{i\uparrow}n_{i\downarrow}\;.$ (8)

Here $U$ is the Hubbard constant and $n_{i\sigma}$ the number operator at site
$i$ with spin $\sigma$. Throughout this paper, we set $\hbar=k_{B}=1$.

###  II.2 Non-interacting band structure

Let us first consider the non-interacting case. Fig. 1(b) depicts the non-
interacting band structure for the $\Gamma$MK$\Gamma$ high-symmetry path with
$t_{2}=0$ and $V_{\rm SOC}=0$ while Fig. 1(d) shows the band structure with
$t_{2}=0$ and $V_{\rm SOC}=0.3$. The bandwidth in this example is almost
$W\sim 10t_{1}$ but can easily be tuned by changing the twist angle between
the layers. The dashed lines indicate half filling and the positions of the
van Hove singularities (VHSs), respectively. Note that for $V_{\rm SOC}=0$ the
bands are largely flat around the time-reversal invariant K-point leading to a
type-I VHS [64], see Fig. 1(c). In the presence of SOC, the saddle points
appear along the ${\bf K-M}$ and ${\bf K-\Gamma}$ paths, away from time-
reversal invariant momenta, leading to type-II VHSs as shown in Fig. 1(e).
These type-II VHSs are expected to lead to a competition of singlet and
triplet pairing in the superconducting phase. We observe that when comparing
the DOS to the experimental data in Ref. 27, the positions of the van Hove
singularities in our model seem to correspond to a non-zero displacement field
in the experiment. It is not clear to us why that is the case.

Spin-orbit coupling leads to a splitting of the Fermi surface into two helical
surfaces as is shown for the half-filled case in Fig. 1(f,g). For a fixed
strength of SOC we can then study the evolution of the Fermi surface with
doping. For fillings $\langle n\rangle=0.6,~{}0.8$ shown in Fig. 1(h,i), both
helical Fermi surfaces consist of hole-like pockets. Upon approaching the VHSs
the topology of the Fermi surface changes and one of the hole pockets becomes
an electron-like pocket, see Figs. 1(j,k). These pockets are a direct result
of the presence of proximity-induced Rashba SOC and are important for the
superconducting properties of the system. Far away from the VHSs (see Fig.
1(l)), both helical Fermi surface become electron-like pockets.

##  III Spin Susceptibilities

Because we will consider spin fluctuations as the potential glue for electron
pairing and superconductivity, we start by exploring the spin susceptibility
as a function of next-nearest neighbor hopping, SOC, and doping.

###  III.1 Green’s functions and RPA susceptibilities

The free electron Matsubara Green’s function for a system in the presence of
ASOC at $U=0$ is expressed by a $2\times 2$ matrix in spin space,

$\hat{\cal G}^{0}_{{\bf k}}(\rm i\nu_{n})=\Big{[}({\rm
i}\nu_{n}-\varepsilon_{{\bf k}})\hat{\mathbb{I}}-V_{\rm SOC}~{}{{\bf g}}_{{\bf
k}}\cdot\hat{\bm{\sigma}}\Big{]}^{-1}\vspace{0.1 cm}$ (9)

where $\nu_{n}=(2n+1)\pi T$ denote the fermionic Matsubara frequencies and
$\hat{\mathbb{I}}$ is the $2\times 2$ identity matrix. The non-interacting
Green’s function in the spin basis can be mapped into the band space using the
relationship

$\hat{\cal G}^{0}_{{\bf k}}(\rm
i\nu_{n})=\frac{1}{2}\sum_{s=\pm}\Big{[}\hat{\mathbb{I}}+s\hat{{\bf g}}_{{\bf
k}}\cdot\hat{\bm{\sigma}}\Big{]}{\cal G}^{0}_{{\bf k},s}(\rm i\nu_{n}),$ (10)

in which $\hat{{\bf g}}_{{\bf k}}={{\bf g}}_{{\bf k}}/|{{\bf g}}_{{\bf k}}|$,
and ${\cal G}^{0}_{{\bf k},s}(\rm i\nu_{n})=[\rm i\nu_{n}-\xi_{{\bf
k},s}]^{-1}$ is the free electron Matsubara Green’s function in band space. To
investigate the spin fluctuations and electron instabilities of the system, we
adopt the standard RPA approach. In the framework of linear response theory,
the spin susceptibility at $U=0$ with momentum ${\bf q}$ and bosonic Matsubara
frequencies $\omega_{m}=2m\pi T$ is defined by [65]

$\displaystyle\begin{aligned} \hat{\chi}^{(0)\sigma_{3}\sigma_{4}}_{{\bf
q},\sigma_{1}\sigma_{2}}&({\rm i}\omega_{m})=\\\ &-\frac{T}{N}\sum_{{\bf
k},{{\rm i}\nu_{n}}}{\cal G}^{0}_{{\bf k},\sigma_{1}\sigma_{2}}({\rm
i}\nu_{n}){\cal G}^{0}_{{\bf k}+{\bf q},\sigma_{3}\sigma_{4}}({\rm
i}\omega_{m}+{\rm i}\nu_{n}).\end{aligned}$ (11)

Using Eq. (10) and performing the summation over the fermionic Matsubara
frequencies ${\rm i\nu_{n}}$, the bare spin susceptibility in momentum-
frequency space is given by

$\displaystyle\begin{aligned} \hat{\chi}^{(0)\sigma_{3}\sigma_{4}}_{{\bf
q},\sigma_{1}\sigma_{2}}({\rm i}\omega_{m})=\frac{1}{4N}\sum_{{\bf
k},ss^{\prime}}\zeta^{{\bf k}{\bf
q},ss^{\prime}}_{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}\frac{n_{f}(\xi_{{\bf
k}+{\bf q},s^{\prime}})-n_{f}(\xi_{{\bf k},s})}{{\rm i}\omega_{m}+\xi_{{\bf
k},s}-\xi_{{\bf k}+{\bf q},s^{\prime}}},\end{aligned}$ (12)

which $n_{f}(\xi_{{\bf k},s})=[1+\exp(\xi_{{\bf k},s}/T)]^{-1}$ denotes the
Fermi distribution function. Here, the weight factor ${\zeta}^{{\bf k}{\bf
q},ss^{\prime}}_{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}$ is obtained by
[66]

$\displaystyle\begin{aligned} \zeta^{{\bf k}{\bf
q},ss^{\prime}}_{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}=\Big{[}\hat{\mathbb{I}}+s\hat{{\bf
g}}_{{\bf
k}}\cdot\hat{\bm{\sigma}}\Big{]}_{\sigma_{1}\sigma_{2}}\Big{[}\hat{\mathbb{I}}+s^{\prime}\hat{{\bf
g}}_{{\bf k}+{\bf
q}}\cdot\hat{\bm{\sigma}}\Big{]}_{\sigma_{3}\sigma_{4}}.\end{aligned}$ (13)

Next, we include the effects of the on-site Hubbard interaction by
perturbatively renormalizing the bare susceptibility. Within RPA, the elements
of the dressed spin susceptibility matrix are given by the following Dyson
equation

$\displaystyle\begin{aligned}
\hat{\chi}^{\sigma_{3}\sigma_{4}}_{\sigma_{1}\sigma_{2}}=\hat{\chi}^{(0)\sigma_{3}\sigma_{4}}_{\sigma_{1}\sigma_{2}}+\sum_{\\{\alpha_{i}\\}}\hat{\chi}^{(0)\alpha_{1}\alpha_{2}}_{\sigma_{1}\sigma_{2}}\hat{U}^{\alpha_{3}\alpha_{4}}_{\alpha_{1}\alpha_{2}}\hat{\chi}^{\sigma_{3}\sigma_{4}}_{\alpha_{3}\alpha_{4}},\end{aligned}$
(14)

where $\alpha_{i}$ represents spin indices on the internal lines of the
Feynman diagrams. Furthermore, $\hat{U}$ is the matrix of the bare electron-
electron interactions in spin space, whose non-zero elements are
$\hat{U}^{\downarrow\downarrow}_{\uparrow\uparrow}=\hat{U}^{\uparrow\uparrow}_{\downarrow\downarrow}=-U$,
and
$\hat{U}^{\uparrow\downarrow}_{\downarrow\uparrow}=\hat{U}^{\downarrow\uparrow}_{\uparrow\downarrow}=+U$.
It is worth to emphasize here again that in the presence of ASOC spin-flip
scattering is present. Thus, both the bubble (screening) and ladder (exchange)
diagrams need to be summed up to obtain the RPA spin susceptibility. The
longitudinal and transverse components of the dressed spin susceptibility are
$\chi^{\perp}=\sum_{\sigma}(\hat{\chi}^{\sigma\sigma}_{\sigma\sigma}-\hat{\chi}^{\sigma\bar{\sigma}}_{\sigma\bar{\sigma}})$,
and
$\chi^{\parallel}=\sum_{\sigma}\hat{\chi}^{\sigma\sigma}_{\bar{\sigma}\bar{\sigma}}$,
respectively. The former corresponds to the bubble diagrams while the latter
is related to the ladder diagrams.

###  III.2 Results

After this brief review of the RPA formalism, we present in the following
results for the static ($\omega=0$) spin susceptibility. Throughout this
section, we keep the interaction strength $U=1.5$ and temperature $T=0.03$
fixed.

In Fig. 2, the influence of the next-nearest neighbor hopping $t_{2}$ on the
spin structure in the absence of spin-orbit coupling is shown.

Figure 2:  Static ($\omega=0$) RPA spin susceptibilities in the absence of SOC
at zero doping. Since $V_{\rm SOC}=0$, the longitudinal and transverse
components of the spin susceptibility are identical.

Note that in this case SU(2) spin-rotational symmetry is preserved. For
$t_{2}=0$, the peak of the spin susceptibility occurs in the $\bm{\Gamma}-{\bf
K}$ path near the ${\bf K}$ point. Increasing the amplitude of the next-
nearest neighbor hopping leads to a warping of the Fermi surface, and the
position of the maximum is slightly shifted towards the ${\bf K}$ point. The
larger $t_{2}$ is, the stronger the spin fluctuations in the $\bm{\Gamma}-{\bf
K}$ path are with the secondary peak shifting towards the ${{\bf{M}}}$ point.
In all cases, the spin fluctuations are incommensurate and for $U=1.5$ the
system is not magnetically ordered. Moreover, for all values of the next-
nearest neighbor hopping, the amplitude of the spin susceptibility at the
$\bm{\Gamma}$ point is almost zero. I.e., pure ferromagnetic spin fluctuations
are absent.

Next, we show in Fig. 3 how the spin orbit coupling $V_{\rm SOC}$ influences
the spin fluctuations in the absence of next-nearest neighbor hopping,
$t_{2}=0$, and at zero doping, $\delta=0$.

Figure 3:  Static ($\omega=0$) longitudinal (upper panel) and transverse
(lower panel) spin susceptibilities in RPA for $t_{1}=1$, and $t_{2}=0$ at
zero doping with $U=1.5$. The red, blue, and green lines represent $V_{\rm
SOC}=0$, $0.3$, and $0.6$, respectively. Note that for $V_{\rm SOC}=0$, the
longitudinal and transverse components of the spin susceptibility are the same
due to SU$(2)$ symmetry.

For $V_{\rm SOC}\neq 0$, the SU(2) spin-rotational symmetry is broken and
transverse and longitudinal components are different. With increasing spin-
orbit coupling, we find that the spin susceptibility is, overall, suppressed.
This reduction can be understood as follows: Rashba SOC leads to spin-flip
processes and breaks the magnetic susceptibility into parallel and opposite
spin channels destabilizing magnetic order. The exception are ferromagnetic
fluctuations which are not affected by spin-flip scattering so that $V_{\rm
SOC}$ can even increase the magnitude of the susceptibility at the $\Gamma$
point. For spin-mediated superconducting pairing we might therefore expect
that triplet pairing is enhanced [67]. This point is consistent with the
results of Ref. [57].

If the magnitude of $t_{2}$ is increased for a fixed strength of SOC, then
magnetic fluctuations are enhanced which ultimately will bring the system for
large $t_{2}$ close to a spin density wave (SDW) state, see Fig. 4.

Figure 4:  Static ($\omega=0$) longitudinal and transverse spin
susceptibilities in RPA for $V_{\rm SOC}=0.3$ at zero doping for different
$t_{2}$.

Ferromagnetic fluctuations are also enhanced which can potentially favour a
spin mediated pairing in the triplet channel. Note, however, that all the
fluctuations remain incommensurate with the lattice.

Finally, we also consider the effect of hole versus electron doping, see Fig.
5.

Figure 5:  Static ($\omega=0$) longitudinal and transverse spin
susceptibilities in RPA for $t_{2}=0$ and $V_{\rm SOC}=0.3$ at different
levels of hole (left column) and electron (right column) doping.

In the hole doped case, an increase of the doping level simply leads to an
overall suppression of spin fluctuations both in the longitudinal and in the
transverse component. The evolution of the susceptibilities is more intricate
in the electron-doped case because of the presence of a van-Hove singularity,
see Fig. 1. For doping levels close to the van-Hove singularity the topology
of the Fermi surface changes and a sharp ferromagnetic peak develops at the
$\bm{\Gamma}$ point. This result is consistent with Ref. [68]. If the strength
of the on-site Hubbard interaction is increased further, the system will
therefore be driven into a ferromagnetically ordered state—at least at the RPA
level of approximation. Close to this ferromagnetic state we might expect the
potential for an enhanced triplet pairing. Away from the van-Hove
singularities the spin fluctuations are always incommensurate and short
ranged.

##  IV Superconductivity

Having obtained the spin susceptibilities we are now in a position to consider
the pairing of electrons and the symmetries of the superconducting gap
function.

###  IV.1 Pairing within RPA spin fluctuation theory

In the channel of Cooper pairs, we consider a pair of electrons with momenta
and spins $({\bf k}^{\prime}\sigma_{3},-{\bf k}^{\prime}\sigma_{4})$ which
interact and scatter with new momenta and spins $({\bf k}\sigma_{1},-{\bf
k}\sigma_{2})$. Within RPA, the interaction between the electrons involved in
the scattering process will be described by

${\cal H}^{\rm RPA}_{\rm int}({\bf k},{\bf
k}^{\prime})=\frac{1}{N}\\!\\!\sum_{{\bf k}{\bf
k}^{\prime},\\{\sigma_{i}\\}}\hat{\Gamma}^{\sigma_{3}\sigma_{4}}_{\sigma_{1}\sigma_{2}}({\bf
k},{\bf k}^{\prime})c^{\dagger}_{{\bf k}\sigma_{1}}c^{\dagger}_{-{\bf
k}\sigma_{2}}c_{-{\bf k}^{\prime}\sigma_{4}}c_{{\bf k}^{\prime}\sigma_{3}}.$
(15)

In the presence of ASOC, the effective vertex
$\hat{\Gamma}^{\sigma_{3}\sigma_{4}}_{\sigma_{1}\sigma_{2}}({\bf k},{\bf
k}^{\prime})$ includes the contributions of longitudinal (screening) and
transverse (exchange) interactions. In order to find the symmetries of the
superconducting instabilities, one can project the superconducting gap
function into separate singlet and triplet channels. In the singlet channel,
the vertex function is symmetric with respect to ${\bf k}$ and ${\bf
k}^{\prime}$, and is defined by

$\hat{\Gamma}^{\rm S}_{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}({\bf k},{\bf
k}^{\prime})=\frac{1}{2}\Big{[}\hat{U}+\frac{3}{2}\hat{U}\hat{\chi}_{{\bf
k}-{\bf k}^{\prime}}\hat{U}+\frac{3}{2}\hat{U}\hat{\chi}_{{\bf k}+{\bf
k}^{\prime}}\hat{U}\Big{]},$ (16)

while, in the triplet channel it is expressed as

$\hat{\Gamma}^{\rm T}_{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}({\bf k},{\bf
k}^{\prime})=-\frac{1}{2}\Big{[}\hat{U}\hat{\chi}_{{\bf k}-{\bf
k}^{\prime}}\hat{U}-\hat{U}\hat{\chi}_{{\bf k}+{\bf
k}^{\prime}}\hat{U}\Big{]},$ (17)

which has been antisymmetrized. It should be mentioned that we only consider
the static form of the vertex function ($\omega=0$); any frequency-dependence
is neglected. Considering only the intra-band Cooper pairings within band
basis, one can use the following transformation to find the effective pairing
interaction between two electrons on the Fermi surface

$\displaystyle V^{{\rm S}/{\rm T}}_{ss^{\prime}}({\bf k},{\bf k}^{\prime})$
(18) $\displaystyle=\sum_{\\{\sigma_{i}\\}}\hat{\Gamma}^{{\rm S}/{\rm
T}}_{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}({\bf k},{\bf
k}^{\prime})\Lambda^{s*}_{{\bf k}\sigma_{1}}\Lambda^{s*}_{-{\bf
k}\sigma_{2}}\Lambda^{s^{\prime}}_{-{\bf
k}^{\prime}\sigma_{4}}\Lambda^{s^{\prime}}_{{\bf k}^{\prime}\sigma_{3}},$

with $\Lambda^{s}_{{\bf k}\sigma}=\langle{\bf k},\sigma|{\bf k},s\rangle$
connecting the state $|{\bf k},\sigma\rangle$ in spin space with $|{\bf
k},s\rangle$ in the band basis.

Using Eq. (18) as the effective interaction within the BCS theory of
superconductivity, the superconducting gap function in the singlet and triplet
channels for band $s$ is then given by [69]

$\Delta_{{\bf k}s}=-\frac{1}{N}\sum_{{\bf k}^{\prime}s^{\prime}}V^{{\rm
S}/{\rm T}}_{ss^{\prime}}\frac{\Delta_{{\bf k}^{\prime}s^{\prime}}}{2E_{{\bf
k}^{\prime}s^{\prime}}}\tanh\Big{(}\frac{E_{{\bf
k}^{\prime}s^{\prime}}}{2k_{B}T}\Big{)},$ (19)

where $E_{{\bf k}s}=\sqrt{\xi^{2}_{{\bf k}s}+\Delta^{2}_{{\bf k}s}}$ defines
the energy dispersion of the superconducting quasiparticles. Close to the
superconducting critical temperature ($T\rightarrow T_{c}$), the quasiparticle
gap can be neglected leading to $E_{{\bf k}s}=|\xi_{{\bf k}s}|$. The
linearized gap equation can then be reduced to an eigenvalue problem to
determine the leading and sub-leading superconducting instabilities of the
system [22, 41]

$\lambda\Delta_{{\bf k}s}=-\frac{1}{(2\pi)^{2}}\sum_{s^{\prime}}\oint_{\rm
FS}\frac{dk^{\prime}_{\parallel}}{v^{\rm F}_{{\bf
k}^{\prime}s^{\prime}}}V^{{\rm S}/{\rm T}}_{ss^{\prime}}\Delta_{{\bf
k}^{\prime}s^{\prime}}.$ (20)

In this equation, $dk^{\prime}_{\parallel}$ is the tangential component of the
momentum differential on the Fermi surface and $v^{\rm F}_{{\bf
k}s}=|\nabla\xi_{{\bf k}s}|$ is the Fermi velocity of band $s$. The largest
eigenvalue of Eq. (20) specifies the leading superconducting order parameter
with $T_{c}\propto\exp(-1/\lambda)$. In addition, the normalized eigenvector
corresponding to the largest eigenvalue determines the momentum dependence of
the superconducting gap function and its nodal structure. Besides, we can
decompose the gap function into an amplitude $\Delta_{0}$ and a dimensionless
symmetry gap function $\phi_{{\bf k}}$ for all possible superconducting
pairings as $\Delta_{{\bf k}}=\Delta_{0}\phi_{{\bf k}}$. The triangular
lattice has $D_{6}$ point group symmetry, containing four one dimensional (1D)
($A_{1}$, $A_{2}$, $B_{1}$, $B_{2}$), and two 2D ($E_{1}$, $E_{2}$)
irreducible representations (Irreps.). $A_{1}$, $A_{2}$, and $E_{2}$ generate
Cooper pairs in the singlet channel, and $B_{1}$, $B_{2}$, $E_{1}$ give rise
to pairing in the triplet channel [70]. The corresponding momentum dependency
of gap functions ($\phi_{{\bf k}}$) for each irrep. are shown in table. 1.

Figure 6:  The normalized superconducting gap functions with (left)
$d_{x^{2}-y^{2}}$, and (right) $d_{xy}$-wave symmetries for $t_{1}=1$,
$t_{2}=0$, and $V_{\rm SOC}=0$, at doping level $\delta=0$. The system is
expected to pick one of the degenerate superconducting ground states
$d_{x^{2}-y^{2}}\pm{\rm i}d_{xy}$ leading to a spontaneously broken time
reversal symmetry.  Table 1: Character table of the superconducting gap
functions for different irreducible representations of the point group $D_{6}$
[70].

Irrep. | Symmetry | $\phi_{{\bf k}}$
---|---|---
$A_{1}$ | ext.$s$-wave | $\cos k_{x}+2\cos\frac{k_{x}}{2}\cos\frac{\sqrt{3}k_{y}}{2}$
$B_{1}$ |  $f_{x(x^{2}-3y^{2})}$-wave | $\sin k_{x}-2\sin\frac{k_{x}}{2}\cos\frac{\sqrt{3}k_{y}}{2}$
$E_{1}$ |  $p_{x}$-wave | $\sin k_{x}+2\sin\frac{k_{x}}{2}\cos\frac{\sqrt{3}k_{y}}{2}$
$E_{1}$ |  $p_{y}$-wave | $\cos\frac{k_{x}}{2}\sin\frac{\sqrt{3}k_{y}}{2}$
$E_{2}$ | $d_{x^{2}-y^{2}}$ | $\cos k_{x}-2\cos\frac{k_{x}}{2}\cos\frac{\sqrt{3}k_{y}}{2}$
$E_{2}$ | $d_{xy}$ | $\sin\frac{k_{x}}{2}\sin\frac{\sqrt{3}k_{y}}{2}$

For each irrep., the eigenvalue $\lambda$ of Eq. (20) forms a $2\times 2$
matrix ($\lambda^{ss^{\prime}}_{l}$). The largest eigenvalue of
$\lambda^{ss^{\prime}}_{l}$ gives [71]

$\lambda^{ss^{\prime}}_{l}=-\frac{\int_{{\rm
FS}_{s}}\frac{dk_{\parallel}}{v^{\rm F}_{{\bf k}s}}\int_{{\rm
FS}_{s^{\prime}}}\frac{dk^{\prime}_{\parallel}}{v^{\rm F}_{{\bf
k}^{\prime}s^{\prime}}}\phi^{l}_{{\bf k}}V^{{\rm S}/{\rm
T}}_{ss^{\prime}}({\bf k},{\bf k}^{\prime})\phi^{l}_{{\bf
k}^{\prime}}}{2\pi^{2}\int_{{\rm
FS}_{s^{\prime}}}\frac{dk^{\prime}_{\parallel}}{v^{\rm F}_{{\bf
k}^{\prime}s^{\prime}}}[\phi^{l}_{{\bf k}^{\prime}}]^{2}},$ (21)

where $l$ is running over the allowed irreps. for point group $D_{6}$. The
diagonal and off-diagonal elements of matrix $\lambda^{ss^{\prime}}_{l}$
relate to the Cooper pairing between to electrons in the same and opposite
helical bands, respectively. The dominant superconducting order parameter in
singlet and triplet channels for irrep. $l$ are given by the largest
eigenvalue of the $\lambda^{ss^{\prime}}_{l}$ matrix, which is denoted as
$\lambda^{\rm eff}_{l}$. It is possible to estimate the superconducting
critical temperature for each irrep. by
$T^{l}_{c}=1.13\omega_{c}\exp(-1/\lambda^{\rm eff}_{l})$, in which
$\omega_{c}$ denotes the cut-off energy and is of the order of the bandwidth.

Figure 7:  The normalized superconducting gap functions in the singlet (upper
row) and triplet (lower row) channels for $t_{1}=1$, $t_{2}=0$, and $V_{\rm
SOC}=0.3$, at doping level $\delta=0$. Rashba SOC violates parity, which
results in parity mixing in the superconducting order parameter. Due to
spontaneously broken time-reversal symmetry, the superconducting gap function
in the singlet and triplet channels are $d_{x^{2}-y^{2}}\pm{\rm i}d_{xy}$ and
$p_{x}\pm{\rm i}p_{y}$, respectively.  Figure 8:  Magnitude of the eigenvalue
(21) for different superconducting channels as a function of doping for (a)
$t_{2}=V_{\rm SOC}=0$ and (b) $t_{2}=0$, $V_{\rm SOC}=0.3$ with $U=2$. The red
region denotes the doping levels where magnetic order competes with
superconductivity.  Figure 9:  Critical line separating an SDW phase from a
superconducting/metallic phase for (a) $t_{2}=V_{\rm SOC}=0$ and (b)
$t_{2}=0$, $V_{\rm SOC}=0.3$. The critical line is determined by finding the
Coulomb interaction $U_{c}$ for each doping level at which the RPA spin
susceptibility diverges.

###  IV.2 Results

First, we consider the symmetry of the gap function. Fig. 6 shows the dominant
superconducting pairing for $t_{2}=0$, $V_{\rm SOC}=0$, and doping $\delta=0$.

Here, the normalized magnitude of the superconducting gap function on the
Fermi surface is depicted. Since Rashba coupling is absent, spin is a good
quantum number and we find that the leading superconducting instability is in
the singlet channel with degenerate $d_{x^{2}-y^{2}}$\- and $d_{xy}$-wave
pairings. This degeneracy originates from the C3v symmetry of the triangular
lattice. Minimization of the energy gives a chiral $d+id$ pairing, which is a
consequence of the geometrical frustration of the triangular lattice [72].
Moreover, we find that this result remains true irrespective of the doping
level. If the system picks one of the degenerate states
$d_{x^{2}-y^{2}}\pm{\rm i}d_{xy}$ then TRS will be spontaneously broken. This
type of chiral pairing belongs to the $C$ class of topological
superconductors, which are characterized by a $\mathbb{Z}$ topological index.

Next, we consider the influence of the SOC on the superconducting gap
function. Since spin is then no longer a good quantum number we expect that
singlet and triplet components are mixed in general. Indeed for $t_{2}=0$,
$V_{\rm SOC}=0.3$, and $U=2$ we find both d-wave and p-wave superconductivity,
see Fig. 7.

As in the case without SOC, the pairing symmetry in the singlet channel has
the chiral form $d_{x^{2}-y^{2}}\pm{\rm i}d_{xy}$. Since Rashba SOC also
breaks parity, we find that the even-parity spin singlet and odd-parity spin-
triplet pairings are mixed together. Our calculations show that the leading
pairing in the triplet channel has both $p_{x}$ and $p_{y}$ symmetries at the
same eigenvalue $\lambda$. In a similar manner to the singlet channel, the
minimization of energy in the triplet channel will thus show a chiral texture
of the form $p_{x}+{\rm i}p_{y}$. This means that the symmetry of Cooper pairs
in the triplet channel belongs to $E_{1}$ irrep. In contrast to the case of a
square lattice, in which the triplet part of the superconducting gap has a
helical characteristic [73], the spontaneous breaking of TRS in the
superconducting gap function for a triangular lattice generates a chiral
triplet order parameter. As a result, the precise form of the superconducting
gap function is $(d_{x^{2}-y^{2}}+{\rm i}d_{xy})+(p_{x}+{\rm i}p_{y})$-wave.
Note, however, that our calculations show that for the parameters chosen in
Fig. 7 the singlet pairing is dominant.

More generally speaking, the magnitude of the eigenvalues in the singlet and
triplet channels will depend on the hopping amplitude $t_{2}$, the strength of
SOC, the Hubbard interaction, and the doping level. As shown in Fig. 8, the
d-wave instability is always dominant in the absence of SOC. With SOC, on the
other hand, we find that a d-wave instability is dominant in the hole-doped
case while a p-wave instability dominates for electron doping close to the SDW
instability. This result is consistent with our earlier results for the spin
susceptibility: For electron doping levels in between the two VHSs, the
triplet pairing is dominant due to the presence of strong ferromagnetic
fluctuations. We stress again that in this case the singlet and triplet
channels will always be mixed and only the relative level of admixture is
changing with doping. Using that $T_{c}\propto\omega_{c}\exp(-1/\lambda_{\rm
eff})$ we can also obtain a rough estimate for the transition temperature
$T_{c}$ where the cutoff frequency $\omega_{c}$ is of the order of the
bandwidth. If we assume a bandwidth in the range of 1 eV we need a
$\lambda_{\rm eff}\geq 0.1$ to get a transition temperature in the range of a
few Kelvin. We see in Fig. 8 that such values for $\lambda_{\rm eff}$ are
indeed reached for the leading superconducting pairings close to the regime
where the system shows an SDW instability.

We can also use the RPA results to determine the ground state phase diagram by
finding the critical Hubbard interaction strength $U_{c}$ at which the spin
susceptibility diverges which indicates an instability towards spin-density
wave (SDW) order. In this picture there is a regime $U<U_{c}$ where
superconductivity is established—note that we are discussing the ground state
phase diagram: the critical temperature $T_{c}$ might be extremely small and
not accessible in experiment which is why we denote this region by SC/metal in
Fig. 9—and a regime $U>U_{c}$ where magnetic order sets in. The critical line
separating the two phases obtained in this manner is shown in Fig. 9.

In both cases, there is a minimum at a filling which corresponds to the van-
Hove singularity in the DOS, see Fig. 1.

##  V Quasiparticle interference

In this section we want to discuss an experimental technique which can
potentially be used to investigate the superconducting state. We note that
quasiparticle interference is an STM technique which does require access to
the surface of the bilayer and is therefore not suitable for a heterostructure
with both a top and a bottom gate. The results presented in this section are
therefore, for now, purely theoretical and applicable for a hypothetical
system where the layer realizing the 2d triangular Hubbard model is accessible
by STM.

Figure 10:  (a) Spectral function of the superconducting state in the absence
of SOC ($V_{\rm SOC}=0$) for $t_{1}=1$, $\Delta_{\rm S}=0.2$, and $\Delta_{\rm
T}=0$. (b) Total QPI for both charge ($\eta=0$), and spin ($\eta=x$, $y$, and
$z$) channels at $\omega=0.1$. Note that because of preserved SU$(2)$
symmetry, all components of the QPI are the same.

###  V.1 Method

First, we briefly review how to calculate the relevant quantities measured
experimentally. The STM-based quasiparticle interference (QPI) technique is
one of the most powerful methods to discriminate between different gap
symmetries in a wide variety of superconducting materials. In the presence of
magnetic or non-magnetic impurities, the local density of states (LDOS) will
be modulated. Consequently, the superposition of incident and scattered
quasiparticles generate interference patterns, carrying information about the
Fermi surface and the superconducting gap function. The QPI pattern is
obtained from the full Green’s function, which is determined by considering
the effect of scattering from random charge and spin impurities at the
surface. The scattering potential is given by [69]

${\cal H}_{\rm imp}=\sum_{{\bf k}{\bf q}\eta}V^{\eta}_{{\bf
q}}\hat{S}_{\eta}\Psi^{\dagger}_{{\bf q}+{\bf k}}\hat{\rho}_{\eta}\Psi_{{\bf
q}},$ (22)

where $\Psi^{\dagger}_{{\bf k}}=(c^{\dagger}_{{\bf
k}\uparrow},c^{\dagger}_{{\bf k}\downarrow},c_{-{\bf k}\uparrow},c_{-{\bf
k}\downarrow})$ represents the Nambu space describing the electron-hole
symmetric basis. Within this space, $\hat{\rho}_{\eta=0,x,y,z}$ represent
$4\times 4$ matrices with
$\\{\hat{\rho}_{\eta}\\}=(\tau_{z}\otimes\sigma_{0},\tau_{0}\otimes\sigma_{x},\tau_{z}\otimes\sigma_{y},\tau_{0}\otimes\sigma_{z})$,
in which $\tau$ and $\sigma$ are Pauli matrices in particle-hole and spin
spaces, respectively. It should be mentioned that $\tau_{0}$ and $\sigma_{0}$
are $2\times 2$ identity matrices. Moreover, $\\{S_{\eta}\\}=(1,{\bf S})$ with
${\bf S}$ denoting the spin of the magnetic impurity. The potential of a non-
magnetic impurity is given by $V^{\eta=0}_{{\bf q}}$, while a magnetic
impurity potential is given by an exchange scattering ($V^{\eta\neq 0}_{{\bf
q}}=V_{\rm ex}$) from an impurity with spin ${\bf S}$. Within the second order
Born approximation, the Fourier transform of the modulation of LDOS at
frequency $\omega$ is [74]

$dN_{\eta\eta^{\prime}}({\bf q},\omega)=-\frac{1}{\pi}V^{\eta^{\prime}}_{{\bf
q}}{\rm Im}\Big{[}\Lambda_{\eta\eta^{\prime}}({\bf q},\omega)\Big{]},$ (23)

in which

$\Lambda_{\eta\eta^{\prime}}({\bf q},\omega)=\frac{1}{N}\sum_{{\bf k}}{\rm
Tr}_{\sigma}\Big{[}\hat{\cal P}\hat{\rho}_{\eta}\hat{G}_{{\bf k}}(\rm
i\omega)\hat{\rho}_{\eta^{\prime}}\hat{G}_{{\bf k}+{\bf q}}(\rm
i\omega)\Big{]}_{\rm i\omega\rightarrow\omega+{\rm i}0^{+}}.$ (24)

In this equation, $\hat{\cal P}=(\tau_{0}+\tau_{z})\otimes\sigma_{0}/2$ is a
projection operator. In addition, $\hat{G}_{{\bf k}}(\rm i\omega)=[\rm
i\omega-\hat{\cal H}^{\rm BdG}_{{\bf k}}]^{-1}$ is the Matsubara Green’s
function of the superconducting state, in which the Boguliubov de-Gennes
Hamiltonian is given by

$\hat{\cal H}^{\rm BdG}_{{\bf k}}=\begin{bmatrix}\hat{\cal H}^{0}_{{\bf
k}}&\hat{\Delta}_{{\bf k}}\\\ \hat{\Delta}^{\dagger}_{{\bf k}}&-\hat{\cal
H}^{0*}_{-{\bf k}}\end{bmatrix}.$ (25)

Here, $\hat{\cal H}^{0}_{{\bf k}}$ represents the non-interacting Hamiltonian
of the normal state. Moreover, the momentum resolved spectral function is
represented by ${\cal A}^{\rm ret}_{{\bf k}}(\omega)=-2{\rm Im}[{\rm
Tr}_{\sigma}\hat{G}_{{\bf k}}({\rm i}\omega)]_{{\rm
i\omega}\rightarrow\omega+{\rm i}0^{+}}$ and $\hat{\Delta}_{{\bf k}}={\rm
i}[\psi_{{\bf k}}\hat{\sigma}_{0}+{\bf d}_{{\bf
k}}\cdot\hat{\bm{\sigma}}]\hat{\sigma}_{y}$ is the matrix of superconducting
gap functions in spin space. Here, $\psi_{{\bf k}}$, and ${\bf d}_{{\bf k}}$
represent the singlet and triplet parts of the superconducting gap function
with amplitudes $\Delta_{\rm S}$ and $\Delta_{\rm T}$, respectively, whose
explicit forms are given by table. 1. Spontaneous breaking of time-reversal
symmetry (TRS) in a triangular lattice generates a degenerate superconducting
ground state for both singlet and triplet channels. The projection of the
superconducting gap function then gives rise to

$\Delta_{{\bf k}s}=\psi_{{\bf k}}+s|{\bf d}_{{\bf k}}|,$ (26)

in which $\Delta_{{\bf k}s}$ denotes the magnitude of the gap function on the
Fermi surface of band $s$. In this picture, the energy dispersion of
quasiparticles in band $s$ is given by $E_{{\bf k}s}=\sqrt{\xi^{2}_{{\bf
k}s}+\Delta^{2}_{{\bf k}s}}$. Using Eq. (24), the QPI patterns for charge and
spin channels can be shown to be given by [75]

$\displaystyle\begin{aligned} \Lambda_{00}({\bf
q},\omega)&=\frac{1}{N}\sum_{{\bf k}ss^{\prime}}[1+ss^{\prime}(\hat{{\bf
g}}_{{\bf k}}\cdot\hat{{\bf g}}_{{\bf k}+{\bf q}})]{\cal K}^{{\bf k}{\bf
q}}_{ss^{\prime}}(\omega),\\\ \Lambda_{uu}({\bf
q},\omega)&=\frac{1}{N}\sum_{{\bf k}ss^{\prime}}[1-ss^{\prime}(\hat{{\bf
g}}_{{\bf k}}\cdot\hat{{\bf g}}_{{\bf k}+{\bf q}}-2\hat{{\bf g}}^{u}_{{\bf
k}}\hat{{\bf g}}^{u}_{{\bf k}+{\bf q}})]{\cal K}^{{\bf k}{\bf
q}}_{ss^{\prime}}(\omega),\end{aligned}$ (27)

where the index $u=(x,y,z)$ is running over the spatial components of the
spin. In addition, the QPI kernel function ${\cal K}^{{\bf k}{\bf
q}}_{ss^{\prime}}(\omega)$ is given by

${\cal K}^{{\bf k}{\bf q}}_{ss^{\prime}}(\omega)=\frac{(\omega+\xi_{{\bf
k}s})(\omega+\xi_{{\bf k}+{\bf q}s^{\prime}})-\Delta_{{\bf k}s}\Delta_{{\bf
k}+{\bf q}s^{\prime}}}{(\omega^{2}-E^{2}_{{\bf k}s})(\omega^{2}-E^{2}_{{\bf
k}+{\bf q}s^{\prime}})}.$ (28)

Figure 11:  (a) Spectral function of the superconducting state in the presence
of SOC ($V_{\rm SOC}=0.3$) for $t_{1}=1$, $\Delta_{\rm S}=0.2$, and
$\Delta_{\rm T}=0.05$ at $\omega=0.1$. The areas with high spectral weight
dominate the quasiparticle scattering events. (b) Total QPI pattern
$\Lambda_{00}({\bf q},\omega)$ in the charge channel carrying the crystal
symmetry. (c,d) QPI pattern in the spin channel due to the scattering by a
magnetic impurity with spins $S_{x}$ and $S_{y}$, respectively. Since SOC
breaks the SU$(2)$ symmetry, $\Lambda_{xx}({\bf q},\omega)$ and
$\Lambda_{yy}({\bf q},\omega)$ do not follow the crystal point group symmetry.
(e) The interference pattern $\Lambda_{zz}({\bf q},\omega)$ corresponding to a
magnetic impurity with spin $S_{z}$ shows again the crystal symmetry.

###  V.2 Results

In Fig. 10(a), an intensity plot for the superconducting spectral function at
$\omega=0.1$ for $t_{2}=0$, $V_{\rm SOC}=0$, $\Delta_{\rm S}=0.2$, and
$\Delta_{\rm T}=0$ is shown.

In this case, the gap function has the chiral $d_{x^{2}-y^{2}}+{\rm
i}d_{xy}$-wave symmetry. The spectral function shows that the areas around the
corners of the hexagon have the largest DOS. Therefore, the scattering wave-
vectors connecting these points are playing the most important role in
determining the QPI pattern. Because of the preserved SU$(2)$ symmetry, all
different components of QPI, including scattering from non-magnetic and
magnetic impurities, have the same shapes. Fig. 10(b) shows the QPI pattern
$\Lambda_{\eta\eta}({\bf q},\omega=0.1)$ ($\eta=\\{0,x,y,z\\}$). The pattern
is consistent with the point-group symmetry of the lattice. Moreover, since
the superconducting gap function has a fully gapped structure, the QPI pattern
resembles the shape of the Fermi surface, including some dominant peaks at the
corners of the hexagon.

Finally, in Figs. 11, the superconducting spectral function and the QPI
patterns generated by non-magnetic and magnetic impurities in the presence of
SOC are reported. Fig. 11(a) shows the spectral function for $t_{2}=0$,
$V_{\rm SOC}=0.3$, $\Delta_{\rm S}=0.2$, and $\Delta_{\rm T}=0.05$ at
frequency $\omega=0.1$. In this situation, the Cooper pairs have the symmetry
of a $(d_{x^{2}-y^{2}}+{\rm i}d_{xy})+(p_{x}+{\rm i}p_{y})$-wave. This plot
carries the information about the Fermi surface together with information
about the spectral weights. The QPI pattern for the scattering of the
superconducting quasiparticles from a non-magnetic impurity is shown in Fig.
11(b). The dominant sharp peaks originate from intra-band scattering. However,
there are also some sub-leading peaks around the $\bm{\Gamma}$ point, which
can be attributed to inter-band scattering. Figs. 11(c,d) show
$\Lambda_{xx}({\bf q},\omega=0.1)$, and $\Lambda_{yy}({\bf q},\omega=0.1)$,
respectively. We observe that these two QPI patterns have more complicated
forms than $\Lambda_{00}({\bf q},\omega=0.1)$. Eq. (27) analytically justifies
this complex form. Since the Rashba SOC only has in-plane components, the QPI
patterns generated by magnetic impurities with spins $S_{x}$ and $S_{y}$ will
not obey the symmetry of the crystal point-group. Lastly, Fig. 11(e) depicts
$\Lambda_{zz}({\bf q},\omega=0.1)$. Based on Eq. (27) it is clear that this
QPI pattern—resulting from the scattering of superconducting quasiparticles by
a magnetic impurity of spin $S_{z}$—shows again the D3 symmetry of the system.

##  VI Summary and Conclusions

We have studied a Rashba-Hubbard Hamiltonian on a triangular lattice as a
model for a twisted bilayer of TMD which has a maximum LDOS at ${\cal
R}^{M}_{M}$ high-symmetry positions. In our model, the Rashba SOC is a
consequence of depositing the 2D thin layers on a thick substrate.

In the absence of second neighbor hopping and Rashba SOC, we find
incommensurate magnetic fluctuations consistent with earlier results in the
literature. We then observe that increasing the amplitude of the next-nearest
neighbor hopping leads to a reduction of the incommensurate spin fluctuations
around the ${\bf K}$ point while amplifying those around the ${{\bf{M}}}$
position.

Looking at electron pairings mediated by these spin fluctuations, we find that
in the case without SOC a chiral spin singlet pairing with
$d_{x^{2}-y^{2}}+{\rm i}d_{xy}$ symmetry dominates. Upon introducing a Rashba
SOC, dominant ferromagnetic fluctuations in the longitudinal susceptibility
can arise which are highly sensitive to the level of filling. These
fluctuations appear when both hole- and electron-like Fermi surfaces exist. In
the case with SOC, the superconducting fluctuations become predominantly of
$(d_{x^{2}-y^{2}}+{\rm i}d_{xy})+(p_{x}+{\rm i}p_{y})$ character. A
superconducting phase with this symmetry belongs to the $C$ class of
topological superconductors which can be characterized by a $\mathbb{Z}$
topological invariant.

In order to experimentally detect gap structures in superconducting materials,
the STM-based quasiparticle interference is a powerful tool. We have therefore
also examined the QPI response both from non-magnetic and magnetic impurities.
The obtained results might serve as a guide for future experiments searching
for exotic superconducting phases in tunable bilayer materials.

## Acknowledgments

We are grateful to A. Greco, Y. Yanase, A. Akbari, S. A. Jafari, A. Romer, M.
Malakhov, and M. N. Najafi for fruitful discussions. J. S. acknowledges
support by the Natural Sciences and Engineering Research Council (NSERC,
Canada) and by the Deutsche Forschungsgemeinschaft (DFG) via Research Unit FOR
2316.

## References

  * Arora  _et al._ [2020] H. S. Arora, R. Polski, Y. Zhang, A. Thomson, Y. Choi, H. Kim, Z. Lin, I. Z. Wilson, X. Xu, J.-H. Chu, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Superconductivity in metallic twisted bilayer graphene stabilized by wse2, Nature 583, 379–384 (2020).
  * Kennes  _et al._ [2021] D. Kennes, M. Claassen, L. Xian, A. Georges, A. Millis, J. Hone, C. Dean, D. Basov, A. Pasupathy, and A. Rubio, Moiré heterostructures as a condensed-matter quantum simulator, Nature Physics 17, 1–9 (2021).
  * Das  _et al._ [2021] I. Das, X. Lu, J. Herzog-Arbeitman, Z.-D. Song, K. Watanabe, T. Taniguchi, B. A. Bernevig, and D. K. Efetov, Symmetry-broken chern insulators and rashba-like landau-level crossings in magic-angle bilayer graphene, Nature Physics 17, 710–714 (2021).
  * Li  _et al._ [2021] H. Li, U. Kumar, K. Sun, and S.-Z. Lin, Spontaneous fractional chern insulators in transition metal dichalcogenide moiré superlattices, Phys. Rev. Research 3, L032070 (2021).
  * Sheffer and Stern [2021] Y. Sheffer and A. Stern, Chiral magic-angle twisted bilayer graphene in a magnetic field: Landau level correspondence, exact wave functions, and fractional chern insulators, Phys. Rev. B 104, L121405 (2021).
  * Wu  _et al._ [2021] S. Wu, Z. Zhang, K. Watanabe, T. Taniguchi, and E. Andrei, Chern insulators, van hove singularities and topological flat bands in magic-angle twisted bilayer graphene, Nature Materials 20, 488–494 (2021).
  * Chen  _et al._ [2021] B.-b. Chen, Y. Liao, Z. Chen, O. Vafek, J. Kang, W. Li, and Z. Meng, Realization of topological mott insulator in a twisted bilayer graphene lattice model., Nature communications 12 1, 5480 (2021).
  * Cao  _et al._ [2020] Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J.-m. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunable correlated states and spin-polarized phases in twisted bilayer–bilayer graphene, Nature 583, 215–220 (2020).
  * Huang  _et al._ [2017] B. Huang, G. Clark, E. Navarro-Moratalla, D. R. Klein, R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A. McGuire, D. H. Cobden, W. Yao, D. Xiao, P. Jarillo-Herrero, and X. Xu, Layer-dependent ferromagnetism in a van der waals crystal down to the monolayer limit, Nature (London) 546, 270–273 (2017).
  * Wu  _et al._ [2019] 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).
  * Zhang  _et al._ [2020] Y. Zhang, Z. Zhan, F. Guinea, J. A. Silva-Guillen, and S. Yuan, Tuning band gaps in twisted bilayer $\mathrm{Mo}{\mathrm{s}}_{2}$, Phys. Rev. B 102, 235418 (2020).
  * Rickhaus  _et al._ [2021] P. Rickhaus, F. K. de Vries, J. Zhu, E. Portoles, G. Zheng, M. Masseroni, A. Kurzmann, T. Taniguchi, K. Watanabe, A. H. MacDonald, T. Ihn, and K. Ensslin, Correlated electron-hole state in twisted double-bilayer graphene, Science 373, 1257–1260 (2021).
  * Cao  _et al._ [2018] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature (London) 556, 80–84 (2018).
  * Song  _et al._ [2019] Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, All magic angles in twisted bilayer graphene are topological, Phys. Rev. Lett. 123, 036401 (2019).
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Origin of magic angles in twisted bilayer graphene, Phys. Rev. Lett. 122, 106405 (2019).
  * Sharma  _et al._ [2020] G. Sharma, M. Trushin, O. P. Sushkov, G. Vignale, and S. Adam, Superconductivity from collective excitations in magic-angle twisted bilayer graphene, Phys. Rev. Research 2, 022040 (2020).
  * Sboychakov  _et al._ [2020] A. O. Sboychakov, A. V. Rozhkov, A. L. Rakhmanov, and F. Nori, Spin density wave and electron nematicity in magic-angle twisted bilayer graphene, Phys. Rev. B 102, 155142 (2020).
  * Wu  _et al._ [2018a] F. Wu, A. H. MacDonald, and I. Martin, Theory of phonon-mediated superconductivity in twisted bilayer graphene, Phys. Rev. Lett. 121, 257001 (2018a).
  * Cao  _et al._ [2018] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43–50 (2018).
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, and D. K. Efetov, Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, Nature (London) 574, 653–657 (2019).
  * Cao  _et al._ [2021] Y. Cao, D. Rodan-Legrain, J. M. Park, N. F. Q. Yuan, K. Watanabe, T. Taniguchi, R. M. Fernandes, L. Fu, and P. Jarillo-Herrero, Nematicity and competing orders in superconducting magic-angle graphene, Science 372, 264–271 (2021).
  * Liu  _et al._ [2018] C.-C. Liu, L.-D. Zhang, W.-Q. Chen, and F. Yang, Chiral spin density wave and $d+id$ superconductivity in the magic-angle-twisted bilayer graphene, Phys. Rev. Lett. 121, 217001 (2018).
  * Kennes  _et al._ [2018] D. M. Kennes, J. Lischner, and C. Karrasch, Strong correlations and $d+\mathit{id}$ superconductivity in twisted bilayer graphene, Phys. Rev. B 98, 241407 (2018).
  * You and Vishwanath [2019] Y.-Z. You and A. Vishwanath, Superconductivity from valley fluctuations and approximate so(4) symmetry in a weak coupling theory of twisted bilayer graphene, npj Quantum Materials 4, 16 (2019).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Tuning superconductivity in twisted bilayer graphene, Science 363, 1059–1064 (2019).
  * Wu  _et al._ [2018b] F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Hubbard model physics in transition metal dichalcogenide moiré bands, Phys. Rev. Lett. 121, 026402 (2018b).
  * Wang  _et al._ [2020] L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. Rhodes, C. Tan, M. Claassen, D. Kennes, Y. Bai, B. Kim, K. Watanabe, T. Taniguchi, X. Zhu, J. Hone, A. Rubio, A. Pasupathy, and C. Dean, Correlated electronic phases in twisted bilayer transition metal dichalcogenides, Nature Materials 19 (2020).
  * Regan  _et al._ [2020] 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–363 (2020).
  * Tang  _et al._ [2020] 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–358 (2020).
  * Scherer  _et al._ [2021] M. M. Scherer, D. M. Kennes, and L. Classen, $\mathcal{N}=4$ chiral superconductivity in moiré transition metal dichalcogenides (2021), arXiv:2108.11406 .
  * Ruiz-Tijerina and Fal’ko [2019] D. A. Ruiz-Tijerina and V. I. Fal’ko, Interlayer hybridization and moiré superlattice minibands for electrons and excitons in heterobilayers of transition-metal dichalcogenides, Phys. Rev. B 99, 125424 (2019).
  * Schrade and Fu [2019] C. Schrade and L. Fu, Spin-valley density wave in moiré materials, Phys. Rev. B 100, 035413 (2019).
  * Zhou  _et al._ [2021] Y. Zhou, D. N. Sheng, and E.-A. Kim, Quantum phases of transition metal dichalcogenide moiré systems (2021), arXiv:2105.07008 .
  * Jung  _et al._ [2014] J. Jung, A. Raoux, Z. Qiao, and A. H. MacDonald, Ab initio theory of moiré superlattice bands in layered two-dimensional materials, Phys. Rev. B 89, 205414 (2014).
  * Zare and Mosadeq [2021] M.-H. Zare and H. Mosadeq, Spin liquid in twisted homobilayers of group-vi dichalcogenides, Phys. Rev. B 104, 115154 (2021).
  * Kiese  _et al._ [2021] D. Kiese, Y. He, C. Hickey, A. Rubio, and D. M. Kennes, Tmds as a platform for spin liquid physics: A strong coupling study of twisted bilayer wse2, arXiv e-prints , arXiv:2110.10179 (2021), arXiv:2110.10179 .
  * Pan  _et al._ [2020] H. Pan, F. Wu, and S. Das Sarma, Band topology, hubbard model, heisenberg model, and dzyaloshinskii-moriya interaction in twisted bilayer ${\mathrm{wse}}_{2}$, Phys. Rev. Research 2, 033087 (2020).
  * Bednorz and Müller [1986] J. G. Bednorz and K. A. Müller, Possible high tc superconductivity in the ba-la-cu-o system, Zeitschrift fur Physik B Condensed Matter 64, 189–193 (1986).
  * Kamihara  _et al._ [2008] Y. Kamihara, T. Watanabe, M. Hirano, and H. Hosono, Iron-based layered superconductor la[o1-xfx]feas (x= 0.05-0.12) with tc = 26 k, Journal of the American Chemical Society 130, 3296–7 (2008).
  * Zare  _et al._ [2017] M.-H. Zare, M. Biderang, and A. Akbari, Mixed-pairing superconductivity in $5d$ mott insulators with antisymmetric exchange: Application to ${\mathrm{sr}}_{2}{\mathrm{iro}}_{4}$, Phys. Rev. B 96, 205156 (2017).
  * Rømer  _et al._ [2021] A. T. Rømer, P. J. Hirschfeld, and B. M. Andersen, Superconducting state of ${\mathrm{sr}}_{2}{\mathrm{ruo}}_{4}$ in the presence of longer-range coulomb interactions, Phys. Rev. B 104, 064507 (2021).
  * Scalapino  _et al._ [1986] D. J. Scalapino, E. Loh, and J. E. Hirsch, $d$-wave pairing near a spin-density-wave instability, Phys. Rev. B 34, 8190–8192 (1986).
  * Sato and Fujimoto [2009a] M. Sato and S. Fujimoto, Topological phases of noncentrosymmetric superconductors: Edge states, majorana fermions, and non-abelian statistics, Phys. Rev. B 79, 094504 (2009a).
  * Scheurer and Schmalian [2015] M. S. Scheurer and J. Schmalian, Topological superconductivity and unconventional pairing in oxide interfaces, Nature Communications 6, 6005 (2015).
  * Sato and Ando [2017] M. Sato and Y. Ando, Topological superconductors: a review, Reports on Progress in Physics 80, 076501 (2017).
  * Chakhalian  _et al._ [2007] J. Chakhalian, J. Freeland, H. Habermeier, G. Cristiani, G. Khaliullin, M. van Veenendaal, and B. Keimer, Orbital reconstruction and covalent bonding at an oxide interface, Science (New York, N.Y.) 318, 1114–7 (2007).
  * Mizukami  _et al._ [2011] Y. Mizukami, H. Shishido, T. Shibauchi, M. Shimozawa, S. Yasumoto, D. Watanabe, M. Yamashita, H. Ikeda, T. Terashima, H. Kontani, and Y. Matsuda, Extremely strong coupling superconductivity in artificial two-dimensional kondo lattices, Nature Physics 7 (2011).
  * Shimozawa  _et al._ [2014] M. Shimozawa, S. K. Goh, R. Endo, R. Kobayashi, T. Watashige, Y. Mizukami, H. Ikeda, H. Shishido, Y. Yanase, T. Terashima, T. Shibauchi, and Y. Matsuda, Controllable rashba spin-orbit interaction in artificially engineered superlattices involving the heavy-fermion superconductor ${\mathrm{cecoin}}_{5}$, Phys. Rev. Lett. 112, 156404 (2014).
  * Naritsuka  _et al._ [2017] M. Naritsuka, T. Ishii, S. Miyake, Y. Tokiwa, R. Toda, M. Shimozawa, T. Terashima, T. Shibauchi, Y. Matsuda, and Y. Kasahara, Emergent exotic superconductivity in artificially engineered tricolor kondo superlattices, Phys. Rev. B 96, 174512 (2017).
  * Naritsuka  _et al._ [2018] M. Naritsuka, P. F. S. Rosa, Y. Luo, Y. Kasahara, Y. Tokiwa, T. Ishii, S. Miyake, T. Terashima, T. Shibauchi, F. Ronning, J. D. Thompson, and Y. Matsuda, Tuning the pairing interaction in a $d$-wave superconductor by paramagnons injected through interfaces, Phys. Rev. Lett. 120, 187002 (2018).
  * Bauer and Sigrist [2012] E. Bauer and M. Sigrist,  _Non-Centrosymmetric Superconductors: Introduction and Overview_ (Springer Berlin Heidelberg, 2012).
  * Nogaki  _et al._ [2021] K. Nogaki, A. Daido, J. Ishizuka, and Y. Yanase, Topological crystalline superconductivity in locally noncentrosymmetric ${\mathrm{cerh}}_{2}{\mathrm{as}}_{2}$, Phys. Rev. Research 3, L032071 (2021).
  * Li and Koshino [2019] Y. Li and M. Koshino, Twist-angle dependence of the proximity spin-orbit coupling in graphene on transition-metal dichalcogenides, Phys. Rev. B 99, 075438 (2019).
  * David  _et al._ [2019] A. David, P. Rakyta, A. Kormányos, and G. Burkard, Induced spin-orbit coupling in twisted graphene–transition metal dichalcogenide heterobilayers: Twistronics meets spintronics, Phys. Rev. B 100, 085412 (2019).
  * Naimer  _et al._ [2021] T. Naimer, K. Zollner, M. Gmitra, and J. Fabian, Twist-angle dependent proximity induced spin-orbit coupling in graphene/transition-metal dichalcogenide heterostructures (2021), arXiv:2108.06126 .
  * Lin  _et al._ [2021] J.-X. Lin, Y.-H. Zhang, E. Morissette, Z. Wang, S. Liu, D. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, and J. I. A. Li, Spin-orbit driven ferromagnetism at half moiré filling in magic-angle twisted bilayer graphene (2021), arXiv:2102.06566 .
  * Yokoyama  _et al._ [2007] T. Yokoyama, S. Onari, and Y. Tanaka, Enhanced triplet superconductivity in noncentrosymmetric systems, Phys. Rev. B 75, 172511 (2007).
  * Greco and Schnyder [2018] A. Greco and A. P. Schnyder, Mechanism for unconventional superconductivity in the hole-doped rashba-hubbard model, Phys. Rev. Lett. 120, 177002 (2018).
  * Bauer  _et al._ [2004] E. Bauer, G. Hilscher, H. Michor, C. Paul, E. W. Scheidt, A. Gribanov, Y. Seropegin, H. Noël, M. Sigrist, and P. Rogl, Heavy fermion superconductivity and magnetic order in noncentrosymmetric ${\mathrm{c}\mathrm{e}\mathrm{p}\mathrm{t}}_{3}\mathrm{S}\mathrm{i}$, Phys. Rev. Lett. 92, 027003 (2004).
  * Sato and Fujimoto [2009b] M. Sato and S. Fujimoto, Topological phases of noncentrosymmetric superconductors: Edge states, majorana fermions, and non-abelian statistics, Phys. Rev. B 79, 094504 (2009b).
  * Biderang  _et al._ [2018] M. Biderang, H. Yavari, M.-H. Zare, P. Thalmeier, and A. Akbari, Edge currents as a probe of the strongly spin-polarized topological noncentrosymmetric superconductors, Phys. Rev. B 98, 014524 (2018).
  * Xian  _et al._ [2019] L. Xian, D. M. Kennes, N. Tancogne-Dejean, M. Altarelli, and A. Rubio, Multiflat bands and strong correlations in twisted bilayer boron nitride: Doping-induced correlated insulator and superconductor, Nano Letters 19, 4934–4940 (2019).
  * Ni  _et al._ [2019] G. X. Ni, H. Wang, B.-Y. Jiang, L. X. Chen, Y. Du, Z. Y. Sun, M. D. Goldflam, A. J. Frenzel, X. M. Xie, M. M. Fogler, and D. N. Basov, Soliton superlattices in twisted hexagonal boron nitride, Nature Communications 10, 4360 (2019).
  * Yao and Yang [2015] H. Yao and F. Yang, Topological odd-parity superconductivity at type-ii two-dimensional van hove singularities, Phys. Rev. B 92, 035132 (2015).
  * Cobo  _et al._ [2016] S. Cobo, F. Ahn, I. Eremin, and A. Akbari, Anisotropic spin fluctuations in ${\mathrm{sr}}_{2}{\mathrm{ruo}}_{4}$: Role of spin-orbit coupling and induced strain, Phys. Rev. B 94, 224507 (2016).
  * Biderang  _et al._ [2019] M. Biderang, M.-H. Zare, and A. Akbari, Superconductivity of mixed parity and frequency in an anisotropic spin-orbit coupling, arXiv e-prints , arXiv:1906.08130 (2019), arXiv:1906.08130 .
  * Onari  _et al._ [2004] S. Onari, R. Arita, K. Kuroki, and H. Aoki, Phase diagram of the two-dimensional extended hubbard model: Phase transitions between different pairing symmetries when charge and spin fluctuations coexist, Phys. Rev. B 70, 094523 (2004).
  * Greco  _et al._ [2020] A. Greco, M. Bejas, and A. P. Schnyder, Ferromagnetic fluctuations in the rashba-hubbard model, Phys. Rev. B 101, 174420 (2020).
  * Mahan [2000] G. D. Mahan,  _Many-Particle Physics_ (Springer-Verlag New York Inc., 2000).
  * Wolf  _et al._ [2018] S. Wolf, T. L. Schmidt, and S. Rachel, Unconventional superconductivity in the extended hubbard model: Weak-coupling renormalization group, Phys. Rev. B 98, 174515 (2018).
  * Rømer  _et al._ [2015] A. T. Rømer, A. Kreisel, I. Eremin, M. A. Malakhov, T. A. Maier, P. J. Hirschfeld, and B. M. Andersen, Pairing symmetry of the one-band hubbard model in the paramagnetic weak-coupling limit: A numerical rpa study, Phys. Rev. B 92, 104505 (2015).
  * Lu  _et al._ [2018] C. Lu, L.-D. Zhang, X. Wu, F. Yang, and J. Hu, $d+id$ chiral superconductivity in a triangular lattice from trigonal bipyramidal complexes, Phys. Rev. B 97, 165110 (2018).
  * Wolf and Rachel [2020] S. Wolf and S. Rachel, Spin-orbit coupled superconductivity: Rashba-hubbard model on the square lattice, Phys. Rev. B 102, 174512 (2020).
  * Byers  _et al._ [1993] J. M. Byers, M. E. Flatté, and D. J. Scalapino, Influence of gap extrema on the tunneling conductance near an impurity in an anisotropic superconductor, Phys. Rev. Lett. 71, 3363–3366 (1993).
  * Akbari and Thalmeier [2013] A. Akbari and P. Thalmeier, Rashba spin-orbit coupling effects in quasiparticle interference of non-centrosymmetric superconductors, EPL (Europhysics Letters) 102, 57008 (2013).