# Reducing the X-ray radiation exposure frequency in cardio-angiography via
deep-learning based video interpolation

Xiao-Lei Yin1, Dong-Xue Liang 1, Lu Wang1, Jing Qiu1, Zhi-Yun Yang3 Jun-Hui
Xing 4 Jian-Zeng Dong 3 and Zhao-Yuan Ma 1  Dong-Xue Liang is with The Futrue
Laboratory, Tsinghua University, Beijing, 10084, China.(E-mail:
<EMAIL_ADDRESS>Xiao-Lei Yin, Lu Wang, Zhao-Yuan Ma and Jing Qiu
are with The Futrue Laboratory, Tsinghua University, Beijing, China. Zhi-Yun
Yang and Jian-Zeng Dong are with Center for Cardiology, Anzhen Hospital,
Capital Medical University, Beijing, China. Jun-Hui Xing is with The First
Affiliated Hospital of Zhengzhou University, Zhengzhou, China.

###### Abstract

Cardiac coronary angiography is a major technology to assist doctors during
cardiac interventional surgeries. Under the exposure of X-ray radiation,
doctors inject contrast agents through catheters to determine the position and
status of coronary vessels in real time. To get a coronary angiography video
with a high frame rate, the doctor needs to increase the exposure frequency
and intensity of the X-ray. This will inevitably increase the X-ray harm to
both patients and surgeons. In this work, we innovatively utilize a deep-
learning based video interpolation algorithm to interpolate coronary
angiography videos. Moreover, we establish a new coronary angiography image
dataset ,which contains 95,039 triplets images to retrain the video
interpolation network model. Using the retrained network we synthesize high
frame rate coronary angiography video from the low frame rate coronary
angiography video. The average peak signal to noise ratio(PSNR) of those
synthesized video frames reaches 34dB. Extensive experiment results
demonstrate the feasibility of using the video frame interpolation algorithm
to synthesize continuous and clear high frame rate coronary angiography video.
With the help of this technology, doctors can significantly reduce exposure
frequency and intensity of the X-ray during coronary angiography.

###### Index Terms:

Coronary angiography, contrast agent, video interpolation, deep learning

##  I Introduction

Cardiac coronary angiography is a major technology to assist doctors during
cardiac interventional surgeries. In Coronary angiography, a radiopaque
contrast agent is injected into a blood vessel and X-rays are taken to produce
detailed images of the blood vessel. Coronary angiography provides information
about the coronary arteries, which supply the heart with oxygen-rich blood.
During insertion, the doctor uses fluoroscopy (a continuous x-ray procedure)
to observe the progress of the catheter as it is threaded into place. Side
effects of radiopaque contrast agents include allergic reactions and kidney
damage. To get a coronary angiography video with a high frame rate, the doctor
needs to increase the exposure frequency and intensity of the X-ray. This will
inevitably increase the X-ray harm to patients and surgeons. The video frame
interpolation is a solution that can synthesize coronary angiography video
with a high frame rate without increasing the exposure frequency and intensity
of the X-ray.

Video frame interpolation is mainly to synthesize several frames in the middle
of two adjacent frames of the original video. Video frame interpolation is a
widely used method in the field of video processing, but never used in
coronary angiography videos. Video frame interpolation can be applied to
generate slow motion video, increase video frame rate, and frame recovery in
video streaming. The traditional video frame interpolation algorithm mainly
includes the following steps: bidirectional motion estimation, motion
interpolation, occlusion inference, and motion compensated frame
interpolation. However, these methods introduce various disturbances, such as
blur, ghost etc. In the past few years, deep-learning methods, especially deep
convolutional network methods[2], [3], [7], [12]-[14] have been effectively
applied and extended in the field of video frame interpolation. These methods
use neural networks to extract kernel functions and optical flow features, and
Apply adjacent frames to compose the middle frame. Although these methods have
been greatly improved, the effect is still not very good when dealing with
large object motion and occlusion.

The DAIN method[2] has good performance in the field of video frame
interpolation. In the method ,a variety of neural network structures are used
to extract image context features, optical flow features, depth map
features[4][5][6], and applies adaptive convolutional layers[3] to synthesize
new video frames. This method can deal with the situation of adapting to
strong motion and large object occlusion.

In this work, we apply video frame interpolation to coronary angiography video
interpolation to increase the frame rate and generate slow motion videos. The
video frame interpolation algorithm we choose is the DAIN method. Moreover, we
establish a new dataset containing 95,039 triplets where each triplet contains
3 consecutive coronary angiographic frames to retrain the network model of
DAIN.Then we utilize the trained network model to interpolate the coronary
angiography videos and analyze the results of interpolation frames. Finally,
we compare the performance of different frame interpolation algorithms in
different periods of the cardiac cycle.

We make the following contributions in this work:

1\. We innovatively apply deep learning-based video frame interpolation
algorithms to coronary angiography videos.

2\. We make a new coronary angiography dataset for video interpolation
algorithm.

##  II Related Work

Recently deep learning and specifically convolutional neural networks(CNNs)
have been successfully applied in computer vision areas, which inspired
various deep learning based frame interpolation methods. The early works on
CNN-based video frame interpolation[1] proposed an CNN architecture that takes
two input frames and directly estimate the intermediate frame. However, these
typical approaches often lead to blurred outputs.

The later methods, instead of directly computing the output pixels, mainly
focused on where to find the pixel from input frames and estimated spatially-
adaptive interpolation kernels to synthesize pixels from a large neighborhood.
In the AdaConv method[13], the adaptive filter is convolved with adjacent
input frame images to synthesize the data of the middle frame. But this method
can only deal with moving objects within a size of 41x41 pixels at most which
can not meet the situation of strenuous exercise and occlusion. A separable
neural network calculation method is applied in the SepConv[12] method, which
can greatly save the use of physical memory during the frame insertion process
and shorten the time used to train the network.

In the SuperSlomo method[7], a CNN is used for bidirectional motion estimate,
then a simplified motion interpolation method is applied, and at the end a
second CNN performs motion estimate refinement and occlusion reasoning. This
work achieved overwhelming quality when applied on videos taken at high frame-
rates. However, it seems that they do not aim at covering a wide range of
motions.

The adaptive deformable neural network[19] is defined in creative deformable
convolution. This method trains a weight kernel at each pixel value position,
and then each kernel weight corresponds to a pixel value with an offset. Any
positive number, so the bi-linear interpolation method is used to synthesize
the pixel value at any position. It performs better under intense movement and
occlusion

In the DAIN method[2], a pretrained model is used to extract the context
feature map[14], depth features map[8], [9], frame interpolation kernel, and
optical flow feature map[11] of the two adjacent frames, and the adaptive
convolution layer is used to synthesize any frame data between two adjacent
frames. The DAIN method can synthesize video frames at arbitrary positions and
has a good performance in dealing with the situation of strong motion and
occlusion of objects.

##  III Method

In this section, we rely on the DAIN method proposed by Bao.w et al.[2]. We
first introduce the DAIN method and describe the main architecture of the
network. And then we redefine the loss function and introduce the
implementation details.

Figure 1:  Architecture of DAIN model[2]. The model consists of the following
submodules: the flow estimation, depth estimation, context extraction, kernel
estimation, Depth-ware flow projection layer, adaptive warping layer, and
frame synthesis networks.

###  III-A Depth-Aware Video Interpolation

The main architecture of the DAIN model is shown in Fig.1. Given two input
images $I_{0}$ and $I_{1}$ and a time $t\in(0,1)$, the goal is to predict the
interpolate image $\hat{I}_{t}$ at time $T=t$. The DAIN model consists of the
following submodules: the flow estimation, depth estimation, context
extraction, kernel estimation, and frame synthesis networks. The DAIN
algorithm uses the proposed depth-aware flow projection layer to obtain
intermediate flows and then warps the input frames,depth maps, and contextual
features within the adaptive warping layer. Finally, the frame synthesis
network generates the output frame with residual learning.

The DAIN model utilizes the pretrained PWC-Net[11] as the optical flow
estimation network. $F_{0\to 1}$ and $F_{1\to 0}$ are the estimated optical
flows by the flow estimation network . $F_{0\to 1}$ and $F_{1\to 0}$ denote
the optical flow from $I_{0}$ to $I_{1}$ and$I_{1}$ to $I_{0}$ ,respectively.
The depth estimation network use the hourglass architecture[9], [16]. $D_{0}$
and $D_{1}$ are the estimated depth map of $I_{0}$ and $I_{1}$. The DAIN
method contains a depth-aware flow projection[2] which combine the optical
flow and the depth map to get the projected flow at time $T=t$ . The projected
flow $F_{t\to 0}$ can be obtained from the flow $F_{0\to 1}$ and the depth map
$D_{0}$. The projected flow $F_{t\to 1}$ can be obtained from the flow
$F_{1\to 0}$ and the depth map $D_{1}$.

The projected flow $F_{t\to 1}$ is defined by:

$F_{t\to
1}(x)=-(1-t)\cdot{\frac{\sum\limits_{y\in\mathcal{S}(x)}\frac{1}{{D}_{1}(y)}\cdot
F_{1\to 0}(y)}{\sum\limits_{y\in\mathcal{S}(x)}\frac{1}{{D}_{1}(y)}}}$ (1)
$\mathcal{S}(x)=\left\\{y:round(y+(1-t)F_{1\to 0}(y))=x,\forall
y\in[H,W]\right\\}$ (2)

Where $\mathcal{S}(x)$ indicates the set of pixels that pass through the
position $x$ at time $t$. The projected flow $F_{t\to 0}$ is defined by:

$F_{t\to
0}(x)=-t\cdot{\frac{\sum\limits_{y\in\mathcal{S}(x)}\frac{1}{{D}_{0}(y)}\cdot
F_{0\to 1}(y)}{\sum\limits_{y\in\mathcal{S}(x)}\frac{1}{{D}_{0}(y)}}}$ (3)

$\mathcal{S}(x)=\left\\{y:round(y+tF_{0\to 1}(y))=x,\forall
y\in[H,W]\right\\}$ (4)

The context extraction model consists a pretrained ResNet[17]. The contextual
features $C_{0}$ and $C_{1}$ are extracted from the input frames $I_{0}$ and
$I_{1}$.The interpolation kernels $K_{0}$ and $K_{1}$ are estimated by a U-Net
architecture[15] network. With the interpolation kernels ($K_{0}$ and $K_{1}$)
and interpolated flows($F_{t\to 0}$ and $F_{t\to 1}$) generated from the
depth-aware flow projection layer, using the adaptive warping layer[3] to warp
the input frames( $I_{0}$ and $I_{1}$), depth maps($D_{0}$ and
$D_{1}$),contextual features($C_{0}$ and $C_{1}$). To generate the output
frame$\widehat{I}_{t}$, the algorithm utilizes a frame synthesis network,
which consists of 3 residual blocks. The algorithm concatenates projected
flows( $\widehat{F}_{t\to 0}$ and $\widehat{F}_{t\to 1}$ ), and interpolation
kernels($K_{0}$ and $K_{1}$),the warped input frames($\widehat{I}_{0}$ and
$\widehat{I}_{1}$), warped depth maps($\widehat{D}_{0}$ and
$\widehat{D}_{1}$), warped contextual features($\widehat{C}_{0}$ and
$\widehat{C}_{1}$) as the input to the frame synthesis network. In addition,
the algorithm linearly blend the two warped frames($\widehat{I}_{0}$ and
$\widehat{I}_{1}$) and enforce the network to predict the residuals between
the ground-truth frame and the blended frame.

###  III-B Implementation Details

Loss Function. We denote the synthesized frame by$\hat{I}_{t}$ and the ground-
truth frame by $I_{t}$. We retrain the DIAN model by optimizing the following
loss function:

$L={\lambda_{1}l_{1}+\lambda_{2}l_{2}}$ (5)

where $\lambda_{1}$ = 0.95 and $\lambda_{2}$ = 0.05.

The loss $l_{1}$ models how good the reconstruction fo the intermediate frame
is:

$l_{1}=\sum\limits_{x}\rho||\hat{I}_{t}(x)-I_{t}(x)||_{1}$ (6)

where $\rho(x)=\sqrt{x^{2}+\epsilon^{2}}$is the Charbonnier Function. We set
the constant $\epsilon$ to $1e-4$.

The purpose of coronary angiography video frame interpolation is to obtain a
clearer coronary artery structure. So we define the structure loss $l_{2}$ as

$\ l_{2}=\sum\limits_{x}\rho||\phi{\hat{I}}_{t}(x)-\phi I_{t}(x)||_{1}$ (7)

where $\phi$ denote the conv4_x features of a pretrained ResNet[17].

Training Dataset. In order to train a better network model conforming to
coronary angiography scene we create a new training dataset. The dataset
contains 950,39 triplets of consecutive frames with a resolution of 480x360.
The triplets are extracted from 31 coronary angiography videos ，the total
duration of which is about 25 hours. The videos are produced by recording
screen images during cardiac intervention. We processed these videos and
selected suitable coronary angiography fragments as the data set. The image
frame rate of coronary angiography is 7-15fps, but the frame rate used for
recording screen images is 25fps, so we preprocess the recorded videos to
remove duplicate video frames. Coronary angiography images contain some
patient’s information, so we occlude and remove these information to ensure
that the patient’s privacy is not leaked. Also we augument the training data
by horizontal and vertical flipping.

Training strategy. Every triplet in the dataset contains 3 consecutive video
frames, which are the previous frame, the middle frame, and the next frame. we
train the network to predict the middle frame that serves as ground truth. We
use a pre-trained model as the initial model for training. The learning rate
is initially 0.0001 and decays half every 20 epochs. The batch size is 4 and
the the network is trained for 100 epochs. We train the network on an NVIDIA
2080Ti GPU card,which takes about 30 hours.

##  IV Experiments and Results

In this section, we first introduce the testing datasets. Second, we compare
the output of the pretrained model with the retrained DAIN model. Then we
compare our results with stat-of-the-art frame interpolation approaches and
analyze the performance of different frame interpolation methods in different
cardiac cycles. Finally we discuss how to further improve the current results.

###  IV-A Testing Dataset

We randomly selected 1,000 triplets from the coronary angiography dataset as
one test dataset which is represented by D1000. At the same time, we extracted
four coronary angiography video clips, each of which contained 30-45
continuously frames, and each video clip includes the process of injection,
diffusion, and disappearance of the contrast agent. We use VC1, VC2, VC3, VC4
to represent these four video clips. For these four clips, we use odd frames
to predict even frame, and then use even frames to predict odd frame, so we
can compare all the predicted frames with the original frames. We utilize the
PSNR and SSIM values as the evaluation criteria for the two-frame image error.

TABLE I:  The experiments results of different datasets. Test data | D1000 | VC1 | VC2
---|---|---|---
| SSIM,PSNR | SSIM,PSNR | SSIM,PSNR
DAIN-Pre | 0.961 40.142 | 0.902 33.839 | 0.900 34.368
DAIN-Re | 0.962 40.523 | 0.904 34.068 | 0.902 34.754

Original frame $I_{0}$

Original frame $I_{1}$

$I^{Pre}_{t=0.5}$,PSNR=34.91

Ground truth

$I^{Re}_{t=0.5}$,PSNR=35.28

Figure 2:  Comparison of the output of the retrained model and the pretrained
model. The $I^{Pre}_{t=0.5}$ is the predicted video frame of pretrained model
and the $I^{Re}_{t=0.5}$ is the predicted video frame of retrained model
$0$$5$$10$$15$$20$$25$$30$$30$$32$$34$$36$$38$$40$$42$Index of
framesPSNR【[dB]DAIN-PreDAIN-Re Figure 3:  PSNR of interpolation frame changes
with frame index, the test data is VC4

###  IV-B Results

We use the testing dataset to compare the results of the pretrained DAIN model
with the retrained model. Fig.2 shows the inputs of two original frames, the
original middle frame as ground truth, the interpolation midldle frame of the
pretrained model, and interpolation midldle frame of the retrained model. The
PSNR of $I^{Re}_{t=0.5}$ is 0.3dB higher than $I^{Pre}_{t=0.5}$. The
comparison of the average values of PSNR and SSIM for D1000 and two video
clips is recorded in TableI. It shows the retrained model has higher PSNR and
SSIM than the pretrained model. Fig.3 shows the PSNR changing with frame index
of one coronary angiography video clip. The retrained model has higher PSNR on
every frame. From Fig.3 we can see that the PSNR value has a similar periodic
change, which is related to the cardiac cycle. We will analyze the reasons for
this phenomenon in the subsequent results.

Ground truth

SepConv-l1

SepConv-lf

SuperSlomo

DAIN

Figure 4: The interpolation Frame of DAIN , SepConv-l1 ,SepConv-lf ,SuperSlomo. It mainly shows the performance of several methods in the case of occlusion and large object movement. As can be seen from the figure, the DAIN method performs better in both aspects.  TABLE II:  The average PSNR value and SSIM value of several frame interpolation algorithms. Test data | VC1 | VC2 | VC3 | VC4 | D1000
---|---|---|---|---|---
| SSIM , PSNR | SSIM , PSNR | SSIM , PSNR | SSIM , PSNR | SSIM , PSNR
SepConv-l1 | 0.901 33.657 | 0.901 34.672 | 0.904 33.649 | 0.902 33.745 | 0.962 39.823
SepConv-lf | 0.878 32.768 | 0.883 34.039 | 0.883 33.003 | 0.883 32.872 | 0.956 39.741
SuperSlomo | 0.901 33.838 | 0.898 34.468 | 0.902 33.768 | 0.901 33.704 | 0.958 38.868
DAIN | 0.904 34.068 | 0.902 34.754 | 0.908 33.977 | 0.908 34.205 | 0.962 40.523

$0$$5$$10$$15$$20$$25$$30$$35$$28$$30$$32$$34$$36$$38$$40$$42$$44$$46$$48$Index
of framesPSNR【[dB]DAINSuperSlomoSepConv-l1SepConv-lf

$0$$5$$10$$15$$20$$25$$30$$35$$40$$45$$28$$30$$32$$34$$36$$38$$40$$42$$44$$46$$48$Index
of framesPSNR [dB]DAINSuperSlomoSepConv-l1SepConv-lf

(a)

(b)

Figure 5:  The PSNR curves of several frame interpolation algorithms on two
continuous coronary angiography clips.(a) shows the results of VC2 and (b)
shows the results of VC3. The vertical red line indicates the beginning of
systole and the vertical black line indicates the beginning of diastole.

###  IV-C Comparison and Discussion

We compare our results with several stat-of-the-art algorithms including
SepConv-l1[12], SepConv-lf[12], SuperSlomo[7]. In order to exclude the
influence of the training dataset, we use the coronary angiography dataset to
retrain these several network models respectively. Then we use these retrained
models to process the testing datasets.

In Table II, it shows the comparisons of different methods on the coronary
angiography testing datasets. The tables summarizes the average PSNR and SSIM
of these experiments. The average PSNR of DAIN method is 0.2dB-1.6dB higher
than the other methods.

In Fig.4, it shows the interpolation frames of several video interpolation
methods in the case of occlusion and large object movement. As shown in the
figure, the DAIN method performs well in both aspects. The image edges of the
interpolated frames are clearer and sharper.

In Fig.5, we plot the PSNR curves of several frame interpolation algorithms on
two coronary angiography video clips. First, the overall PSNR curve of DAIN is
higher than other methods. Second, all the PSNR curves has a periodic change.
When the systolic and diastolic phases are switched in the cardiac cycle,
lower PSNR values appear. We analyze the reason for this result is that when
the phase of the cardiac cycle is switched, the optical flow characteristics
between two adjacent frames are not obvious enough, and the possible motion
vectors of the coronary arteries cannot be accurately predicted. At the same
time, we can see that the PSNR value decreases with the index of the frame.
The reason for this problem is that the vascular motion area of the coronary
angiography image continues to expand with the diffusion of the contrast agent
in the coronary arteries. As shown in Figure 5 (b), after the contrast agent
has diffused to the end of the coronary artery, the PSNR no longer changes
significantly between the 35th and 45th frames.

In order to further improve the quality of coronary angiography interpolation
frame, we can use the trained coronary artery segmentation [18]network model
to extract the coronary artery segmentation results from the synthetic frame
and the original frame, respectively. We take the error of coronary artery
segmentation in two frames as part of the loss function. This method may also
improve the problem that the interpolation frames have lower PSNR values when
the cardiac cycle phase switched.

##  V Conclusion

In this paper, we innovatively apply the deep learning-based video frame
interpolation algorithm to coronary angiography. We establish a new coronary
angiography dataset to retrain the network of DAIN method. The retrained model
has a better performance in the application scenarios of coronary angiography.
Moreover, we retrain several other deep learning-based algorithms and compare
the results of these frame interpolation algorithms. The retained DAIN model
has better performance than other methods in the case of occlusion and large
object movement. Extensive experiment results demonstrate the feasibility of
using the video frame interpolation algorithm to synthesize continuous and
clear high frame rate coronary angiography video. With the help of this
technology, doctors can significantly reduce exposure frequency and intensity
of the X-ray during coronary angiography.

## References

  * [1] G. Long et al., ”Learning image matching by simply watching video”, Proc. Eur. Conf. Comput. Vis, pp. 434–450, 2016.
  * [2] W. Bao et al., ”Depth-Aware Video Frame Interpolation”, Proc. IEEE Conf. Comput. Vis. Pattern Recog., pp. 3703-3712, 2019.
  * [3] W. Bao et al., ”MEMC-Net: Motion Estimation and Motion Compensation Driven Neural Network for Video Interpolation and Enhancement”, IEEE TPAMI, Sep. 2019.
  * [4] W. Chen et al., ”Single-image depth perception in the wild”, Proc. Adv. Neural Inform. Process. Syst., pp. 730-738, 2016.
  * [5] D. Eigen et al., ”Depth map prediction from a single image using a multi-scale deep network”, Proc. Adv. Neural Inform. Process. Syst., pp. 2366-2374, 2014. 2
  * [6] H. Fu et al., ”Deep ordinal regression network for monocular depth estimation”, Proc. IEEE Conf. Comput. Vis. Pattern Recog., pp.2002-2011, 2018.
  * [7] H. Jiang et al., ”Super SloMo: High Quality Estimation of Multiple Intermediate Frames for Video Interpolation”, Proc. IEEE Conf. Comput. Vis. Pattern Recog., pp. 9000-9008, 2018.
  * [8] Y. Kuznietsov et al., ”Semi-supervised deep learning for monocular depth map prediction”, Proc. IEEE Conf. Comput. Vis. Pattern Recog., pp. 2215-2223, 2017.
  * [9] Z. Li and N. Snavely, ”Megadepth: Learning single-view depth prediction from internet photos”, Proc. IEEE Conf. Comput. Vis. Pattern Recog., pp. 2041-2050, 2018.
  * [10] Y. Zou et al., ”DF-Net: Unsupervised Joint Learning of Depth and Flow using Cross-Task Consistency”, Proc. Eur. Conf. Comput. Vis., pp.38-55, 2018.
  * [11] D. Sun et al., ”Pwc-net: Cnns for optical flow using pyramid, warping, and cost volume”, Proc. IEEE Conf. Comput. Vis. Pattern Recog., pp. 8934-8943, 2018.
  * [12] S. Niklaus et al., ”Video frame interpolation via adaptive separable convolution”, Proc. Int. Conf. Comput. Vis., 2017.
  * [13] S. Niklaus et al., ”Video frame interpolation via adaptive convolution”, Proc. IEEE Conf. Comput. Vis. Pattern Recog., pp. 2270-2279, 2017.
  * [14] S. Niklaus and F. Liu., ”Context-aware synthesis for video frame interpolation”, Proc. IEEE Conf. Comput. Vis. Pattern Recog., pp. 1701-1710, 2018.
  * [15] O. Ronneberger et al., ”U-Net: Convolutional networks for biomedical image segmentation”, Proc. Int. Conf. Med. Image Comput. Comput.-Assisted Intervention, pp. 234-241, 2015.
  * [16] W. Chen et al., ”Single-image depth perception in the wild”,Proc. Adv. Neural Inform. Process. Syst., pp.730-738, 2016.
  * [17] K. He et al., ”Deep residual learning for image recognition”, Proc. IEEE Conf. Comput. Vis. Pattern Recog., pp. 770-778, 2016.
  * [18] K. Jo et al., ”Segmentation of the Main Vessel of the Left Anterior Descending Artery Using Selective Feature Mapping in Coronary Angiography”, IEEE Access, vol. 7, pp. 919-930, 2019.
  * [19] H. Lee et al., ”Learning Spatial Transform for Video Frame Interpolation”, arXiv:1907.10244, 2019.

> >


# Type-II topological phase transitions of topological skyrmion phases

Reyhan Ay  Max Planck Institute for Chemical Physics of Solids, Nöthnitzer
Straße 40, 01187 Dresden, Germany  Max Planck Institute for the Physics of
Complex Systems, Nöthnitzer Straße 38, 01187 Dresden, Germany  Izmir Institute
of Technology, Gülbahçe Kampüsü, 35430 Urla Izmir, Türkiye  Joe H. Winter  Max
Planck Institute for Chemical Physics of Solids, Nöthnitzer Straße 40, 01187
Dresden, Germany  Max Planck Institute for the Physics of Complex Systems,
Nöthnitzer Straße 38, 01187 Dresden, Germany  SUPA, School of Physics and
Astronomy, University of St. Andrews, North Haugh, St. Andrews KY16 9SS, UK
A. M. Cook  Max Planck Institute for Chemical Physics of Solids, Nöthnitzer
Straße 40, 01187 Dresden, Germany  Max Planck Institute for the Physics of
Complex Systems, Nöthnitzer Straße 38, 01187 Dresden, Germany

###### Abstract

We present minimal toy models for topological skyrmion phases of matter, which
generically realize type-II topological phase transitions in effectively non-
interacting systems, those which occur without closing of the minimum direct
bulk energy gap. We study the bulk-boundary correspondence in detail to show
that a non-trivial skyrmion number yields a rich bulk-boundary correspondence.
We observe gapless edge states, which are robust against disorder, due to non-
trivial skyrmion number. Edge states corresponds to bands, which do not
traverse the bulk gap, instead yielding gaplessness due to their overlap in
energy and exponential localization on opposite edges of the system. These
gapless boundary modes can occur for total Chern number zero, and furthermore
correspond to rich real-space spin textures with strong polarization of spin
along the real-space edge. By introducing toy models generically exhibiting
type-II topological phase transitions and characterizing the bulk-boundary
correspondence due to non-trivial skyrmion number in these models, we lay the
groundwork for understanding consequences of the quantum skyrmion Hall effect.

The quantum Hall effect (QHE), in which a two-dimensional electron gas
subjected to an out-of-plane external magnetic field exhibits Hall
conductivity quantized to rational numbers in units of $e^{2}/h$ [1, 2],
serves as the foundation for much of our understanding of topological phases
of matter [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. In particular, it
is the foundation for those topological states defined by mappings to the
space of projectors onto occupied states as in the ten-fold way classification
scheme [17, 18, 19, 20, 21, 22, 23, 24, 25]. Recently-introduced topological
skyrmion phases (TSPs) of matter [26], however, instead derive from a much
larger set of mappings to the space of myriad observable expectation values,
$\langle\mathcal{O}\rangle$. Such mappings can also be divided into
topologically-distinct sectors and understood as lattice counterparts of a
generalization of the QHE to the quantum skyrmion Hall effect (QSkHE) [27].
Ultimately, this relates to the generalization of the concept of a particle:
in the QSkHE, point charges of the QHE generalize to point-like quantum
skyrmions forming in myriad observable fields, with corresponding
generalizations of gauge fields, and $p$-dimensional charges are more
generally topological textures in underlying fields. This physics serves as
the link between topological states descending from the QHE, and extended
magnetic skyrmions [28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41].

Figure 1: Phase diagrams characterizing bulk topology of Hamiltonian Eq. 1 for
half-filling, including a) total Chern number $\mathcal{C}$, b) skyrmion
number $\mathcal{Q}$, c) minimum direct bulk energy gap $\Delta
E_{\mathrm{min}}$, and d) minimum ground state spin magnitude
$|\boldsymbol{S}|_{\mathrm{min}}$, with each of these four quantities plotted
vs. atomic spin-orbit coupling (SOC) strength $\lambda$ and mass term $M$.
Other parameters are fixed at $\beta=0.5$ and $\Delta=0.6$.

A variety of TSPs have been identified [26, 27, 42, 43], but the significance
of these implications strongly motivates efforts to better understand this
physics. A fundamental signature of TSPs is the type-II topological phase
transition [26], in which a topological invariant changes in effectively non-
interacting systems, while maintaining fixed occupancy of the bands and while
respecting the symmetries protecting the topological phase, without the
closing of the minimum direct bulk energy gap. This occurs, for instance, when
topological skyrmion phases are realized for mappings to the spin expectation
value of occupied states and spin is not conserved due to non-negligible
atomic spin-orbit coupling [26]. In this work, we therefore make a key
contribution to understanding this broader set of topological phases
associated with the QSkHE, by investigating the type-II topological phase
transition in minimal models. We show that type-II topological phase
transitions are generic to these models, and also characterize a bulk-boundary
correspondence specific to topological skyrmion phases, distinguished by
topologically-robust, gapless boundary modes with distinctive spin textures.

Minimal models for type-II topological phase transitions—We construct toy
models generically realizing type-II topological phase transitions inspired
very directly by tight-binding Bogoliubov de Gennes Hamiltonians for
superconducting Sr2RuO4 in which topological skyrmion phases were first
discovered [26, 44, 45], which we also provide in the Supplementary Materials
(Supp. Mat.), Section I, “Tight-binding model for mirror subsectors of
superconducting Sr2RuO4 Bogoliubov de Gennes Hamiltonian and associated spin
operators”. Our results are therefore very relevant to transition metal oxide
superconductors.

We consider minimal Bloch Hamiltonians realizing topological skyrmion phases
of the following form, with three essential terms,

$\mathcal{H}(\boldsymbol{k})=\boldsymbol{d}(\boldsymbol{k})\cdot\boldsymbol{\tilde{S}}+\mathcal{H}_{\mathrm{pair}}(\boldsymbol{k})+\mathcal{H}_{\mathrm{SOC}},$
(1)

and basis vector

$\Psi^{\dagger}_{\boldsymbol{k}}=\left(c_{\boldsymbol{k},xy,\uparrow},c_{\boldsymbol{k},yz,\downarrow},c_{\boldsymbol{k},xz,\downarrow},c^{\dagger}_{\boldsymbol{-k},xy,\uparrow},c^{\dagger}_{\boldsymbol{-k},yz,\downarrow},c^{\dagger}_{\boldsymbol{-k},xz,\downarrow}\right),$
(2)

where $c_{\boldsymbol{k},\ell,\sigma}$ annihilates a fermion with momentum
$\boldsymbol{k}$ in orbital $\ell$ with spin $\sigma$, and $\\{xy,yz,xz\\}$
correspond to the t2g orbital degree of freedom (dof) and
$\\{\uparrow,\downarrow\\}$ correspond to a spin $1/2$ dof. We note that the
unconventional ordering of annihilation and creation operators results from
rotating the basis of the Sr2RuO4 tight-binding model to that for which the
mirror operator taking $z\rightarrow-z$ is diagonal [44, 26]. Here,
$\boldsymbol{d}(\boldsymbol{k})\cdot\boldsymbol{\tilde{S}}$ is the term
required to produce a non-trivial skyrmionic texture in the ground state spin
expectation value, where
$\boldsymbol{\tilde{S}}=\mathrm{diag}\left(\boldsymbol{S},-\boldsymbol{S}^{*}\right)$,
and $\boldsymbol{S}$ is the vector of spin operators, $\mathcal{S}_{x}$,
$\mathcal{S}_{y}$, and $\mathcal{S}_{z}$, for the particle sector. This normal
state spin representation is provided in the Supp. Mat., Section I, “Tight-
binding model for mirror subsector I of superconducting Sr2RuO4 Bogoliubov de
Gennes Hamiltonian”. The second term,
$\mathcal{H}_{\mathrm{pair}}(\boldsymbol{k})$, is a pairing term required to
couple the generalized particle-hole conjugate sectors (for a Bogoliubov de
Gennes Bloch Hamiltonian, $\mathcal{H}_{\mathrm{pair}}(\boldsymbol{k})$ is the
superconducting pairing term). Finally, $\mathcal{H}_{\mathrm{SOC}}$ is the
atomic spin-orbit coupling (SOC) term. We find that the first two terms yield
only the type-I topological phase transition, and the final SOC term is
required to produce the type-II transition as expected from past work.

We focus here on study of six-band models, for which the type-II topological
phase transition has already been observed in work introducing the topological
skyrmion phases, for a tight-binding model previously used to study Sr2RuO4
with spin-triplet pairing [45, 44, 26]. To construct a concrete toy model, we
use the spin operators, atomic spin-orbit coupling term, and pairing term of
mirror subsector I of the previously-studied Bogoliubov de Gennes Bloch
Hamiltonian for Sr2RuO4 with spin-triplet pairing. The form of this mirror
subsector Bloch Hamiltonian is reviewed in the Supp. Mat., Section I, “Tight-
binding model for mirror subsector I of superconducting Sr2RuO4 Bogoliubov de
Gennes Hamiltonian”. The spin-orbit coupling term of the normal state
$h_{SOC}$ and the superconducting gap function of the pairing term
$\Delta_{1}(\boldsymbol{k})$ take the following forms, respectively,

$\displaystyle h_{SOC}=\begin{pmatrix}0&-i\lambda&-\lambda\\\
i\lambda&0&i\lambda\\\ -\lambda&-i\lambda&0\end{pmatrix},\hskip
5.69054pt\Delta_{1}(\boldsymbol{k})=\mathrm{diag}\left(\delta_{-},\delta_{+},\delta_{+}\right),$
(3)

where $\lambda$ is atomic spin-orbit coupling strength,
$\delta_{\mp}=\left(\Delta_{0}/2\right)\left[i\sin(k_{x})\mp\sin(k_{y})\right]$,
and $\Delta_{0}$ is superconducting pairing strength.

We choose a $\boldsymbol{d}(\boldsymbol{k})=\langle
d_{x}(\boldsymbol{k}),d_{y}(\boldsymbol{k}),d_{z}(\boldsymbol{k})\rangle$
vector previously used to describe one of the simplest two-band models for a
Chern insulator on a square lattice, known as the QWZ model [46],
$d_{x}(\boldsymbol{k})=\sin(k_{x})$, $d_{y}(\boldsymbol{k})=\sin(k_{y})$,
$d_{z}(\boldsymbol{k})=\beta(2+\mathcal{M}-\cos(k_{x})-\cos(k_{x}))$.

Bulk topology— For the chosen model given by Eq. 1, we first compute bulk
phase diagrams as a function of the mass parameter $\mathcal{M}$ of the QWZ
model and the atomic spin-orbit coupling strength $\lambda$ in parallel to
earlier work [26]. We assume half-filling of the six bands and compute the
total Chern number $\mathcal{C}$, the skyrmion number $\mathcal{Q}$, the
minimum direct bulk energy gap $\Delta E_{\mathrm{min}}$, and the minimum spin
magnitude over the Brillouin zone, $|\boldsymbol{S}|_{\mathrm{min}}$, as shown
in Fig. 1 a), b), c) and d), respectively. Additional characterization of bulk
topology for a second choice of $\boldsymbol{d}(\boldsymbol{k})$, previously
used by Sticlet  _et al._ [47] to realize two-band Chern insulator phases with
Chern number $\mathcal{C}=\pm 2$ for the lower band, is detailed in the Supp.
Mat., Section II, “Additional characterization of bulk topology for a second
choice of $\boldsymbol{d}(\boldsymbol{k})$ vector”.

The total Chern number $\mathcal{C}$ phase diagram, Fig. 1 a), exhibits a rich
variety of non-trivial regions, including a central set of four rhombus-like
regions formed by overlapping stripe-like regions. The total Chern number is
trivial for large $|\lambda|$ and $\mathcal{M}$ near $-2$. In constrast, the
skyrmion number $\mathcal{Q}$ phase diagram, Fig. 1 b) is simpler, with two
non-trivial regions of $\mathcal{Q}=\pm 1$, respectively. Considering these
two phase diagrams together, we find a variety of regions with $\mathcal{C}$
non-zero and $\mathcal{Q}$ zero and vice versa.

Examining the minimum direct bulk energy gap and minimum spin magnitude over
the Brillouin zone shown in Fig. 1 c) and d), respectively, we find a great
variety of both type-I ($\Delta E_{\mathrm{min}}$ goes to zero) and type-II
($\Delta E_{\mathrm{min}}$ remains finite while
$|\boldsymbol{S}|_{\mathrm{min}}$ goes to zero) topological phase transitions
[26]. We find such variety in values of the bulk topological invariants and in
the topological phase transitions is generic for models of the form Eq. 1,
making them invaluable for understanding the interplay between the total Chern
number $\mathcal{C}$ and skyrmion number $\mathcal{Q}$.

Figure 2: Slab spectra for Hamiltonian Eq. 1 with open boundary conditions
(OBCs) in the $\hat{x}$-direction and periodic boundary conditions (PBCs) in
the $\hat{y}$-direction. Overlapping edge states for non-trivial skyrmion
number $\mathcal{Q}$ are present as shown in a) and b), which become critical
as shown in c) during the type-II topological phase transition, before finally
pulling apart fully yielding a gap as shown in d). Parameter values are a)
$\mathcal{M}=-2.7$ and $\lambda=0.8$, b) $\mathcal{M}=-2.6$ and $\lambda=0.8$,
c) $\mathcal{M}=-2.558$ and $\lambda=0.8$, d) $\mathcal{M}=-1.9$ and
$\lambda=0.8$.

Bulk-boundary correspondence—We now characterize bulk-boundary correspondence
for the specific realization of Eq. 1 characterized in the bulk in Fig. 1,
focusing on better understanding the bulk-boundary correspondence associated
with changing $\mathcal{Q}$. To do so, we first consider four points in phase
space along a cut through the phase diagram Fig. 1 labeled ‘Transition A’, for
which $\mathcal{C}=0$ but $\mathcal{Q}$ changes from -1 to +1, with three
other type-II topological phase transitions B, C and D also shown in the Supp.
Mat., Section III, “Additional slab spectra through type-II topological phase
transitions”.

We plot the slab spectrum of the Hamiltonian for open boundary conditions
(OBCs) in the $\hat{x}$-direction and periodic boundary conditions (PBCs) in
the $\hat{y}$-direction along Transition A in Fig. 2. Within the
$\mathcal{Q}\neq 0$ region, we observe gaplessness at the edge due to
overlapping of two bands in the gap, exponentially-localized on opposite
edges. Edge states extend from the bulk valence (conduction) to bulk valence
(conduction) bands, similarly to edge states of TSPs in three-band models
[27]. During the type-II topological phase transition, the overlap of these
in-gap states shrinks and is finally lost. Gaplessness is lost in the vicinity
of the point at which $|\boldsymbol{S}_{\mathrm{min}}|=0$: the separation in
phase space between the point at which gaplessness is lost and the point at
which $\mathcal{Q}$ changes in the bulk tends to increase with increasing
$|\lambda|$, as shown in the Supp. Mat., Section IV, “Effect of atomic spin-
orbit coupling strength on transition A”. This is a remarkable result in
combination with the robustness of these gapless edge states against disorder
presented in the next section: although we gain considerable understanding of
the spin topology through observable-enriched entanglement [48] and the view
of non-trivial $\mathcal{Q}$ yielding an effective Chern insulator of the spin
sector upon tracing out non-spin degrees of freedom, the entanglement due to
non-negligible atomic spin-orbit coupling enriches the physics considerably
and requires further investigation. These results may correspond to one of the
more general scenarios of the QSkHE, corresponding to topological transport of
quantum skyrmion-like textures generalizing the point charges of the QHE,
which are not well-approximated by coarse-graining to point charges.

While gapless edge states associated with non-trivial $\mathcal{Q}$ are
present in Fig. 2 a) and lost over type-II topological phase transition A as
shown in Fig. 2 a) to d), gapless edge states are not observed over this cut
through the phase diagram for OBCs instead in the $\hat{y}$-direction and PBCs
in the $\hat{x}$-direction. While Transitions B and C are similarly
anisotropic, Transition D is more isotropic, with gapless edge states observed
for each set of boundary conditions in the vicinity of Transition D. These
additional slab spectra are shown in the Supp. Mat., Section III, “Additional
slab spectra through type-II topological phase transitions”.

Figure 3:  Further characterization of Transition A: a) Slab spectrum for
Hamiltonian Eq. 1 with parameters $\mathcal{M}=-2.8$ and $\lambda=0.8$.
$x$-component, $y$-component, and $z$-component of spin expectation value,
$\langle S_{x}(k_{y})\rangle$, $\langle S_{y}(k_{y})\rangle$, and $\langle
S_{z}(k_{y})\rangle$ for the two overlapping edge states in a) are shown in
b), c) and d), respectively. Filled (unfilled) symbols correspond to occupied
states at energy $E<0$ (unoccupied states at energy $E>0$). Color bar shows
localization of edge states as difference in probability density of a given
edge state at the left edge, $|\psi_{L}|^{2}$, and right edge,
$|\psi_{R}|^{2}$.

To further characterize the bulk-boundary correspondence, we also compute the
spin texture for the gapless edge states, as shown in Fig. 3. Despite rich
physics of the bulk Hamiltonian including non-negligible atomic spin-orbit
coupling, the gapless edge states exhibit very strong polarization of spin
along the edge. The polarization is chiral, as shown in Fig. 3 c): There is a
sign difference in $\langle S_{y}(k_{y})\rangle$ between the two edges.

Edge state spin textures for each of Transition A to D are shown in the Supp.
Mat. Section V, “Additional edge state spin textures through type-II
topological phase transitions”. Edge state spin polarization along sample
edges weakens with loss of gaplessness but remains pronounced.

Figure 4: a) Spectrum vs. eigenvalue index $i$ for system of size
$N_{x}=N_{y}=40$ and parameter set $\mathcal{M}=-2.6$ and $\lambda=0.8$ (along
Transition A cut) for different values of disorder strength $\rho$, for system
with OBCs in each of the $\hat{x}$\- and $\hat{y}$-directions except for PBC-x
(OBC in $\hat{y}$-direction, PBCs in $\hat{x}$-direction), or PBC-y (OBC in
$\hat{x}$-direction, PBC in $\hat{y}$-direction). Inset highlights low-energy
spectrum. b) Disorder-averaged probability density for eigenstate number i =
4800 from a) for $\rho=0.20$, with 100 disorder realizations. c) Log of
disorder-averaged probability density for eigenstate with eigenvalue -0.00121
and i = 4800 vs. position in the $\hat{x}$-direction for $y=20$, shown for
different values of $\rho$, with system size and parameter set as in a). d)
Disorder-averaged spin texture for eigenstate with eigenvalue -0.00121 and
$i=4800$ over real-space for $\rho=0.20$ as in b), with system size and
parameter set as in a).

Robustness of gapless boundary modes against disorder—We now study the
robustness of the gapless edge states due to non-zero $\mathcal{Q}$ against
disorder, by Fourier-transforming the Bloch Hamiltonian and studying the
system with OBCs in each of the $\hat{x}$\- and $\hat{y}$-directions, along
Transition A studied in Fig. 2 and 3. We introduce a spatially-varying on-site
potential term that is $\mathcal{C}^{\prime}$-invariant at each site $(x,y)$
$\mathcal{H}_{\mathrm{dis}}(x,y)=\rho\varepsilon(x,y)\mathrm{diag}\left(\mathbf{I}_{3},-\mathbf{I}_{3}\right)$,
where $\rho$ is disorder strength, $\varepsilon(x,y)$ is chosen from a uniform
random distribution over the interval $\left[-1,1\right]$, and
$\mathbf{I}_{3}$ is the $3\times 3$ identity matrix. The spectrum of the clean
system for PBCs and OBCs, as well as the disorder-averaged spectrum for OBCs,
is shown in Fig. 4 a). Eigenenergy $E$ increases steeply with index $i$ for
low-energy states for $\rho=0$, corresponding to edge states, and this slope
persists to $\rho=0.2$. At $\rho=0.4$, the slope near $E=0$ is more gradual
and bulk-like, indicating collapse of the bulk energy gap. For $\rho=0.2$, we
show the disorder-averaged probability distribution over the sample, for the
state just below zero in energy of each disorder realization. The state is
strongly localized on the $x=0$ and $x=39$ edges. Fig. 4 c) shows a
corresponding cut through b) for $y=20$, revealing exponential decay of the
probability density into the bulk for $\rho=0.2$, which is lost at $\rho=0.4$,
in correspondence with changes observed in Fig 4 a).

We also examine the robustness of edge state spin texture against disorder, as
shown in Fig. 4 d). We see the strong polarization of the spin along the edge
for $\alpha_{0}=0.2$, but also see a subdominant helicity in the spin
components perpendicular to the edge. Additionally, there is some evidence of
counter-polarized spin texture just away from the edge. These features are
present before disorder-average as shown in the Supp. Mat., Section VI,
“Boundary mode real-space spin textures”, which also, notably, exhibit a real-
space, extended skyrmionic texture in the bulk. This indicates that the
ingredients required to realize topological skyrmion phases, even in
relatively simple toy models, also yield extended real-space magnetic
skyrmions typically studied in the continuum or in far more complex lattice
models [49]. We defer a more complete study of these rich real-space magnetic
textures to a future work.

Discussion and conclusion— We introduce a class of minimal Bloch Hamiltonians
generically realizing type-II topological phase transitions of topological
skyrmion phases [26, 42, 43] and the quantum skyrmion Hall effect [27]. We
also present evidence of gapless edge states due to non-trivial skyrmion
number corresponding to bands which do not traverse the bulk gap, and instead
yield gaplessness due to their overlap in energy and localization on opposite
edges of the system. This gaplessness is topologically-robust against
disorder, and corresponds to distinctive edge state spin textures, including a
dominant polarization of spin along the edge, and a subdominant feature that
appears helical in nature. These arc-like gapless boundary modes also occur
generically in other models [27, 48], indicating results presented here are
broadly-relevant to topological skyrmion phases and the quantum skyrmion Hall
effect.

###### Acknowledgements.

Acknowledgements—This research was supported in part by the National Science
Foundation under Grants No.NSF PHY-1748958 and PHY-2309135, and undertaken in
part at Aspen Center for Physics, which is supported by National Science
Foundation grant PHY-2210452.

## References

  * Klitzing  _et al._ [1980] K. v. Klitzing, G. Dorda, and M. Pepper, New method for high-accuracy determination of the fine-structure constant based on quantized $\mathrm{H}$all resistance, Phys. Rev. Lett. 45, 494 (1980).
  * Tsui  _et al._ [1982] D. C. Tsui, H. L. Stormer, and A. C. Gossard, Two-dimensional magnetotransport in the extreme quantum limit, Phys. Rev. Lett. 48, 1559 (1982).
  * Laughlin [1983a] R. B. Laughlin, Anomalous quantum hall effect: An incompressible quantum fluid with fractionally charged excitations, Phys. Rev. Lett. 50, 1395 (1983a).
  * Laughlin [1983b] R. B. Laughlin, Anomalous quantum hall effect: An incompressible quantum fluid with fractionally charged excitations, Phys. Rev. Lett. 50, 1395 (1983b).
  * Kallin and Halperin [1984] C. Kallin and B. Halperin, Excitations from a filled landau level in the two-dimensional electron gas, Physical Review B 30, 5655 (1984).
  * Zhang  _et al._ [1989] S. C. Zhang, T. H. Hansson, and S. Kivelson, Effective-field-theory model for the fractional quantum hall effect, Phys. Rev. Lett. 62, 82 (1989).
  * Halperin  _et al._ [1993] B. I. Halperin, P. A. Lee, and N. Read, Theory of the half-filled landau level, Physical Review B 47, 7312 (1993).
  * Halperin [1982] B. I. Halperin, Quantized hall conductance, current-carrying edge states, and the existence of extended states in a two-dimensional disordered potential, Physical Review B 25, 2185 (1982).
  * Halperin [1984] B. I. Halperin, Statistics of quasiparticles and the hierarchy of fractional quantized hall states, Physical Review Letters 52, 1583 (1984).
  * Wen [1991] X.-G. Wen, Gapless boundary excitations in the quantum hall states and in the chiral spin states, Physical Review B 43, 11025 (1991).
  * Jain [1989] J. K. Jain, Composite-fermion approach for the fractional quantum hall effect, Phys. Rev. Lett. 63, 199 (1989).
  * Bernevig  _et al._ [2002] B. A. Bernevig, C.-H. Chern, J.-P. Hu, N. Toumbas, and S.-C. Zhang, Effective field theory description of the higher dimensional quantum hall liquid, Annals of Physics 300, 185 (2002).
  * Nayak  _et al._ [2008] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Non-abelian anyons and topological quantum computation, Rev. Mod. Phys. 80, 1083 (2008).
  * Sarma  _et al._ [2015] S. D. Sarma, M. Freedman, and C. Nayak, Majorana zero modes and topological quantum computation, npj Quantum Information 1, 15001 (2015).
  * Freedman  _et al._ [2002] M. H. Freedman, M. J. Larsen, and Z. Wang, The Two-Eigenvalue Problem and Density¶of Jones Representation of Braid Groups, Communications in Mathematical Physics 228, 177 (2002).
  * Kitaev [2001] A. Y. Kitaev, Unpaired majorana fermions in quantum wires, Physics-Uspekhi 44, 131 (2001).
  * Ryu  _et al._ [2010] S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. W. Ludwig, Topological insulators and superconductors: tenfold way and dimensional hierarchy, New Journal of Physics 12, 065010 (2010).
  * Schnyder  _et al._ [2008] A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Classification of topological insulators and superconductors in three spatial dimensions, Phys. Rev. B 78, 195125 (2008).
  * Qi  _et al._ [2008] X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Topological field theory of time-reversal invariant insulators, Phys. Rev. B 78, 195424 (2008).
  * Kitaev  _et al._ [2009] A. Kitaev, V. Lebedev, and M. Feigel’man, Periodic table for topological insulators and superconductors, in  _AIP Conference Proceedings_ (AIP, 2009).
  * Altland and Zirnbauer [1997] A. Altland and M. R. Zirnbauer, Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures, Phys. Rev. B 55, 1142 (1997).
  * Chiu  _et al._ [2016] C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Classification of topological quantum matter with symmetries, Rev. Mod. Phys. 88, 035005 (2016).
  * Hsieh  _et al._ [2008] D. Hsieh, D. Qian, L. Wray, Y. Xia, Y. S. Hor, R. J. Cava, and M. Z. Hasan, A topological dirac insulator in a quantum spin hall phase, Nature 452, 970 (2008).
  * Chang  _et al._ [2013] C.-Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang, M. Guo, K. Li, Y. Ou, P. Wei, L.-L. Wang, Z.-Q. Ji, Y. Feng, S. Ji, X. Chen, J. Jia, X. Dai, Z. Fang, S.-C. Zhang, K. He, Y. Wang, L. Lu, X.-C. Ma, and Q.-K. Xue, Experimental observation of the quantum anomalous hall effect in a magnetic topological insulator, Science 340, 167 (2013), https://www.science.org/doi/pdf/10.1126/science.1234414 .
  * Else  _et al._ [2021] D. V. Else, R. Thorngren, and T. Senthil, Non-fermi liquids as ersatz fermi liquids: General constraints on compressible metals, Phys. Rev. X 11, 021005 (2021).
  * Cook [2023a] A. M. Cook, Topological skyrmion phases of matter, Journal of Physics: Condensed Matter 35, 184001 (2023a).
  * Cook [2023b] A. M. Cook, Quantum skyrmion Hall effect (2023b), arXiv:2305.18626 [cond-mat.str-el] .
  * Bogdanov and Yablonskii [1989] A. N. Bogdanov and D. Yablonskii, Thermodynamically stable “vortices” in magnetically ordered crystals. the mixed state of magnets, Zh. Eksp. Teor. Fiz 95, 178 (1989).
  * Yu  _et al._ [2010] X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa, and Y. Tokura, Real-space observation of a two-dimensional skyrmion crystal, Nature 465, 901 (2010).
  * Mühlbauer  _et al._ [2009] S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, Skyrmion lattice in a chiral magnet, Science 323, 915 (2009), https://www.science.org/doi/pdf/10.1126/science.1166767 .
  * Schulz  _et al._ [2012] T. Schulz, R. Ritz, A. Bauer, M. Halder, M. Wagner, C. Franz, C. Pfleiderer, K. Everschor, M. Garst, and A. Rosch, Emergent electrodynamics of skyrmions in a chiral magnet, Nature Physics 8, 301 (2012).
  * Parkin  _et al._ [2008] S. S. P. Parkin, M. Hayashi, and L. Thomas, Magnetic domain-wall racetrack memory, Science 320, 190 (2008), https://www.science.org/doi/pdf/10.1126/science.1145799 .
  * Sampaio  _et al._ [2013] J. Sampaio, V. Cros, S. Rohart, A. Thiaville, and A. Fert, Nucleation, stability and current-induced motion of isolated magnetic skyrmions in nanostructures, Nature Nanotechnology 8, 839 (2013).
  * Donnelly  _et al._ [2020] C. Donnelly, S. Finizio, S. Gliga, M. Holler, A. Hrabec, M. Odstrčil, S. Mayr, V. Scagnoli, L. J. Heyderman, M. Guizar-Sicairos, and J. Raabe, Time-resolved imaging of three-dimensional nanoscale magnetization dynamics, Nature Nanotechnology 15, 356 (2020).
  * Zhang  _et al._ [2020] X. Zhang, Y. Zhou, K. M. Song, T.-E. Park, J. Xia, M. Ezawa, X. Liu, W. Zhao, G. Zhao, and S. Woo, Skyrmion-electronics: writing, deleting, reading and processing magnetic skyrmions toward spintronic applications, Journal of Physics: Condensed Matter 32, 143001 (2020).
  * Maccariello  _et al._ [2018] D. Maccariello, W. Legrand, N. Reyren, K. Garcia, K. Bouzehouane, S. Collin, V. Cros, and A. Fert, Electrical detection of single magnetic skyrmions in metallic multilayers at room temperature, Nature Nanotechnology 13, 233 (2018).
  * Jiang  _et al._ [2017] W. Jiang, X. Zhang, G. Yu, W. Zhang, X. Wang, M. Benjamin Jungfleisch, J. Pearson, X. Cheng, O. Heinonen, K. L. Wang, Y. Zhou, A. Hoffmann, and S. te Velthuis, Direct observation of the skyrmion Hall effect, Nature Physics 13, 162 (2017).
  * Zang  _et al._ [2011] J. Zang, M. Mostovoy, J. H. Han, and N. Nagaosa, Dynamics of skyrmion crystals in metallic thin films, Phys. Rev. Lett. 107, 136804 (2011).
  * Neubauer  _et al._ [2009] A. Neubauer, C. Pfleiderer, B. Binz, A. Rosch, R. Ritz, P. G. Niklowitz, and P. Böni, Topological hall effect in the $a$ phase of mnsi, Phys. Rev. Lett. 102, 186602 (2009).
  * Bruno  _et al._ [2004] P. Bruno, V. Dugaev, and M. Taillefumier, Topological hall effect and berry phase in magnetic nanostructures, Physical review letters 93, 096806 (2004).
  * Lee  _et al._ [2009] M. Lee, W. Kang, Y. Onose, Y. Tokura, and N. P. Ong, Unusual hall effect anomaly in mnsi under pressure, Phys. Rev. Lett. 102, 186601 (2009).
  * Liu  _et al._ [2023] S.-W. Liu, L.-k. Shi, and A. M. Cook, Defect bulk-boundary correspondence of topological skyrmion phases of matter, Phys. Rev. B 107, 235109 (2023).
  * Flores-Calderon and Cook [2023] R. Flores-Calderon and A. M. Cook, Time-reversal invariant topological skyrmion phases (2023), arXiv:2306.14204 [cond-mat.mes-hall] .
  * Ueno  _et al._ [2013] Y. Ueno, A. Yamakage, Y. Tanaka, and M. Sato, Symmetry-protected majorana fermions in topological crystalline superconductors: Theory and application to $\mathrm{Sr}_{2}\mathrm{RuO}_{4}$, Phys. Rev. Lett. 111, 087002 (2013).
  * Ng and Sigrist [2000] K. K. Ng and M. Sigrist, The role of spin-orbit coupling for the superconducting state in sr2ruo4, Europhysics Letters (EPL) 49, 473 (2000).
  * Qi  _et al._ [2006] X.-L. Qi, Y.-S. Wu, and S.-C. Zhang, Topological quantization of the spin hall effect in two-dimensional paramagnetic semiconductors, Phys. Rev. B 74, 085308 (2006).
  * Sticlet  _et al._ [2012] D. Sticlet, F. Piéchon, J.-N. Fuchs, P. Kalugin, and P. Simon, Geometrical engineering of a two-band chern insulator in two dimensions with arbitrary topological index, Phys. Rev. B 85, 165456 (2012).
  * Winter  _et al._ [2023] J. H. Winter, R. Ay, B. Braunecker, and A. M. Cook, Observable-enriched entanglement (2023), (in progress) .
  * Peters  _et al._ [2023] R. Peters, J. Neuhaus-Steinmetz, and T. Posske, Quantum skyrmion hall effect in $f$-electron systems, Phys. Rev. Res. 5, 033180 (2023).

Supplemental material for “Type-II topological phase transitions of
topological skyrmion phases”
Reyhan Ay1,2,3, Joe H. Winter1,2,4, and Ashley M. Cook1,2,∗

1Max Planck Institute for Chemical Physics of Solids, Nöthnitzer Strasse 40,
01187 Dresden, Germany
2Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Strasse
38, 01187 Dresden, Germany
3Izmir Institute of Technology, Gülbahçe Kampüsü, 35430 Urla Izmir, Türkiye
4SUPA, School of Physics and Astronomy, University of St. Andrews, North
Haugh, St. Andrews KY16 9SS, UK ∗Electronic address<EMAIL_ADDRESS>
##  S1 I. Tight-binding model for mirror subsectors of superconducting Sr2RuO4
Bogoliubov de Gennes Hamiltonian and associated spin operators

A two-dimensional model for Sr2RuO4 [44] in the $x-y$ plane consists of a
Bogoliubov de Gennes Hamiltonian with a spin half degree of freedom, three-
fold $t_{2g}$ orbital degree of freedom, and particle-hole degree of freedom,
corresponding to $12\times 12$ matrix representation. This Hamiltonian is
invariant under a mirror operation taking $z$ to $-z$, and may be block-
diagonalized by going to the basis in which the mirror operator matrix
representation is diagonal. We present the Bloch Hamiltonians of these mirror
subsectors here, each with $6\times 6$ matrix representation. We consider the
case where each mirror subsector Bloch Hamiltonian itself possesses particle-
hole symmetry [44, 26] corresponding to invariance under charge conjugation
$\mathcal{C}$ ($\mathcal{C}^{2}=-1$), as well as a generalized particle-hole
symmetry corresponding to invariance under operation
$\mathcal{C}^{\prime}=\mathcal{C}\mathcal{I}$, where $\mathcal{I}$ is spatial
inversion.

One mirror subsector Bloch Hamiltonian ($6\times 6$ matrix representation)
takes the following form:

$\displaystyle\mathcal{H}_{1}(\boldsymbol{k})=\begin{pmatrix}h_{1}(\boldsymbol{k})&\Delta_{1}(\boldsymbol{k})\\\
\Delta_{1}^{\dagger}(\boldsymbol{k})&-h_{1}^{*}(-\boldsymbol{k})\end{pmatrix}.$
(S1)

with basis vector
$\Psi^{\dagger}_{\boldsymbol{k}}=\left(c_{\boldsymbol{k},xy,\uparrow},c_{\boldsymbol{k},yz,\downarrow},c_{\boldsymbol{k},xz,\downarrow},c^{\dagger}_{\boldsymbol{-k},xy,\uparrow},c^{\dagger}_{\boldsymbol{-k},yz,\downarrow},c^{\dagger}_{\boldsymbol{-k},xz,\downarrow}\right)$
and

$\displaystyle
h_{1}(\boldsymbol{k})=\begin{pmatrix}-\mu_{B}H_{z}+\mu^{\prime}+4t_{3}\cos(k_{x})\cos(k_{y})+2t_{2}\left[\cos(k_{x})+\cos(k_{y})\right]&-i\lambda&-\lambda\\\
i\lambda&\mu_{B}H_{z}+\mu+2t_{1}\cos(k_{x})&i\lambda+4t_{4}\sin(k_{x})\sin(k_{y})\\\
-\lambda&-i\lambda+4t_{4}\sin(k_{x})\sin(k_{y})&\mu_{B}H_{z}+\mu+2t_{1}\cos(k_{y})\end{pmatrix}$
(S2)

and

$\displaystyle\Delta_{1}(\boldsymbol{k})=\begin{pmatrix}{\Delta_{0}\over
2}\left[i\sin(k_{x})-\sin(k_{y})\right]&0&0\\\ 0&{\Delta_{0}\over
2}\left[i\sin(k_{x})+\sin(k_{y})\right]&0\\\ 0&0&{\Delta_{0}\over
2}\left[i\sin(k_{x})+\sin(k_{y})\right]\end{pmatrix}.$ (S3)

Here, $\mu_{B}H_{z}$ is strength of applied Zeeman field in the
$\hat{z}$-direction, $\mu^{\prime}$ is an energy offset between $xy$ and
$xz$/$yz$ orbitals, $\mu$ is chemical potential, $t_{1}$, $t_{2}$, $t_{3}$,
$t_{4}$ are hopping integrals, and $\lambda$ is the atomic spin-orbit coupling
strength.

The basis of this mirror subsector may be used to identify a spin
representation provided in past work on topological skyrmion phases [26]. The
spin representation is
$\boldsymbol{\tilde{S}}=\mathrm{diag}\left(\boldsymbol{S},-\boldsymbol{S}^{*}\right)$,
where $\boldsymbol{S}$ is the vector of spin operators, $\mathcal{S}_{x}$,
$\mathcal{S}_{y}$, and $\mathcal{S}_{z}$, for the particle sector, taken to be

$\displaystyle S_{x}={1\over\sqrt{2}}\begin{pmatrix}0&1&1\\\ 1&0&1\\\
1&1&0\end{pmatrix},\hskip
14.22636ptS_{y}={1\over\sqrt{2}}\begin{pmatrix}0&-i&-i\\\ i&0&-i\\\
i&i&0\end{pmatrix},\hskip
14.22636ptS_{z}={1\over\sqrt{2}}\begin{pmatrix}2&0&0\\\ 0&-1&0\\\
0&0&-1\end{pmatrix}.$ (S4)

For completeness, we also list the second mirror subsector Bloch Hamiltonian
although we do not use it as the basis for toy models in this work. This Bloch
Hamiltonian takes the form

$\displaystyle\mathcal{H}_{2}(\boldsymbol{k})=\begin{pmatrix}h_{2}(\boldsymbol{k})&\Delta_{2}(\boldsymbol{k})\\\
\Delta_{2}^{\dagger}(\boldsymbol{k})&-h_{2}^{*}(-\boldsymbol{k})\end{pmatrix}.$
(S5)

with basis vector
$\Phi^{\dagger}_{\boldsymbol{k}}=\left(c_{\boldsymbol{k},xy,\downarrow},c_{\boldsymbol{k},yz,\uparrow},c_{\boldsymbol{k},xz,\uparrow},c^{\dagger}_{\boldsymbol{-k},xy,\downarrow},c^{\dagger}_{\boldsymbol{-k},yz,\uparrow},c^{\dagger}_{\boldsymbol{-k},xz,\uparrow}\right)$
and

$\displaystyle
h_{2}(\boldsymbol{k})=\begin{pmatrix}\mu_{B}H_{z}+\mu^{\prime}+4t_{3}\cos(k_{x})\cos(k_{y})+2t_{2}\left[\cos(k_{x})+\cos(k_{y})\right]&-i\lambda&\lambda\\\
i\lambda&-\mu_{B}H_{z}+\mu+2t_{1}\cos(k_{x})&-i\lambda+4t_{4}\sin(k_{x})\sin(k_{y})\\\
\lambda&i\lambda+4t_{4}\sin(k_{x})\sin(k_{y})&-\mu_{B}H_{z}+\mu+2t_{1}\cos(k_{y})\end{pmatrix}$
(S6)

and

$\displaystyle\Delta_{2}(\boldsymbol{k})=\begin{pmatrix}{\Delta_{0}\over
2}\left[i\sin(k_{x})+\sin(k_{y})\right]&0&0\\\ 0&{\Delta_{0}\over
2}\left[i\sin(k_{x})-\sin(k_{y})\right]&0\\\ 0&0&{\Delta_{0}\over
2}\left[i\sin(k_{x})-\sin(k_{y})\right]\end{pmatrix}.$ (S7)

Note that $h_{1}(-\boldsymbol{k})=h_{1}(\boldsymbol{k})$ and
$h_{2}(-\boldsymbol{k})=h_{2}(\boldsymbol{k})$. The spin representation for
this Hamiltonian is related by symmetry to the representation for
$\mathcal{H}_{1}(\boldsymbol{k})$ and also provided in earlier work [26].

##  S2 II. Additional characterization of bulk topology for a second choice of
$\boldsymbol{d}(\boldsymbol{k})$ vector

Here, we present a second $\boldsymbol{d}(\boldsymbol{k})$-vector for use in
the six-band Bloch Hamiltonian toy model for topological skyrmion phases of
matter, and generate phase diagrams using this second
$\boldsymbol{d}(\boldsymbol{k})$-vector similar to those shown in Fig. 1 in
the main text. We consider the $\boldsymbol{d}(\boldsymbol{k})$-vector of a
well-known two-band Bloch Hamiltonian for a Chern insulator, which realizes
non-trivial Chern number of $\pm 2$ for the two-band model [47],

$\boldsymbol{d}_{s}(\boldsymbol{k})=\langle
d_{x}(\boldsymbol{k}),d_{y}(\boldsymbol{k}),d_{z}(\boldsymbol{k})\rangle,$
(S8)

and

$\displaystyle d_{x}(\boldsymbol{k})$ $\displaystyle=\alpha\cos(k_{x}),$ (S9)
$\displaystyle d_{y}(\boldsymbol{k})$ $\displaystyle=\alpha\cos(k_{y}),$
$\displaystyle d_{z}(\boldsymbol{k})$ $\displaystyle=\beta\cos(k_{x}+k_{y})$

Figure S5: Phase diagrams characterizing bulk topology of Hamiltonian Eq. 1
with $\boldsymbol{d}(\boldsymbol{k})$-vector Eq. S8 for half-filling,
including a) total Chern number $\mathcal{C}$, b) skyrmion number
$\mathcal{Q}$, c) minimum direct bulk energy gap $\Delta E_{\mathrm{min}}$,
and d) minimum ground state spin magnitude $|\boldsymbol{S}|_{\mathrm{min}}$,
with each of these four quantities plotted vs. atomic spin-orbit coupling
(SOC) strength $\lambda$ and hopping parameter $\alpha$. Other parameters are
fixed at $\Delta=0.6$, and $\beta=1$

##  S3 III. Additional slab spectra through type-II topological phase
transitions

Here, we present additional results on evolution of the spectrum and bulk-
boundary correspondence through type-II topological phase transitions for open
boundary conditions (OBC) in the $\hat{x}$\- ($\hat{y}$)-direction and
periodic boundary conditions (PBC) in the $\hat{y}$\- ($\hat{x}$)-direction.
Transition A for OBC in the $\hat{x}$-direction is shown in the main text,
Fig. 2, so here we show additional results in this case for Transitions B to D
in Figs. S6, S7, S8, respectively. We then show evolution of slab spectra for
OBC in the $\hat{y}$-direction for all Transitions A, B, C, and D in Figs. S9,
S10, S11, and S12, respectively.

###  S3.1 Transitions B to D, OBC in x

Here, we show evolution of slab spectra for OBC in the $\hat{x}$-direction for
Transitions B, C, and D in S6, S7, and S8, respectively.

Figure S6: Slab spectra for open boundary conditions in the
$\hat{x}$-direction, periodic boundary conditions in the $\hat{y}$-direction
across Transition B for parameters a) $\mathcal{M}=0.3625$ and
$\lambda=-0.75$, b) $\mathcal{M}=0.125$ and $\lambda=-0.7$, c)
$\mathcal{M}=-0.08875$ and $\lambda=-0.655$, d) $\mathcal{M}=-0.35$ and
$\lambda=-0.60$. Figure S7: Slab spectra for open boundary conditions in the
$\hat{x}$-direction, periodic boundary conditions in the $\hat{y}$-direction
across Transition C for parameters a) $\mathcal{M}=-5.10$ and
$\lambda=0.9875$, b) $\mathcal{M}=-4.7$ and $\lambda=0.8875$, c)
$\mathcal{M}=-4.35$ and $\lambda=0.8$, d) $\mathcal{M}=-4.0$ and
$\lambda=0.7125$. Figure S8: Slab spectra for open boundary conditions in the
$\hat{x}$-direction, periodic boundary conditions in the $\hat{y}$-direction
across Transition D for parameters a) $\mathcal{M}=-0.9$ and $\lambda=0.35$,
b) $\mathcal{M}=-1.06$ and $\lambda=0.35$, c) $\mathcal{M}=-1.10$ and
$\lambda=0.35$, d) $\mathcal{M}=-1.90$ and $\lambda=0.35$.

###  S3.2 Transitions A to D, OBC in y

Here, we show evolution of slab spectra for OBC in the $\hat{y}$-direction for
all Transitions A, B, C, and D in Figs. S9, S10, S11, and S12, respectively.
As shown in Figs. S9, S10, and S11, Transitions A-C slab spectra for these
boundary conditions are qualitatively unchanged across the type-II transition
for parameter sets over which topological-robust gaplessness is gained/lost
for OBCs instead in the $\hat{x}$-direction, while the slab spectrum for
Transition D exhibits additional band-crossings at zero-energy across the
type-II transition also for these boundary conditions as shown in Fig. S12.

Figure S9: Slab spectra for open boundary conditions in the
$\hat{y}$-direction, periodic boundary conditions in the $\hat{x}$-direction
across Transition A for parameters a) $\mathcal{M}=-2.7$ and $\lambda=0.8$, b)
$\mathcal{M}=-2.6$ and $\lambda=0.8$, c) $\mathcal{M}=-2.558$ and
$\lambda=0.8$, d) $\mathcal{M}=-1.9$ and $\lambda=0.8$. Figure S10: Slab
spectra for open boundary conditions in the $\hat{y}$-direction, periodic
boundary conditions in the $\hat{x}$-direction across Transition B for
parameters a) $\mathcal{M}=0.3625$ and $\lambda=-0.75$, b) $\mathcal{M}=0.125$
and $\lambda=-0.7$, c) $\mathcal{M}=-0.08875$ and $\lambda=-0.655$, d)
$\mathcal{M}=-0.35$ and $\lambda=-0.60$. Figure S11: Slab spectra for open
boundary conditions in the $\hat{y}$-direction, periodic boundary conditions
in the $\hat{x}$-direction across Transition C for parameters a)
$\mathcal{M}=-5.10$ and $\lambda=0.9875$, b) $\mathcal{M}=-4.7$ and
$\lambda=0.8875$, c) $\mathcal{M}=-4.35$ and $\lambda=0.8$, d)
$\mathcal{M}=-4.0$ and $\lambda=0.7125$. Figure S12: Slab spectra for open
boundary conditions in the $\hat{y}$-direction, periodic boundary conditions
in the $\hat{x}$-direction across Transition D for parameters a)
$\mathcal{M}=-0.9$ and $\lambda=0.35$, b) $\mathcal{M}=-1.06$ and
$\lambda=0.35$, c) $\mathcal{M}=-1.10$ and $\lambda=0.35$, d)
$\mathcal{M}=-1.90$ and $\lambda=0.35$.

###  S3.3 Additional results for Transition D with open boundary conditions in
$\hat{y}$-direction

We show additional results highlighting particular features of the gapless
boundary states associated with the type-II topological phase transition in
Fig. S13. In Fig. S13 a), states are present within the bulk gap in the
vicinity of $k_{x}=0$, but do not overlap to yield gaplessness. With even more
negative $\mathcal{M}$, gaplessness occurs through crossing of these in-gap
bands at a single $\boldsymbol{k}$-point as shown in Fig. S13 b). At even more
negative $\mathcal{M}$, gaplessness persists via crossing of the in-gap bands
around $k_{x}=0$ at two separate $\boldsymbol{k}$-points, while gaplessness
near $k_{x}=-\pi/a$ is lost.

Figure S13: Slab spectra for open boundary conditions in the
$\hat{y}$-direction, periodic boundary conditions in the $\hat{x}$-direction
highlighting edge state band connectivity across Transition D for parameters
a) $\mathcal{M}=-1.30$ and $\lambda=0.35$, b) $\mathcal{M}=-1.60$ and
$\lambda=0.35$, c) $\mathcal{M}=-1.90$ and $\lambda=0.35$

##  S4 IV. Effect of atomic spin-orbit coupling strength on transition A

Here, we show the bulk skyrmion number phase diagram Fig. 1 b) in Fig. S14 a),
highlighting four points in phase space with corresponding slab spectra for
OBC in the $\hat{x}$-direction shown in Fig. S14 b), c), d) and e),
respectively. Each slab spectrum depicts edge states which just touch rather
than overlap, with gaplessness lost through fine-tuned increase in mass
parameter $\mathcal{M}$. In this region, the separation in phase space between
when the skyrmion number $\mathcal{Q}$ changes in the bulk and loss of
gaplessness at the edge increases with increasing atomic spin-orbit coupling
strength $\lambda$.

Figure S14: Onset of gaplessness with varying atomic spin-orbit coupling
strength. a) depicts the skyrmion number $\mathcal{Q}$ vs. atomic spin orbit
coupling strength $\lambda$ and mass term $\mathcal{M}$ also shown in Fig. 1
b), with phase space coordinates for four slab spectra computed with open
boundary conditions in the $\hat{x}$-direction and periodic boundary
conditions in the $\hat{y}$-direction shown in b) for $\mathcal{M}=-2.136$,
$\lambda=0.68$, c) $\mathcal{M}=-2.069$, $\lambda=0.66$, d)
$\mathcal{M}=-2.001$, $\lambda=0.64$, e) $\mathcal{M}=-1.934$, $\lambda=0.62$.

##  S5 V. Additional edge state spin textures through type-II topological
phase transitions

In this section, we show evolution of the edge state spin textures through
four type-II topological phase transitions A to D shown in Figs. S15, S16,
S17, S18, respectively.

###  S5.1 Transition A

Figure S15: Edge state spin textures, with associated slab spectra, for open
boundary conditions in the $\hat{x}$-direction, periodic boundary conditions
in the $\hat{y}$-direction, for Transition A. Parameter sets are a)
$\mathcal{M}=-2.7$ and $\lambda=0.8$, b) $\mathcal{M}=-2.6$ and $\lambda=0.8$,
c) $\mathcal{M}=-2.558$ and $\lambda=0.8$, d) $\mathcal{M}=-2.4$ and
$\lambda=0.8$.

###  S5.2 Transition B

Figure S16: Edge state spin textures, with associated slab spectra, for open
boundary conditions in the $\hat{x}$-direction, periodic boundary conditions
in the $\hat{y}$-direction, for Transition B. Parameter sets are a)
$\mathcal{M}=0.3625$ and $\lambda=-0.75$, b) $\mathcal{M}=0.125$ and
$\lambda=-0.7$, c) $\mathcal{M}=-0.08875$ and $\lambda=-0.655$, d)
$\mathcal{M}=-0.35$ and $\lambda=-0.60$.

###  S5.3 Transition C

Figure S17: Edge state spin textures, with associated slab spectra, for open
boundary conditions in the $\hat{x}$-direction, periodic boundary conditions
in the $\hat{y}$-direction, for Transition C. Parameter sets are a)
$\mathcal{M}=-5.10$ and $\lambda=0.9875$, b) $\mathcal{M}=-4.7$ and
$\lambda=0.8875$, c) $\mathcal{M}=-4.35$ and $\lambda=0.8$, d)
$\mathcal{M}=-4.0$ and $\lambda=0.7125$.

###  S5.4 Transition D

Figure S18: Edge state spin textures, with associated slab spectra, for open
boundary conditions in the $\hat{x}$-direction, periodic boundary conditions
in the $\hat{y}$-direction, for Transition D. Parameter sets are a)
$\mathcal{M}=-0.9$ and $\lambda=0.35$, b) $\mathcal{M}=-1.06$ and
$\lambda=0.35$, c) $\mathcal{M}=-1.10$ and $\lambda=0.35$, d)
$\mathcal{M}=-1.90$ and $\lambda=0.35$.

##  S6 VI. Boundary mode real-space spin textures

Here, we show real-space spin textures for individual eigenstates of
Hamiltonian Eq. 1 without disorder ($\rho=0$) for OBCs in each of the
$\hat{x}$\- and $\hat{y}$-directions in correspondence with Fig. 4 in the main
text.

Figure S19: Real-space spin textures of eigenstates number a) i = 4800 at
eigenvalue -0.00121, b) i = 4799 at eigenvalue -0.01217, c) i = 4798 at
eigenvalue -0.02288, and d) i = 4797 at eigenvalue -0.02486.