# The topological magnon bands in the Flux state in Sashtry-Sutherland lattice

Dhiman Bhowmick  Pinaki Sengupta  School of Physical and Mathematical
Sciences, Nanyang Technological University, 21 Nanyang Link, Singapore 637371,
Singapore

###### Abstract

We investigate low energy magnon excitations above the non-collinear flux
state and non-coplanar canted flux state in a Heisenberg anti-ferromagnet with
Dzyaloshinskii-Moriya interaction (DMI) on a Sashtry-Sutherland lattice. While
previous studies have shown the presence of topological magnetic excitation in
the dimer and ferromagnetic phases on the Shastry-Sutherland lattice, our
results establish the non-trivial topology of magnons in the anti-
ferromagnetic flux and canted flux states. Our results uncover the existence
of a multitude of topological phase transitions in the magnon sector –
evidenced by the changing Chern numbers of the single magnon bands – as the
Hamiltonian parameters are varied, even when the ground state remains
unchanged. The thermal Hall conductivity is calculated and its derivative is
shown to exhibit a logarithmic divergence at the phase transitions,
independent of the type of band touching involved. This may provide a useful
means to identify the energy at which the transition occurs. Finally, we
propose the way to realize the studied model in a practical material.

PACS numbers

###### pacs:

Valid PACS appear here

##  I Introduction

The study of topological phases of matter has gained widespread interest
during the past decade. While topological phases are realized in both
fermionicKane and Mele (2005); Zhang  _et al._ (2005) and bosonic
systemsOwerre (2016a); Karzig  _et al._ (2015); Kim  _et al._ (2016a); Zhang
_et al._ (2018a), much of the advancement (theoretical investigations and
experimental realizations of topological phases) has been confined to
fermionic systems; study of bosonic topological phases have attracted
widespread interest relatively recently. This is partly due to the fact that
topological character of ground state phases of (non-interacting) fermions is
readily identified from the properties of the energy bands and there exist
well developed experimental probes to detect them. On the other hand, the
ground state of bosonic systems is often a condensate and topological
character is manifested in low lying excitationsOwerre (2016a); Kim  _et al._
(2016a); Zhang  _et al._ (2018a).

Quantum magnets have served as a versatile test bed for realizing novel
bosonic phases, including bosonic topological phases. The topological
character of the magnetic phase is manifested through the behavior of magnons.
Magnons are charge neutral quasi-particle excitations in insulating magnetic
systems. Analogous to electrons in standard topological insulators, magnons in
magnetic insulators exhibit thermal Hall effectShuichi and Akihiro (2017),
spin-Nernst effectCheng  _et al._ (2016a); Zyuzin and Kovalev (2016a); Shiomi
_et al._ (2017); Kim  _et al._ (2016b), and magnon-driven spin Seebeck
effectJin  _et al._ (2015). The interest in these systems is driven by both
fundamental reasons and potential for technological applications. The recent
use of Skyrmions in spintronics for efficient magnetic storage and read/write
devices with minimal Joule heating effectFert  _et al._ (2013); Jonietz  _et
al._ (2010) underscores the potential practical applications of topologically
no-trivial magnetic states. The wide range of quantum magnets with varying
interactions and lattice structures as well as the ability to control the
number of quantized excitations with an external magnetic field make them
ideal for exploring novel magnetic phases. Geometrically frustrated quantum
magnets are particulalry promising in realizing and controlling topologically
non-trivial spin texturesLee  _et al._ (2008); Kim  _et al._ (2018); Shahzad
and Sengupta (2017); Tomiyasu  _et al._ (2012); Sagayama  _et al._ (2013);
Disseler (2014). The interplay between competing interactions, geometric
frustration and external magnetic field result in a wide variety of magnetic
phases that are not commonly observed in their non-frustrated counterparts. In
most cases, topological excitations in quantum magnets are driven
Dzyaloshinskii-Moriya interaction (DMI)Owerre (2016b); Zyuzin and Kovalev
(2016b); Cheng  _et al._ (2016b); Owerre (2017a); Lee  _et al._ (2018),
although topological magnon bands can exist without DMI as well Owerre
(2017b); Kim  _et al._ (2019) due to non-coplaner chiral spin-texture. Most
strikingly, the change in spin-texture by changing the parameters in the non-
coplanar spin systems, gives rise to variety of topological phases in the same
systemKim  _et al._ (2019); Laurell and Fiete (2018).

The Shastry-Sutherland model is a paradigmatic model for the study of
frustrated magnetism. Since the degree of frustration can be tuned by varying
the ratio of the diagonal and axial bonds, the model exhibits a wide range of
novel magnetic phasesMiyahara and Ueda (1999); Koga and Kawakami (2000);
Läuchli  _et al._ (2002). The existence of a number of materials with
underlying SS geometry of the magnetic ions offers the prospect of observing
theoretically predicted phases and phenomena in real materialsKageyama  _et
al._ (1999); Siemensmeyer  _et al._ (2008); Wierschem and Sengupta (2013); Ye
_et al._ (2017). Since DMI is ubiquitous in all of these materials, it is
natural to supplement the canonical SS model with DMI. This results in an even
richer variety of magnetic orderings including colinear, coplanar and non-
coplanar spin configurations, several of which host topological magnonsShahzad
and Sengupta (2017, 2017). Previous studies of topological magnons for the SS
lattice were restricted to the dimer Romhányi  _et al._ (2015a) and the ferro-
magnetic phases Malki and Uhrig (2019). In recent past there has been a
growing interest in studying topological magnons in non-collinear spin
configurations in frustrated lattices. Here we present the results of our
investigation of topological magnons in the recently proposed flux state, that
is stabilised by DMI perpendicular to the lattice plane in the SS
latticeShahzad and Sengupta (2017).In presence of in plane DMI, this evolves
to the canted flux state – the resulting magnon bands in an external
longitudinal magnetic field, carry non-zero Chern numbers that determine the
topological character of the magnon bands. Varying the different components of
the DMI result in a sequence of topological phase transitions where the Chern
number for the magnon bands change over a wide range of possible values.

The topological phase transitions in the system can be detected from the first
derivative of the thermal Hall conductivity which exhibits a logarithmic
divergence at the transitionKim  _et al._ (2019); Zhang  _et al._ (2010). The
peak height of the logarithmic divergence increases with temperature following
an algebraic relation. The interpolation of the first derivative of the
thermal Hall conductance as a function of temperature yields information on
the nature of band touching (gap closing) at the phase transition.

##  II Model Hamiltonian and Method

Figure 1: (color online) The Heisenberg and Dzyaloshinskii–Moriya(DM)
interactions on Shastry-Sutherland lattice.

The model Hamiltonian is given by,

$\displaystyle\pazocal{H}=$ $\displaystyle J\sum_{\left\langle\left\langle
i,j\right\rangle\right\rangle}\mathbb{S}_{i}\cdot\mathbb{S}_{j}+J^{\prime}\sum_{\left\langle
i,j\right\rangle}\mathbb{S}_{i}\cdot\mathbb{S}_{j}$
$\displaystyle+\mathbb{D}\cdot\sum_{\left\langle\left\langle
i,j\right\rangle\right\rangle}\left(\mathbb{S}_{i}\times\mathbb{S}_{j}\right)+\mathbb{D^{\prime}}\cdot\sum_{\left\langle
i,j\right\rangle}\left(\mathbb{S}_{i}\times\mathbb{S}_{j}\right)$
$\displaystyle-B\sum_{i}\mathbb{S}_{i}^{z},$ (1)

where $J$ and $J^{\prime}$ ($>0$) are anti-ferromagnetic Heisenberg spin-
exchange on the axial and diagonal bonds of the SS lattice respectively
(henceforth referred to as SS bonds). $\mathbb{D}$ and $\mathbb{D^{\prime}}$
are DM-vectors on the SS bonds as shown in Fig. 1. The nature of DM
interactions are chosen according to the symmetry constraints of the SS
lattice (Romhányi  _et al._ , 2015a; McClarty  _et al._ , 2017; Malki and
Schmidt, 2017). The origin of the DM interactions are further discussed in the
section Sec.IV.

We start with a classical ground state and investigate quantized low energy
excitations, magnons, focusing on identifying any topological character. The
classical ground state of the Hamiltonian Eq.1 is derived by replacing the
local spin moments by classical vectors of unit magnitude. The state of each
spin is specified by the polar co-ordinates,

$\mathbb{S}=S(\sin(\theta)\cos(\phi),\sin(\theta)\sin(\phi),\cos(\theta)).$
(2)

The ground state spin configuration is obtained by minimizing the energy of
the Hamiltonian w.r.t the angles $\theta$ and $\phi$. For the ground state
phases of interest, viz., the flux and the canted flux states, the magnetic
unit cell is of the same size as the unit cell of the SS lattice and consists
of four sites as shown in Fig.1.

The classical phases are further discussed in Sec.III. Since we are interested
in primarily the flux and the (in plane DMI induced) canted flux states, we
start by identifying the parameter ranges where these are realized. The ground
state phases of the Hamiltonian Eq.(1) for classical spins has been
investigated in Ref.[Shahzad and Sengupta, 2017] – with no in-plane component
of DMI, the flux state (Fig.2(a)) is stabilized above a critical value of the
normal component of the DMI along the axial bonds, $D_{\perp}$(Fig.1) (the DMI
on the diagonal bonds are constrained by the symmetry of the lattice to lie on
the plane of the lattice). The continuous U(1) symmetry of the Hamiltonian is
spontaneously broken in the flux state and Fig.2(a) shows one of the
degenerate ground states. Interestingly, the flux state state is also realized
in the square lattice, but for a much stronger DMI. The geometric frustration
of the SS lattice facilitates the appearance of the flux state for a more
moderate (and realistic) strength of DM interaction. However, the symmetry of
the SS lattice allows for in-plane DMI components of DMI, denoted in this work
by $D$ on the diagonal bonds and $D_{||,s}$, $D_{||,ns}$ on the axial bonds
(see Fig.1). Any non-zero in-plane DMI tilts the spins out of plane keeping
the in plane spin component of nearest-neighbour sites perpendicular to each
other. The spins on the two distinct diagonal bonds cant in opposite direction
– for one of the diagonals, the spins cant out of the plane, whereas for the
other diagonal, they cant into the plane of the lattice. The in-plane
components are aligned along the diagonal bonds, as depicted in Fig.3(a). This
ground state spin configuration is referred to as the canted-flux state. In
the presence of the in-plane components of DMI, the U(1) symmetry of the
Hamiltonian is explicitly broken and there is no spontaneous breaking of U(1)
symmetry in the in-plane DMI driven canted flux state.

To study the excitations above magnetic ground state, we have used the
linearized Holstein-Primakoff transformationOwerre (2016a, 2017a). The
Holstein-Primakoff transformation is a versatile and extensively used approach
to study low energy magnon excitations above magnetically ordered ground state
phases in quantum magnets. Previous implementations of this method have
largely been restricted to collinear magnetic orderings Owerre (2016a); Malki
and Uhrig (2019); Kim  _et al._ (2016a); Owerre (2017a); Cheng  _et al._
(2016b); Zhang  _et al._ (2018b); Lee  _et al._ (2018) although there have
been recent attempts to extend it to 120∘ non-collinear magnetic order in
triangular lattices Kim  _et al._ (2018). In this work, we have extended the
Holstein-Primakoff approach to study magnon excitations above complex magnetic
orders with longer periodicity. Here we present a brief discussion the method.
First, the local co-ordinate axis at each site of the lattice is rotated such
that the $S_{z}$ axis is aligned along the local spin direction. For low
temperature excitations the linearized Holstein-Primakoff transformation is
given by,
$\hat{S}^{\prime+}_{i,a}=\sqrt{2S}\hat{a}_{i},\;\hat{S}^{\prime-}_{i,a}=\sqrt{2S}\hat{a}^{\dagger}_{i},\;\hat{S}^{z}_{i,a}=S-\hat{a}^{\dagger}_{i}\hat{a}_{i},$
where, we consider $\hbar=1$ and $\hat{a}_{i}^{\dagger}(\hat{a}_{i})$
represent creation (annihilation) operators for quantized excitations above
the magnetic ground state at site $i$. These obey bosonic commutation
relations, $[\hat{a}_{i},\hat{a}_{j}^{\dagger}]I=\delta_{i,j}$ and
$[\hat{a}_{i},\hat{a}_{j}]=0=[\hat{a}_{i}^{\dagger},\hat{a}_{j}^{\dagger}]$.
Since the unit cell consists of 4 sites, there are four species of bosons
corresponding to each inequivalent lattice site. In the next section, the
detailed classical ground state and corresponding magnon bands and their
topological properties are discussed.

##  III  Results and Discussion

###  III.1 Flux State

Figure 2: (color online) (a) The flux state in Shastry-Sutherland lattice. (b)
The magnon band structure for $J=1.0,J^{\prime}=1.1,D_{\perp}=0.8,B_{z}=0.0$.
The inset of the figure shows the Brillouin zone and high symmetry lines. (c)
The magnon band structure for $J=1.0,J^{\prime}=1.1,D_{\perp}=0.8,B_{z}=0.5$.
Inset of the figure shows the magnified magnon band structure to show the
lifting of four-fold degeneracy at M-point. (d) The non-abelian Berry
curvature of lowest and second-lowest bands for
$J=1.0,J^{\prime}=1.1,D_{\perp}=0.8,B_{z}=0.5$. For convenience the negative
Berry-curvature which is concentrated at M-point is divided by 100, to
increase the visibility of the Berry-curvature distribution throughout the
Brillouin zone.

The flux state is stabilized as the ground state when the component of the DMI
normal to the plane of the lattice on the axial bonds, $D_{\perp}$, exceeds a
critical magnitude (e.g. when $J\approx J^{\prime}$, flux state is stable for
$D_{\perp}\gtrapprox 0.6J$) and the in-plane components vanish for all DMI
(the DMI on the diagonal bonds are constrained by symmetry requirements to be
strictly in-plane). The flux state is comprised of the nearest-neighbour spins
aligned perpendicular to each other and parallel to the plane of the lattice
(Fig.2(a)), which is energetically favored by the perpendicular DM-component
$D_{\perp}$. The state is characterized by the spontaneous breaking of
continuous U(1) spin-rotation symmetry about the z-axis. In the bosonic
(magnon) language, the ground state (flux state) is the vacuum and the
Holstein-Primakoff bosons represent quantized low energy excitations above
this ground state. The magnon bands at zero magnetic field are shown in the
Fig.2(b). At the $\Gamma$ point the lowest band becomes gapless revealing the
presence of Goldstone-mode associated with U(1) symmetry breaking. The bands
along line-$\overline{\text{MX}}$ are twofold degenerate – these can be
understood in terms of Kramer’s degeneracyLee  _et al._ (2018); Young and Kane
(2015). The operator
$\hat{m}_{2}=\left\\{\tilde{M}_{2}\tau|\bm{\delta}_{2}\right\\}$ commutes with
the Hamiltonian where $\tilde{M}_{2}$ is the reflection operator along the
$P_{2}$ axis and $\bm{\delta}_{2}$ is the translation by have lattice
parameter along the same axis, as shown in Fig.2(a); $\tau$ is the time-
reversal operator. But, $\hat{m}_{2}$ is not the symmetry operator for the
classical ground state shown in the Fig.2(a). Instead, symmetry operator for
the ground state is given by
$\hat{m}^{\prime}_{2}=\left\\{\tilde{M}_{2}e^{i\pi\hat{S}^{y}}\tau|\bm{\delta}_{2}\right\\}$,
which contains an additional rotation of spin by $\pi$ about the $y-$axis. On
the line-$\overline{\text{MX}}$ in the Brillouin zone,
$(\hat{m}^{\prime}_{2})^{2}=e^{ik_{x}}=-1$. Hence $\hat{m}^{\prime}_{2}$ is
anti-unitary operator with a squared value of $-1$, which in turn, results in
the Kramer’s degeneracy. $\hat{m}^{\prime}_{2}$ maps one Kramer’s degenerate
wavefunction along the $\overline{\text{MX}}$-line to the other Kramer’s
degenerate wavefunction along the $\overline{\text{M'X}}$-line. Further, the
symmetry operation
$\hat{c}_{2,1}=\left\\{\tilde{C}_{2,1}e^{i\pi\hat{S}^{y}}|\bm{\delta}_{1}\right\\}$,
maps the Kramer’s degenerate state from $\overline{\text{M'X}}$ to
$\overline{\text{MX}}$, where $\tilde{C}_{2,1}$ is the two-fold rotation about
the $P_{1}$ axis and $\bm{\delta}_{1}$ as shown in the Fig.2(a). Thus, the
band degeneracy along $\overline{\text{MX}}$-line is protected by symmetries
$\hat{m}_{2}$ and $\hat{c}_{2,1}$. Additionally, there is four-fold degeneracy
at the M-point due to the presence of the symmetry
$\hat{m}_{y}=\left\\{\tilde{M}_{y}\tau|\mathbb{0}\right\\}$, where
$\tilde{M}_{y}$ is the reflection operator about the $y-$axis as shown in
Fig.2(a). The symmetry $\hat{m}_{y}$ maps one pair of Kramer’s degenerate
state to the other pair of Kramer’s degenerate state at M-point. The Berry-
curvature is not well defined for bands, because the bands are degenerate.

The band structure in the presence of an external longitudinal magnetic field
is shown in Fig.2(c). Application of magnetic field produces a finite out of
plane spin component, which breaks the $\hat{m}_{y}$-symmetry. Thus, the four-
fold degeneracy reduces to two fold degeneracy at the M-point, as shown in the
inset of Fig.2(c). But the band sticking along $\overline{\text{MX}}$-line and
Goldstone modes at $\Gamma$-point are preserved. The degeneracy of the bands
prevents the calculation of the standard Berry curvature and Chern numbers of
the bands, although there is no symmetry constraints to make Berry curvature
zero. One can, however, calculate the non-abelian Berry-curvature and non-
abelian Chern number of the degenerate bands following the procedure outlined
in Ref.Matsuoka  _et al._ , 2005(see Appendix.C for details). The results for
the non-abelian Berry curvature of the lowest and second-lowest band is shown
in Fig.2(d). The negative part of Berry-curvature is highly concentrated at
the M-point and the positive berry curvature is distributed in the remaining
Brillouin zone. The inversion symmetry $\hat{i}$ of the system with an
inversion center at point-c in Fig.2(a) is also reflected in the Berry-
curvature. Our results show that the negative and positive contribution of the
non-abelian Berry curvature cancels, yielding a vanishing non-abelian Chern
number for the lower (or the upper) pair of bands.

###  III.2 Canted flux state

####  III.2.1 Topological magnon bands and topological phase diagram

Figure 3: (color online) (a) The spin configuration in 3Q-ordered phase.
Circle with dot and the circle with cross represent spin component out of
plane upward and downward directions respectively. (b) The magnon band
structure for
$J=1.0,J^{\prime}=1.1,D_{\perp}=0.8,D=0.2,D_{||,s}=0.05,D_{||,ns}=0.1,B_{z}=0.0$.
The inset of the figure shows the Brillouin zone and high symmetry lines and
high symmetry points. (c) The magnon band structure for
$J=1.0,J^{\prime}=1.1,D_{\perp}=0.8,D=0.2,D_{||,s}=0.05,D_{||,ns}=0.1,B_{z}=0.3$.
The Chern numbers of the bands are also shown in the Figure. (d) The Berry
curvature third Band for
$J=1.0,J^{\prime}=1.1,D_{\perp}=0.8,D=0.2,D_{||,s}=0.05,D_{||,ns}=0.1,B_{z}=0.3$.

When in-plane components of the DM-interactions (both along the axial as well
as the diagonal bonds), a canted-flux state, as shown in the Fig.3(a), is
realized as the ground state. The in-plane components of spins are directed
along the diagonal bonds; additionally, the spins acquire an out-of-plane
component. There is two-fold degeneracy in the ground state configuration –
one of them is shown in Fig.3(a) with the spins pointing inwards along the
diagonal bonds; the other degenerate state is obtained by flipping the spins
so that they point outward along the diagonal bonds and the out-of-plane spin
components are also flipped. The spins are canted away of the plane of lattice
at an angle given by ,

$\theta=\frac{\pi}{2}+\frac{1}{2}\tan^{-1}\left(\frac{4D+8(D_{||,s}-D_{||,ns})}{4J+8D_{\perp}-8J^{\prime}}\right),$
(3)

where the spins along one diagonal are canted out of the plane while those
along the other diagonal are canted into the plane.

The magnon bands for the canted-flux state is shown in the Fig.3(b), where the
bands are observed to be degenerate at $\Gamma$, M and X-point in the
Brillouin zone. These degeneracies are protected by symmetries(See Appendix.A
for details).

In the presence of an external longitudinal magnetic field, the canting angles
of the two pairs of spins are no longer identical. The energy of the canted
flux state in a magnetic field is given by,

$\displaystyle\frac{E_{cl}}{NS^{2}}=$ $\displaystyle
J(cos(2\theta_{1})+\cos(2\theta_{2}))+J^{\prime}\cos(\theta_{1})\cos(\theta_{2})$
$\displaystyle+D(\sin(2\theta_{1})-\sin(2\theta_{2}))-8D_{\perp}\sin(\theta_{1})\sin(\theta_{2})$
$\displaystyle+4(D_{||,s}-D_{||,ns})\sin(\theta_{2}-\theta_{1})$
$\displaystyle-\frac{2B}{S}(\cos(\theta_{1})+\cos(\theta_{2})),$ (4)

where $\theta_{1}$ and $\theta_{2}$ are the (different) canting angles angles
made by the spins into and out of the plane of the lattice. The energy is
minimized for $\theta_{1}\neq\theta_{2}$ – the application of magnetic field
not only lifts the non-symmorphic symmetries but also renders the magnetic
symmetry group of the ground state trivial. As a consequence, all the four
bands are gapped out, as shown in Fig.3(c). The magnon Hamiltonian in presence
of all DM-interactions and magnetic field can be found in the Appendix.B.

The lifting of the symmetry-protected degeneracies in the energy bands allows
us to calculate the Berry curvature and Chern number associated with each band
separately in the usual manner(see Appendix.C). The results reveal that the
energy bands acquire topological character for the parameter set chosen in
Fig.3. A representative Berry-curvature distribution for the third band is
shown in Fig.3(d). The Berry curvature is concentrated near the M and M’
points of the Brillouin zone. The existence of the two-fold rotational
symmetry at the center of the diagonal bond is reflected in the Berry-
curvature. Unlike the flux state, the Berry curvature across the entire
Brillouin zone do not cancel and this results in a non-zero Chern number for
three of the four bands (Fig.3(c))(see Appendix.C for details).

Figure 4: (color online) (a) The topological phases of magnon in the parameter
space $D$ and $D_{||,ns}$, with $D_{||,s}=0.1,B_{z}=0.2$. (b)The topological
phases of magnon in the parameter space $D_{||,s}$ and $D$, with
$D_{||,ns}=0.1,B_{z}=0.2$. (c) The topological phases of magnon in the
parameter space $D_{||,s}$ and $B_{z}$, with $D_{||,ns}=0.1,D=0.05$. (d) The
topological phases of magnon in the parameter space $D$ and $B_{z}$, with
$D_{||,s}=0.05,D_{||,ns}=0.1$. For all the plots the residual parameters are
taken as $J=1,J^{\prime}=1.1,D_{\perp}=0.8$. The four digits from left to
right denotes Chern-number from lower band to upper band and the digit with
bar denotes negative Chern-number. The topological phase transition denoted by
the lines red, blue and black boxed-lines. The red, blue and black denotes the
band gap closing between upper, middle and lower pairs of bands respectively.
The thermal conductivity is plotted in Fig.5 along the purple lines in (c) and
(d). The purple lines in (c) corresponds to $D=0.5$. The line in (d)
corresponds to $D=0.2$ .

The interplay between competing Heisenberg and DM interactions, together with
geometric frustration and external magnetic field results in topologically
ordered energy eigenstates. As the relative strengths of the different
competing interactions are varied, the energy levels shift and the bands cross
/ touch in pairs at different points in the Brillouin zone. The accompanying
phase transitions are topological in nature as they are characterized by the
change in Chern number of the pair of bands involved. By identifying the state
of the system with the band topologyKim  _et al._ (2018), we find a wide
variety of topological phases in different parameter regimes, shown in Fig.4.
The four numbers in the figure represent the four Chern numbers from lower to
upper magnon bands. The bar above the number denotes the negative Chern
number. The color of the phase boundaries identify the pair of bands involved
in the transition. Most strikingly, tuning the strength of the different
components of the DMI and applied magnetic field over a small range result in
multitude of topologically distinct set of single magnon bands, even though
the ground state remains the unaltered (canted flux state). While this is
driven by the non-coplanarity of the ground state spin configuration, the
exact mechanism of the change in geometry of the magnon bands, or their
robustness against interaction effects is not clear.Kim  _et al._ (2018).

Topological phase transition occurs due to band reopening after closing at the
the high symmetry points $\Gamma$, X, M and points along line
$\overline{\Gamma\text{M}}$ and $\overline{\Gamma\text{M'}}$. Except for the
$\Gamma$-point, all other k-points in the Brillouin zone can be mapped into
another k-point using two-fold rotational symmetry. Thus, the Chern number of
the bands changes by $\pm 1$, if band touching happens at the $\Gamma$-point.
It is noticeable that this kind of phase transition happen in the upper right
region of Fig.4(b). Otherwise, the the Chern numbers changes by $\pm 2$,
because accidental band touching take places at two points in the Brillouin
zone due to two-fold rotational symmetry of the system. Most of the phase
transition is associated with change in Chern number $\pm 2$ in Fig.4.

####  III.2.2 Thermal Hall conductance and its derivative

The non-trivial topology of magnon bands give rise to thermal Hall effect in
the magnetic system. The expression of reduced thermal Hall conductivity is
given byMatsumoto and Murakami (2011a, b),

$\kappa^{\prime}_{xy}=\frac{\kappa_{xy}\hbar}{k_{B}}=\frac{T^{\prime}}{(2\pi)^{2}}\sum_{n}\int_{BZ}c_{2}(\rho_{n,\mathbb{k}})\Omega^{n}_{xy}(\mathbb{k})d^{2}k,$
(5)

where, $\kappa_{xy}$ is the thermal Hall conductivity, $T^{\prime}$ is the
scaled temperature, $T^{\prime}=k_{B}T$, and
$\rho_{n,\mathbb{k}}=1/(\exp(\epsilon^{n}(\mathbb{k})/T^{\prime})-1)$ is the
Bose-Einstein distribution function with $\epsilon^{n}(\mathbb{k})$ as the
energy of the $n$-th magnon band at $\mathbb{k}$-point in Brillouin zone. The
reduced temperature $T^{\prime}$, magnon energies $\epsilon^{n}(\mathbb{k})$
as well as reduced-thermal Hall conductivity are normalized in unit of $JS$.
In Fig.5(a) and 5(b), the results for the reduced thermal Hall conductivity is
plotted as a function of Magnetic field $B_{z}$ and $D_{||,s}$ respectively,
along the purple lines in Fig.4(b) and Fig.4(d). The different coloured
regions in Fig5(a) and 5(b) denote distinct topological regions along the
purple lines in Fig.4(d) and Fig.4(b) respectively. At the boundary of the
topological regions the band gap closes as shown in the inset of the figures.
Generically, band closing occurs at the Dirac point, but sometime a semi-Dirac
point is encountered. The type of semi-Dirac point at the boundary between
green and purple topological regions in Fig.5(b) is also reported in the
reference Ref.[Wu, 2014].

Figure 5: (color online) (a)-(b) The reduced thermal conductivity as a
function of Magnetic field $B$ and DM-term $D_{||,s}$ along the purple lines
in Fig.4(d) and Fig.4(b) respectively. Inset of the figure shows the band
touching points at the boundaries of the topological region. The first, second
and third letters denote the band touching point, type of band touching, and
the energy wise ordering of bands.For example, ”$\Gamma$M’, TD, UB”, band
touches at $\Gamma$M’ line , the type of band touching is tilted Dirac-type
and Upper bands touch respectively. ”SD” means the Semi-Dirac point and ”UB”
means Upper band touches and ”MB” means middle band touches etc. The magnetic
field is a reduced quantity and connected with experimental magnetic field as
$h_{z}=\frac{BJ}{g\mu_{B}}$, where $g$ is Lande-g factor and $\mu_{B}$ is the
Bohr magneton. (c)-(d) The derivative of the thermal Hall conductivity w.r.t.
magnetic field and DM-term $D_{||,s}$ as a function of Magnetic field $B$ and
DM-term $D_{||,s}$ respectively. The numbers in the figures are Chern numbers
with convention given in Fig.4.

Figs.5(c) and 5(d) present the derivative of reduced thermal Hall conductivity
with respect to magnetic field $B$ and DMI $D_{||,s}$ as a function of $B$ and
the symmetric component of the in-plane DMI $D_{||,s}$. At the boundary
between two distinct topological phases, the derivative in the thermal Hall
conductance has a logarithmic divergence. The origin of the divergence on the
basis of Weyl-point also discussed in Ref.[Kim  _et al._ , 2019]. It is
observed from the figures that the logarithmic divergence is universal and
independent of the type of band touching. We also have been shown analytically
in AppendixD that the nature of divergence is same for tilted Dirac point and
Semi-Dirac point. Furthermore, it also can be seen from the figure that the
peak height of divergence grows faster, if the band touching happens at the
lower pair of bands, due to the larger contribution to the thermal Hall
conductivity from magnons in the lower bands. Finally, the sign of the
divergence is positive (negative) if Chern number of lower band increases
(decreases) at the topological phase transition.

Figure 6: (color online) (a)-(b) The derivative of conductivity as a function
of parameters $B$ and $D_{||,s}$ respectively, for different temperature. The
method of plotting is described in the main text. (c)-(d) The derivative of
conductance as a function of temperature is fitted using Eq.6, at different
$B$ and $D_{||,s}$ respectively. The circled points correspond to the
calculated derivative of thermal Hall conductivity and the straight line
represents the fitted curve. The fitting parameters
$\left\\{A,\epsilon_{0}\right\\}$ for blue, red and black curves in (c) are
$\left\\{4.6886,2.9243\right\\}$, $\left\\{10.4323,2.7539\right\\}$ and
$\left\\{5.5643,2.8666\right\\}$. The fitting parameters
$\left\\{A,\epsilon_{0}\right\\}$ for blue, red and black curves in (d) are
$\left\\{20.0830,4.2135\right\\}$, $\left\\{103.5474,4.1025\right\\}$ and
$\left\\{15.4258,4.2656\right\\}$.

We found that the logarithmic divergence follows simple analytical expression
as a function of temperature given as,

$\frac{\partial\kappa^{\prime}_{xy}}{\partial
p}=A\left[\ln\left(\frac{1+\rho_{0}}{\rho_{0}}\right)\right]^{2}\exp(\frac{\epsilon_{0}}{T})\rho_{0}^{2},$
(6)

where, $\epsilon_{0}$ is the band touching point during topological phase
transition, $\rho_{0}=1/(\exp(\epsilon_{0}/T)-1)$, $A$ is a constant
independent of temperature and proportional to $\log(p)$, $p$ is a parameter
of system(e.g. magnetic field etc.). Moreover, at a temperature, lower
compared with the energy of band touching point, the equation 6 transforms
into,

$\frac{\partial\kappa^{\prime}_{xy}}{\partial
p}=A\epsilon_{0}^{2}\exp(\frac{-\epsilon_{0}}{T}).$ (7)

To demonstrate the validity of Eq.6, we have chosen the topological phase
transition points near $B=0.276$ and $D_{||,s}=0.3937$ of Fig.5(c) and
Fig.5(d) respectively. The numerically calculation Berry-curvature around the
transition point is computationally expensive and in-accurate. So, to
correctly calculate the derivative of thermal Hall-conductivity, first the
thermal Hall conductivity has been calculated near the transition point and
fitted using the expression,

$\kappa^{\prime}_{xy}=m\ln(|p-p_{0}|)+m_{0}+m_{1}(p-p_{0})+m_{2}(p-p_{0})^{2}+m_{3}(p-p_{0})^{3},$
(8)

where $p_{0}$ is the critical point. Then the derivative of the expression has
been plotted and shown in Fig.6(a) and Fig.6(b). The divergent peak and nearby
points increase with the temperature. In Fig.6(c) and 6(d), the derivative of
conductivity is plotted and fitted as a function of temperature using the
Eq.6, considering $A$ and $\epsilon_{0}$ as fitting parameter. The band
touches at $\epsilon_{0}=2.65$ for the phase transition at $B=0.276$ in
Fig.6(a) and at $\epsilon_{0}=4.079$ for the phase transition at
$D_{||,s}=0.3937$ in Fig.6(b). The values of fitting parameter $\epsilon_{0}$
described in Fig.6(c) and Fig.6(d) is quite near the band touching points.

##  IV Material Realizaion

Figure 7: The symmetry of the Shastry-Sutherland lattice corresponding to the
allowed DM-interactions in Fig.1

We have carefully restricted our choice of Hamiltonian parameters to realistic
ranges. Nearly equal Heisenberg exchange interactions on the diagonal and
axial bonds($J\approx J^{\prime}$) in SS-lattice is observed in the rare earth
tetraborides (RB4, R=Er,Tm)Siemensmeyer  _et al._ (2008); Wierschem and
Sengupta (2013); Ye  _et al._ (2017). The nature of DM-interaction for the
study(Fig.1(b)) is based on the symmetry of CuO layers of the canonical
Shastry Sutherland compound, SrCu2(BO3)2 in the low temperature phaseRomhányi
_et al._ (2011); Sparta  _et al._ (2001); Romhányi  _et al._ (2015b), where
the bonds along the two inequivalent diagonals of the unit cell are shifted
out of the plane in opposite directions as shown in Fig.7. Although the nature
and strength of DM-interaction in the native rare earth tetraborides has not
been investigated experimentally, the possibility of forming Van der Waals
heterostructures with heavy metals offers the ability to induce and tune DM
interactions over a extended range. The presence of geometric frustration
reduces the critical $D_{\perp}$ to achieve flux-state in SS-lattice,
facilitating the realization of the canted flux state and associated
topological magnon excitations in these materials. Finally, DM interactions
can also be induced by incident circularly polarized optical waveOwerre
(2017b), further enhancing the possibility .

##  V Conclusion

In conclusion we consider the flux state of Shastry-Sutherland lattice and
showed that in presence of in plane DMI and magnetic field, the system gives
rise to non-trivial topological magnon bands. The canted flux state is a non-
coplanar spin structure. This leads to a various topologically distinct magnon
band structure. Again, we observed the nature of first derivative of thermal
Hall conductance is logarithmic divergent at the topological phase transition,
independent of the type of band touching. We have presented a simple
temperature dependent parametric relation for the thermal Hall conductance,
which might be useful to extract the energy of band touching during
topological phase transition. Finally, we have suggested an experimental
realization of the model studied. In the present work, we have assumed a
dilute gas of magnons without any interaction. At finite temperatures, as the
density of thermally excited magnons increase, effects of interaction gain
importance. Interaction between magnons further re-normalizes the bands and
impart a finite life-time, which in turn can change the topological phase
diagram obtained in this study. The study of the topological magnon bands for
this model in presence of interaction is planned for the future.

Financial support from the Ministry of Education, Singapore, in the form of
grant MOE2016-T2-1-065 is gratefully acknowledged.

##  Appendix A Symmetry protected degeneracy in canted-flux state

The symmetry protection of the band sticking at X-point in the Brillouin zone
in the canted flux state without magnetic field has been explained using
Herring’s methodHerring (1942); BolinLiao (2012).

The symmetry operators which keep the X-point in the Brillouin invariant or
change it by a reciprocal lattice vector are,

$\displaystyle\left\\{M_{2}\middle|\bm{\delta}_{2}\right\\},\;\left\\{M_{1}\middle|\bm{\delta}_{1}\right\\},$
$\displaystyle\left\\{C_{2z}\middle|\bm{\delta}_{y}\right\\},$

where,

$\displaystyle
M_{2}=\tilde{M}_{2}e^{i\pi\hat{S}_{z}}e^{i\pi\hat{S}_{1}},\;M_{1}=\tilde{M}_{1}e^{i\pi\hat{S}_{z}}e^{i\pi\hat{S}_{2}},$
$\displaystyle\;C_{2z}=\tilde{C}_{2z}e^{i\pi\hat{S}_{z}}e^{i\pi\hat{S}_{1}}e^{i\pi\hat{S}_{z}}e^{i\pi\hat{S}_{2}}.$

The symmetry operators $\tilde{M}_{2}$, $\tilde{M}_{1}$ and $\tilde{C}_{2z}$
are mirror reflection along axis-$P_{1}$, reflection along axis-$P_{2}$ and
twofold rotation around z-axis at the sublattice-c respectively(Fig.2(a)).

The set of translational operators makes the invariant subgroup,

$\pazocal{T}_{\mathbb{k}}=\left\\{E\middle|m\bm{\alpha}_{1}+n\bm{\alpha}_{2}\right\\},\quad
m\in\text{even},\;n\in\text{integer},$ (9)

where the translational operator follows the constraint
$exp(i\mathbb{k}\cdot\mathbb{t})=1$ at
X-point($\mathbb{t}=\left(\frac{\pi}{a},0\right)$).

Then the factor group $\pazocal{G}_{\mathbb{k}}/\pazocal{T}_{\mathbb{k}}$ can
be obtained by deriving the coset of invariant subgroup
$\pazocal{T}_{\mathbb{k}}$, where $\pazocal{G}_{\mathbb{k}}$ is the symmetry
operators which keeps the X-point invariant or change it by reciprocal lattice
vector. The factor group $\pazocal{G}_{\mathbb{k}}/\pazocal{T}_{\mathbb{k}}$
is given by,

$\displaystyle
e=\left\\{E\middle|m\bm{\alpha}_{1}+n\bm{\alpha}_{2}\right\\},\quad
e^{\prime}=\left\\{E\middle|p\bm{\alpha}_{1}+q\bm{\alpha}_{2}\right\\}$
$\displaystyle
m_{2}=\left\\{M_{2}\middle|m\bm{\alpha}_{1}+n\bm{\alpha}_{2}+\bm{\delta}_{1}\right\\},\quad
m^{\prime}_{2}=\left\\{M_{2}\middle|p\bm{\alpha}_{1}+q\bm{\alpha}_{2}+\bm{\delta}_{1}\right\\}$
$\displaystyle
m_{1}=\left\\{M_{1}\middle|m\bm{\alpha}_{1}+n\bm{\alpha}_{2}+\bm{\delta}_{2}\right\\},\quad
m^{\prime}_{1}=\left\\{M_{1}\middle|p\bm{\alpha}_{1}+q\bm{\alpha}_{2}+\bm{\delta}_{2}\right\\}$
$\displaystyle
c_{2}=\left\\{C_{2z}\middle|m\bm{\alpha}_{1}+n\bm{\alpha}_{2}+\bm{\delta}_{x}\right\\},\quad
c^{\prime}_{2}=\left\\{C^{\prime}_{2z}\middle|p\bm{\alpha}_{1}+q\bm{\alpha}_{2}+\bm{\delta}_{x}\right\\},$

where $m\in even$, $p\in odd$ and $n,q\in integer$.

Next, we have derived the character table using package named ”GAP”GAP , and
the derived character table is given by,

Table 1: Character Table of unitary subgroup | $e$ | $e^{\prime}$ | $\left\\{m_{2},m^{\prime}_{2}\right\\}$ | $\left\\{m_{1},m^{\prime}_{1}\right\\}$ | $\left\\{c_{2},c^{\prime}_{2}\right\\}$
---|---|---|---|---|---
$\Gamma_{1}$ | 1 | 1 | 1 | 1 | 1
$\Gamma_{2}$ | 1 | 1 | -1 | -1 | 1
$\Gamma_{3}$ | 1 | 1 | 1 | -1 | -1
$\Gamma_{4}$ | 1 | 1 | -1 | 1 | -1
$\Gamma_{5}$ | 2 | -2 | 0 | 0 | 0

The only valid representation is $\Gamma_{5}$, since the translational
operators $e=\left\\{E\middle|p\bm{\alpha}_{1}+q\bm{\alpha}_{2}\right\\}$
should follow the relation,

$\displaystyle exp(i\mathbb{k}\cdot\mathbb{t})=-\pazocal{I},$

where, $\pazocal{I}$ is the identity matrix with a dimension equals to the
dimension of representation.

Because the only valid representation is $\Gamma_{5}$ which is two
dimensional, the bands are doubly degenerate at the X-point.

In the similar manner, it can be shown that the degeneracy at the $\Gamma$ and
M-points are also symmetry protected.

##  Appendix B Magnon Hamiltonian and diagonalization

The quadratic magnon Hamiltonian is given by,

$\displaystyle\pazocal{H}=S$
$\displaystyle\sum_{\mathbb{k}}\left[Q_{\mathbb{k}}(\delta_{1})\hat{a}^{\dagger}_{\mathbb{k}}\hat{c}_{\mathbb{k}}+Q_{\mathbb{k}}(-\delta_{2})^{*}\hat{a}^{\dagger}_{\mathbb{k}}\hat{d}_{\mathbb{k}}\right.$
$\displaystyle+Q_{\mathbb{k}}(\delta_{2})^{*}\hat{b}^{\dagger}_{\mathbb{k}}\hat{c}_{\mathbb{k}}+Q_{\mathbb{k}}(-\delta_{1})\hat{b}^{\dagger}_{\mathbb{k}}\hat{d}_{\mathbb{k}}$
$\displaystyle-\frac{1}{2}(J+\epsilon_{A})e^{-i\mathbb{k}\cdot\mathbb{x}}\hat{a}^{\dagger}_{\mathbb{k}}\hat{b}_{\mathbb{k}}$
$\displaystyle\left.-\frac{1}{2}(J-\epsilon_{B})e^{i\mathbb{k}\cdot\mathbb{y}}\hat{c}^{\dagger}_{\mathbb{k}}\hat{d}_{\mathbb{k}})\right]+\mathtt{h.c.}$
$\displaystyle+S$
$\displaystyle\sum_{\mathbb{k}}\left[P_{\mathbb{k}}(\delta_{1})^{*}\hat{a}_{\mathbb{k}}\hat{c}_{-\mathbb{k}}+P_{\mathbb{k}}(-\delta_{2})\hat{a}_{\mathbb{k}}\hat{d}_{-\mathbb{k}}\right.$
$\displaystyle+P_{\mathbb{k}}(\delta_{2})\hat{b}_{\mathbb{k}}\hat{c}_{-\mathbb{k}}+P_{\mathbb{k}}(-\delta_{1})^{*}\hat{b}_{\mathbb{k}}\hat{d}_{-\mathbb{k}}$
$\displaystyle+\frac{1}{2}\left(J-\epsilon_{A}\right)e^{i\mathbb{k}\cdot\mathbb{x}}\hat{a}_{\mathbb{k}}\hat{b}_{-\mathbb{k}}$
$\displaystyle\left.+\frac{1}{2}\left(J+\epsilon_{B}\right)e^{-i\mathbb{k}\cdot\mathbb{y}}\hat{c}_{\mathbb{k}}\hat{d}_{-\mathbb{k}}\right]+\mathtt{H.c.}$
$\displaystyle+S$
$\displaystyle\sum_{\mathbb{k}}\left[\left(2E-\epsilon_{A}+\frac{B}{S}\cos(\theta_{1})\right)\hat{a}^{\dagger}_{\mathbb{k}}\hat{a}_{\mathbb{k}}\right.$
$\displaystyle+\left(2E-\epsilon_{A}+\frac{B}{S}\cos(\theta_{1})\right)\hat{b}^{\dagger}_{\mathbb{k}}\hat{b}_{\mathbb{k}}$
$\displaystyle+\left(2E+\epsilon_{B}+\frac{B}{S}\cos(\theta_{2})\right)\hat{c}^{\dagger}_{\mathbb{k}}\hat{c}_{\mathbb{k}}$
$\displaystyle\left.+\left(2E+\epsilon_{B}+\frac{B}{S}\cos(\theta_{2})\right)\hat{d}^{\dagger}_{\mathbb{k}}\hat{d}_{\mathbb{k}}\right],$
(10)

where,

$\displaystyle E$
$\displaystyle=-2J^{\prime}C^{X}_{\theta}+2D_{\perp}S^{X}_{\theta}-D_{||,ns}(\zeta_{\theta_{1}}-\zeta_{\theta_{2}})$
$\displaystyle\quad-D_{||,s}(\zeta_{\theta_{2}}-\zeta_{\theta_{1}})$
$\displaystyle\epsilon_{A}$
$\displaystyle=J\cos(2\theta_{1})+D\sin(2\theta_{1})$
$\displaystyle\epsilon_{B}$
$\displaystyle=-J\cos(2\theta_{2})+D\sin(2\theta_{2})$ $\displaystyle
Q_{\mathbb{k}}(\delta)$
$\displaystyle=\lambda_{J^{\prime}}^{+}(\mathbb{k},\delta)-\lambda_{D_{\perp}}^{+}(\mathbb{k},\delta)$
$\displaystyle-
D_{||,ns}(\xi_{\theta_{1}}e^{i\mathbb{k}\cdot\mathbb{\delta}}-\xi_{\theta_{2}}e^{-i\mathbb{k}\cdot\mathbb{\delta}})$
$\displaystyle+D_{||,s}(\xi_{\theta_{1}}e^{-i\mathbb{k}\cdot\mathbb{\delta}}-\xi_{\theta_{2}}e^{i\mathbb{k}\cdot\mathbb{\delta}})$
$\displaystyle P_{\mathbb{k}}(\delta)$
$\displaystyle=\lambda_{J^{\prime}}^{-}(\mathbb{k},\delta)-\lambda_{D_{\perp}}^{-}(\mathbb{k},\delta)$
$\displaystyle-
D_{||,ns}(\xi_{\theta_{1}}^{*}e^{i\mathbb{k}\cdot\mathbb{\delta}}-\xi_{\theta_{2}}e^{-i\mathbb{k}\cdot\mathbb{\delta}})$
$\displaystyle+D_{||,s}(\xi_{\theta_{1}}^{*}e^{-i\mathbb{k}\cdot\mathbb{\delta}}-\xi_{\theta_{2}}e^{i\mathbb{k}\cdot\mathbb{\delta}})$
$\displaystyle\lambda_{J^{\prime}}^{\pm}(\mathbb{k},\delta)$
$\displaystyle=J^{\prime}(iC^{\pm}_{\theta}+S^{X}_{\theta})\cos(\mathbb{k}\cdot\delta)$
$\displaystyle\lambda_{D_{\perp}}^{\pm}$
$\displaystyle=D_{\perp}(C^{X}_{\theta}\pm 1)\cos(\mathbb{k}\cdot\delta)$
$\displaystyle\xi_{\theta_{i}}$
$\displaystyle=\frac{i}{2}\sin(\theta_{i})+\frac{1}{2}\cos(\theta_{i})\sin(\theta_{\bar{i}})$
$\displaystyle\zeta_{\theta_{i}}$
$\displaystyle=\sin(\theta_{i})\cos(\theta_{\bar{i}}),\;\;(i,\bar{i})=(1,2)\;\text{or}\;(2,1)$
$\displaystyle C^{X}_{\theta}$
$\displaystyle=\cos(\theta_{1})\cos(\theta_{2})$ $\displaystyle
C^{\pm}_{\theta}$ $\displaystyle=\cos(\theta_{1})\pm\cos(\theta_{2})$
$\displaystyle S^{X}_{\theta}$
$\displaystyle=\sin(\theta_{1})\sin(\theta_{2})$

The Hamiltonian can be re-written in the form of Bogoliubov de Gennes matrix,

$\pazocal{H}=\sum_{\mathbb{k}}\Psi^{\dagger}(\mathbb{k})\pazocal{H}(\mathbb{k})\Psi(\mathbb{k}),$
(13)

where
$\Psi^{\dagger}(\mathbb{k})=(\hat{a}^{\dagger}_{\mathbb{k}},\hat{b}^{\dagger}_{\mathbb{k}},\hat{c}^{\dagger}_{\mathbb{k}},\hat{d}^{\dagger}_{\mathbb{k}},\hat{a}_{-\mathbb{k}},\hat{b}_{-\mathbb{k}},\hat{c}_{-\mathbb{k}},\hat{d}_{-\mathbb{k}})^{T}$.
The matrix $\pazocal{H}(\mathbb{k})$ is diagonalized using the para-unitary
transformation,

$T^{\dagger}(\mathbb{k})\pazocal{H}(\mathbb{k})T(\mathbb{k})=\epsilon(\mathbb{k}),$
(14)

where,

$\epsilon(\mathbb{k})=\begin{pmatrix}E_{1}(\mathbb{k})&&&&\\\ &...&&&&\\\
&&E_{N}(\mathbb{k})&&&\\\ &&&E_{1}(-\mathbb{k})&&\\\ &&&&...&\\\
&&&&&E_{N}(-\mathbb{k})\end{pmatrix}.$ (15)

The para-unitary transformation preserves the Bosonic commutation relation of
the eigenvectors obtained after diagonalization.

The para-unitary matrices follow the relation,

$T^{\dagger}(\mathbb{k})\sigma_{3}T(\mathbb{k})=T(\mathbb{k})\sigma_{3}T^{\dagger}(\mathbb{k})=\sigma_{3},$
(16)

where $\sigma_{3}$ is the 3rd component of Pauli matrix. Using Eq.16, we can
re-write Eq.14 as,

$\sigma_{3}\pazocal{H}(\mathbb{k})T(\mathbb{k})=T(\mathbb{k})\sigma_{3}\epsilon(\mathbb{k})$
(17)

So the problem of para-unitary transformation become the problem of
diagonalization of the matrix $\sigma_{3}\pazocal{H}(\mathbb{k})$.

##  Appendix C Calculation of Chern number and Berry curvature

The presence of non-zero Chern numbers of the bands denote non-trivial
topology of the magnon bands. The expression of the Chern number for n-th band
is given by,

$C_{n}=\frac{1}{2\pi}\int_{\text{BZ}}\Omega^{n}_{xy}(\mathbb{k})dk_{x}dk_{y},$
(18)

where the Berry curvatureMatsumoto  _et al._ (2014),

$\Omega^{n}_{xy}(\mathbb{k})=i\epsilon_{\mu\nu}\left[\sigma_{3}\frac{\partial
T^{\dagger}_{\mathbb{k}}}{\partial k_{\mu}}\sigma_{3}\frac{\partial
T_{\mathbb{k}}}{\partial
k_{\nu}}\right]_{nn}\qquad\left(n=1,2,3,...,2N\right),$ (19)

Where $N$ is the number of sublatices and $T_{\mathbb{k}}$ is the para-unitary
matrix which diagonalizes the Bogoliubov Hamiltonian(See Appendix.B). Each
column of the para-unitary matrix $T_{\mathbb{k}}$ corresponds to the wave-
functions at $\mathbb{k}$. The Eq.19 is not suitable for numerical
calculation, since numerically Gauge degrees of freedom of the wave-function
is not well controlled and so the derivative of the wave-function is ill
defined. The Berry curvature in Eq.19 can be re-structured as(See Appendix.C),

$\displaystyle\Omega^{n}_{xy}(\mathbb{k})=\sum_{m\neq
n}i(\sigma_{3})_{nn}(\sigma_{3})_{mm}$ $\displaystyle\left[\frac{\left\langle
u^{n}(\mathbb{k})\middle|\frac{\partial\pazocal{H}_{\mathbb{k}}}{\partial
k_{x}}\middle|u^{m}(\mathbb{k})\right\rangle\left\langle
u^{m}(\mathbb{k})\middle|\frac{\partial\pazocal{H}_{\mathbb{k}}}{\partial
k_{y}}\middle|u^{n}(\mathbb{k})\right\rangle-(k_{x}\leftrightarrow
k_{y})}{((\sigma_{3}\epsilon_{\mathbb{k}})_{nn}-(\sigma_{3}\epsilon_{\mathbb{k}})_{mm})^{2}}\right],$
(20)

where, $\left|u^{m}(\mathbb{k})\right\rangle$ is the m-th column of the para-
unitary matrix $T(\mathbb{k})$ and $\left\langle u^{m}(\mathbb{k})\right|$ is
the m-th row of the para-unitary matrix $T^{\dagger}(\mathbb{k})$.

To calculate the non-Abelian Berry curvature, we followed the method described
in Ref.Matsuoka  _et al._ , 2005. For the sake of completeness, here we
describe the calculation briefly. The expression of Chern number is given by,

$\tilde{F}_{12}(k_{l})=\frac{1}{i}\ln(U_{1}(k_{l})U_{2}(k_{l}+\hat{1})U_{1}(k_{l}+\hat{2})^{-1}U_{2}(k_{l})^{-1}),$
(21)

where, $k_{l}$ is the discrete k-point in the Brillouin zone and if $N_{i}$ is
the number of lattice points in the $i$-th direction(where $i=x$ or $y$) then,
$l=1,2,...,N_{1}N_{2}$. Again $k_{l}+\hat{i}$ denotes the k-point next to
$k_{l}$ along the $i$-th direction. Again,

$U_{\mu}(k_{l})=\text{det}(\psi^{\dagger}(k_{l})\psi(k_{l}+\hat{\mu}))/N_{\mu}(k_{l}),$
(22)

where,
$N_{\mu}(k_{l})=|\text{det}(\psi^{\dagger}(k_{l})\psi(k_{l}+\hat{\mu}))|$ and
$\psi(k_{l})=(T_{\mathbb{k}})_{mn}$. The indices $m$ and $n$ denote the band-
indices for which the non-Abelian Berry-curvature is calculated.

##  Appendix D Logarithmic divergence in derivative of thermal Hall
conductivity near phase transition

###  D.1 Tilted Dirac point or Generalized Weyl point

The Hamiltonian corresponds to the tilted Dirac point or generalized Weyl
point is given by,

$\hat{H}=w^{0}_{x}k_{x}\mathbb{I}+w^{0}_{y}k_{y}\mathbb{I}+w_{x}k_{x}\sigma_{x}+w_{y}k_{y}\sigma_{y}+p\sigma_{z}+\epsilon\mathbb{I},$
(23)

where $k_{x}$, $k_{y}$ are the momentum with respect to the band touching
point. $\sigma_{i}$ ($i=x,y,z$) are the Pauli’s matrices. $p$ is a
perturbation to open the gap. $\epsilon$ denotes the energy of band touching
point. The Hamiltonian is general Hamiltonian for any linear dispersions.

The energies correspond to the Hamiltonian,

$\displaystyle E^{\pm}(\mathbb{k})$
$\displaystyle=\epsilon+k_{x}w^{0}_{x}+k_{y}w^{0}_{y}\pm\sqrt{p^{2}+w^{2}_{x}k^{2}_{x}+w^{2}_{y}k^{2}_{y}}$
$\displaystyle=\epsilon_{0\mathbb{k}}\pm\omega_{\mathbb{k}},$ (24)

where $\epsilon_{0\mathbb{k}}=\epsilon_{0}+k_{x}w^{0}_{x}+k_{y}w^{0}_{y}$ and
$\omega_{\mathbb{k}}=\sqrt{p^{2}+w^{2}_{x}k^{2}_{x}+w^{2}_{y}k^{2}_{y}}$ the
corresponding eigenvectors are,

$x^{\pm}=\frac{1}{N_{\pm}}\begin{pmatrix}\frac{p\pm\sqrt{p^{2}+w^{2}_{x}k^{2}_{x}+w^{2}_{y}k^{2}_{y}}}{w_{x}k_{x}+iw_{y}k_{y}}\\\
1\end{pmatrix},$ (25)

where $N_{\pm}$ are the normalization constant.

The expression of the Berry curvature of the lower band is given by,

$\displaystyle\Omega_{1}$
$\displaystyle=\frac{1}{(E^{+}-E^{-})^{2}}\text{Im}\left(\left\langle
x_{1}\middle|\frac{\partial\pazocal{H}}{\partial
k_{x}}\middle|x_{2}\right\rangle\left\langle
x_{2}\middle|\frac{\partial\pazocal{H}}{\partial
k_{y}}\middle|x_{1}\right\rangle\right),$ (26)

where Im$(z)$ denotes imaginary part of $z$.

Using Eq.23, Eq.24 and Eq.25,

$\left\langle x_{1}\middle|\frac{\partial\pazocal{H}}{\partial
k_{x}}\middle|x_{2}\right\rangle=\frac{2w_{x}}{\sqrt{N_{1}N_{2}}}\left[\frac{w_{x}k_{x}p-iw_{y}k_{y}\sqrt{p^{2}w^{2}_{x}k^{2}_{x}+w^{2}_{y}k^{2}_{y}}}{w^{2}_{x}k^{2}_{x}+w^{2}_{y}k^{2}_{y}}\right]$
(27)

Similarly,

$\left\langle x_{2}\middle|\frac{\partial\pazocal{H}}{\partial
k_{y}}\middle|x_{1}\right\rangle=\frac{2w_{y}}{\sqrt{N_{1}N_{2}}}\left[\frac{w_{y}k_{y}p-iw_{x}k_{x}\sqrt{p^{2}w^{2}_{x}k^{2}_{x}+w^{2}_{y}k^{2}_{y}}}{w^{2}_{x}k^{2}_{x}+w^{2}_{y}k^{2}_{y}}\right]$
(28)

Using Eq.27 and Eq.28 in Eq. 26, the Berry curvature of the lower band,

$\Omega_{1}=-\frac{2w_{x}w_{y}p}{(p^{2}+w^{2}_{x}k^{2}_{x}+w^{2}_{y}k^{2}_{y})^{3/2}}$
(29)

Similarly, for the upper band,

$\Omega_{2}=+\frac{2w_{x}w_{y}p}{(p^{2}+w^{2}_{x}k^{2}_{x}+w^{2}_{y}k^{2}_{y})^{3/2}}$
(30)

The Thermal Hall conductivity expression is given by,

$\kappa_{xy}=\frac{k_{B}^{2}T}{(2\pi)^{2}\hbar}\sum_{n}\int_{BZ}c_{2}(\rho_{n\mathbb{k}})\Omega_{n}(\mathbb{k})d^{2}k,$
(31)

$\displaystyle\sum_{n}\Omega_{n}(\mathbb{k})c_{2}(\rho_{n\mathbb{k}})$
$\displaystyle=\frac{2w_{x}w_{y}p}{(p^{2}+w^{2}_{x}k^{2}_{x}+w^{2}_{y}k^{2}_{y})^{3/2}}\left[c_{2}(\rho_{2\mathbb{k}})-c_{2}(\rho_{1\mathbb{k}})\right]$
$\displaystyle\propto\frac{2w_{x}w_{y}p}{(p^{2}+w^{2}_{x}k^{2}_{x}+w^{2}_{y}k^{2}_{y})^{3/2}}\omega_{\mathbb{k}}\quad$
$\displaystyle\left[where,c_{2}(\epsilon_{0\mathbb{k}}+\right.$
$\displaystyle\left.\omega_{\mathbb{k}})-c_{2}(\epsilon_{0\mathbb{k}}-\omega_{\mathbb{k}})\propto\omega_{\mathbb{k}}\text{for
small $p$ and $\mathbb{k}$ }\right]$
$\displaystyle\propto\frac{p}{(p^{2}+w^{2}_{x}k^{2}_{x}+w^{2}_{y}k^{2}_{y})}$
(32)

From Eq.31, integrating around the band touching points
$k=\sqrt{w^{2}_{x}k^{2}_{x}+w^{2}_{y}k^{2}_{y}}<k_{c}$, we get,

$\displaystyle\frac{\partial\kappa_{xy}}{\partial p}$
$\displaystyle\propto\frac{\partial}{\partial
p}\int_{k<k_{c}}\frac{p}{(p^{2}+w^{2}_{x}k^{2}_{x}+w^{2}_{y}k^{2}_{y})}d^{2}k$
$\displaystyle\propto\int_{k<k_{c}}\frac{2\pi
kdk}{(p^{2}+k^{2})}-\int_{k<k_{c}}\frac{4p^{2}\pi kdk}{(p^{2}+k^{2})^{2}}$
$\displaystyle\propto\pi\ln\left(\frac{p^{2}+k^{2}_{c}}{p^{2}}\right)+2p^{2}\pi\left[\frac{1}{p^{2}+k^{2}_{c}}-\frac{1}{p^{2}}\right]$
$\displaystyle\propto\ln(|p|)\quad\left[\text{Near phase transtion $k_{c}\gg
p$}\right]$ (33)

###  D.2 Semi-Dirac point

The dispersion of semi-Dirac point,

$E^{\pm}(\mathbb{k})=\pm\sqrt{(w^{2}_{x}k^{2}_{x})^{2}+w^{2}_{y}q^{2}_{y}}$
(34)

Out of many possibilities, two different possible Hamiltonians are,

$\pazocal{H}_{1}=\begin{pmatrix}0&w_{x}q_{x}+iw_{y}q_{y}\\\
w_{x}q_{x}-iw_{y}q_{y}&0\end{pmatrix}$

$\pazocal{H}_{2}=\begin{pmatrix}w_{x}q_{x}&iw_{y}q_{y}\\\
-iw_{y}q_{y}&-w_{x}q_{x}\end{pmatrix},$ (35)

The Berry curvatures of the both Hamiltonian is given by,

$\Omega^{\pm}_{1}=\pm\frac{4pq_{x}w_{x}w_{y}}{(p^{2}+w^{4}_{x}q^{4}_{x}+w^{2}_{y}q^{2}_{y})^{3/2}},\Omega^{\pm}_{2}=0,$
(36)

Thus, the second Hamiltonian in Eq.35 does not produce any Berry curvature and
so topological phase transition can not happen with this kind of Hamiltonian.
Similarly, it can be proved that,

$\frac{\partial\kappa_{xy}}{\partial p}\propto\ln(|p|)$ (37)

Particularly, in our study, we got semi-Dirac points with dispersion,

$E^{\pm}(\mathbb{k})=w_{0}^{2}k_{y}^{2}\pm\sqrt{w_{x}^{2}k_{x}^{2}+w_{y}^{2}k_{y}^{4}}$
(38)

Again particularly for phase transition point, near $D_{||,s}=0.4$ in Fig.5(d)
is associated with Semi-Dirac dispersion with $w_{0}=w_{y}$, which is also
encountered in referenceWu (2014).

##  Appendix E Temperature dependence of derivative of thermal Hall
conductivity near phase transition

The reduced thermal Hall conductivity is given as,

$\kappa^{\prime}_{xy}=\frac{\kappa_{xy}\hbar}{k_{B}}=\frac{T^{\prime}}{(2\pi)^{2}}\sum_{n}\int_{BZ}c_{2}(\rho_{n,\mathbb{k}})\Omega^{n}_{xy}(\mathbb{k})d^{2}k.$
(39)

Near a topological phase transition, the derivative of the of Thermal hall
conductivity contributes to a logarithmic divergence. So, the main
contribution of derivative of thermal Hall conductance comes from the band
touching point at energy $\epsilon_{0}$. Near this point at phase transition
the Berry curvature is equal and opposite for the two
bands$(\Omega^{1}_{xy}=-\Omega^{2}_{xy})$. Near the band touching point,

$\displaystyle\frac{\partial\kappa^{\prime}_{xy}}{\partial p}$
$\displaystyle=\frac{T^{\prime}}{2\pi^{2}}\frac{\partial}{\partial
p}\left[\Omega^{1}_{xy}(\mathbb{k})\left(c_{2}(\rho(\epsilon_{0}+\epsilon_{\mathbb{k}})-c_{2}(\rho(\epsilon_{0}-\epsilon_{\mathbb{k}}))\right)\right]$
$\displaystyle=\frac{T^{\prime}}{2\pi^{2}}\frac{\partial}{\partial
p}\left[\Omega^{1}_{xy}(\mathbb{k})\;\left(2\left.\frac{dc_{2}}{d\rho}\right|_{\rho_{0}}\left.\frac{d\rho}{d\epsilon}\right|_{\epsilon_{0}}\right)\epsilon_{\mathbb{k}}\right]$
$\displaystyle=A\rho_{0}^{2}\exp(\frac{\epsilon_{0}}{T^{\prime}})\left[\ln\left(\frac{1+\rho_{0}}{\rho_{0}}\right)\right]^{2},\;$
(40)

where, $\rho_{0}=1/(\exp(\frac{\epsilon_{0}}{T^{\prime}})-1)$ and $A$ is a
temperature independent parameter.

## References

  * Kane and Mele (2005) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005).
  * Zhang  _et al._ (2005) Y. Zhang, Y.-W. Tan, H. L. Stormer, and P. Kim, Nature 438, 201–204 (2005).
  * Owerre (2016a) S. A. Owerre, Journal of Physics Condensed Matter 28, 386001 (2016a), arXiv:1602.06772 [cond-mat.str-el] .
  * Karzig  _et al._ (2015) T. Karzig, C.-E. Bardyn, N. H. Lindner, and G. Refael, Phys. Rev. X 5, 031001 (2015).
  * Kim  _et al._ (2016a) S. K. Kim, H. Ochoa, R. Zarzuela, and Y. Tserkovnyak, Physical Review Letters 117 (2016a), 10.1103/physrevlett.117.227201.
  * Zhang  _et al._ (2018a) Y. Zhang, S. Okamoto, and D. Xiao, Physical Review B 98 (2018a), 10.1103/physrevb.98.035424.
  * Shuichi and Akihiro (2017) M. Shuichi and O. Akihiro (2017).
  * Cheng  _et al._ (2016a) R. Cheng, S. Okamoto, and D. Xiao, Phys. Rev. Lett. 117, 217202 (2016a).
  * Zyuzin and Kovalev (2016a) V. A. Zyuzin and A. A. Kovalev, Phys. Rev. Lett. 117, 217203 (2016a).
  * Shiomi  _et al._ (2017) Y. Shiomi, R. Takashima, and E. Saitoh, Phys. Rev. B 96, 134425 (2017).
  * Kim  _et al._ (2016b) S. K. Kim, H. Ochoa, R. Zarzuela, and Y. Tserkovnyak, Phys. Rev. Lett. 117, 227201 (2016b).
  * Jin  _et al._ (2015) H. Jin, S. R. Boona, Z. Yang, R. C. Myers, and J. P. Heremans, Phys. Rev. B 92, 054436 (2015).
  * Fert  _et al._ (2013) A. Fert, V. Cros, and J. Sampaio, Nature nanotechnology 8 3, 152 (2013).
  * Jonietz  _et al._ (2010) F. Jonietz, S. Mühlbauer, C. Pfleiderer, A. Neubauer, W. Münzer, A. Bauer, T. Adams, R. Georgii, P. Böni, R. Duine, K. Everschor-Sitte, M. Garst, and A. Rosch, Science (New York, N.Y.) 330, 1648 (2010).
  * Lee  _et al._ (2008) S. Lee, A. Pirogov, M. Kang, K.-H. Jang, M. Yonemura, T. Kamiyama, S.-W. Cheong, F. Gozzo, N. Shin, H. Kimura, Y. Noda, and J.-G. Park, Nature 451, 805 (2008).
  * Kim  _et al._ (2018) K.-S. Kim, K. H. Lee, S. B. Chung, and J.-G. Park, arXiv e-prints , arXiv:1811.11201 (2018), arXiv:1811.11201 [cond-mat.mes-hall] .
  * Shahzad and Sengupta (2017) M. Shahzad and P. Sengupta, Phys. Rev. B 96, 224401 (2017), arXiv:1701.04297 [cond-mat.str-el] .
  * Tomiyasu  _et al._ (2012) K. Tomiyasu, K. Matsuhira, K. Iwasa, M. Watahiki, S. Takagi, M. Wakeshima, Y. Hinatsu, M. Yokoyama, K. Ohoyama, and K. Yamada, Journal of the Physical Society of Japan 81, 034709 (2012), https://doi.org/10.1143/JPSJ.81.034709 .
  * Sagayama  _et al._ (2013) H. Sagayama, D. Uematsu, T. Arima, K. Sugimoto, J. J. Ishikawa, E. O’Farrell, and S. Nakatsuji, Phys. Rev. B 87, 100403 (2013).
  * Disseler (2014) S. M. Disseler, Phys. Rev. B 89, 140413 (2014).
  * Owerre (2016b) S. A. Owerre, Journal of Physics Condensed Matter 28, 386001 (2016b), arXiv:1602.06772 [cond-mat.str-el] .
  * Zyuzin and Kovalev (2016b) V. A. Zyuzin and A. A. Kovalev, Phys. Rev. Lett. 117, 217203 (2016b).
  * Cheng  _et al._ (2016b) R. Cheng, S. Okamoto, and D. Xiao, Phys. Rev. Lett. 117, 217202 (2016b).
  * Owerre (2017a) S. A. Owerre, Journal of Applied Physics 121, 223904 (2017a), arXiv:1608.00545 [cond-mat.mes-hall] .
  * Lee  _et al._ (2018) K. H. Lee, S. B. Chung, K. Park, and J.-G. Park, Phys. Rev. B 97, 180401 (2018).
  * Owerre (2017b) S. A. Owerre, Journal of Physics Condensed Matter 29, 03LT01 (2017b), arXiv:1608.08605 [cond-mat.str-el] .
  * Kim  _et al._ (2019) K.-S. Kim, K. H. Lee, S. B. Chung, and J.-G. Park, Phys. Rev. B 100, 064412 (2019).
  * Laurell and Fiete (2018) P. Laurell and G. A. Fiete, Phys. Rev. B 98, 094419 (2018), arXiv:1804.09783 [cond-mat.str-el] .
  * Miyahara and Ueda (1999) S. Miyahara and K. Ueda, Phys. Rev. Lett. 82, 3701 (1999).
  * Koga and Kawakami (2000) A. Koga and N. Kawakami, Phys. Rev. Lett. 84, 4461 (2000).
  * Läuchli  _et al._ (2002) A. Läuchli, S. Wessel, and M. Sigrist, Phys. Rev. B 66, 014401 (2002).
  * Kageyama  _et al._ (1999) H. Kageyama, K. Yoshimura, R. Stern, N. V. Mushnikov, K. Onizuka, M. Kato, K. Kosuge, C. P. Slichter, T. Goto, and Y. Ueda, Phys. Rev. Lett. 82, 3168 (1999).
  * Siemensmeyer  _et al._ (2008) K. Siemensmeyer, E. Wulf, H.-J. Mikeska, K. Flachbart, S. Gabáni, S. Mat’aš, P. Priputen, A. Efdokimova, and N. Shitsevalova, Phys. Rev. Lett. 101, 177201 (2008).
  * Wierschem and Sengupta (2013) K. Wierschem and P. Sengupta, Phys. Rev. Lett. 110, 207207 (2013).
  * Ye  _et al._ (2017) L. Ye, T. Suzuki, and J. G. Checkelsky, Physical Review B 95, 174405 (2017), arXiv:1606.04003 [cond-mat.mtrl-sci] .
  * Shahzad and Sengupta (2017) M. Shahzad and P. Sengupta, Journal of Physics: Condensed Matter 29, 305802 (2017).
  * Romhányi  _et al._ (2015a) J. Romhányi, K. Penc, and R. Ganesh, Nature Communications 6, 6805 (2015a), arXiv:1406.1163 [cond-mat.str-el] .
  * Malki and Uhrig (2019) M. Malki and G. S. Uhrig, Phys. Rev. B 99, 174412 (2019), arXiv:1811.08981 [cond-mat.str-el] .
  * Zhang  _et al._ (2010) L. Zhang, J. Ren, J.-S. Wang, and B. Li, Phys. Rev. Lett. 105, 225901 (2010).
  * McClarty  _et al._ (2017) P. A. McClarty, F. Krüger, T. Guidi, S. F. Parker, K. Refson, A. W. Parker, D. Prabhakaran, and R. Coldea, Nature Physics 13, 736 (2017), arXiv:1609.01922 [cond-mat.str-el] .
  * Malki and Schmidt (2017) M. Malki and K. P. Schmidt, Phys. Rev. B 95, 195137 (2017).
  * Owerre (2017a) S. A. Owerre, Journal of Applied Physics 121, 223904 (2017a).
  * Zhang  _et al._ (2018b) Y. Zhang, S. Okamoto, and D. Xiao, Physical Review B 98 (2018b), 10.1103/physrevb.98.035424.
  * Lee  _et al._ (2018) K. H. Lee, S. B. Chung, K. Park, and J.-G. Park, Phys. Rev. B 97, 180401 (2018), arXiv:1712.09801 [cond-mat.mes-hall] .
  * Young and Kane (2015) S. M. Young and C. L. Kane, Phys. Rev. Lett. 115, 126803 (2015).
  * Matsuoka  _et al._ (2005) E. Matsuoka, K. Hayashi, A. Ikeda, K. Tanaka, T. Takabatake, and M. Matsumura, Journal of the Physical Society of Japan 74, 1382 (2005), arXiv:cond-mat/0503172 [cond-mat.mes-hall] .
  * Matsumoto and Murakami (2011a) R. Matsumoto and S. Murakami, Phys. Rev. B 84, 184406 (2011a).
  * Matsumoto and Murakami (2011b) R. Matsumoto and S. Murakami, Phys. Rev. Lett. 106, 197202 (2011b).
  * Wu (2014) Y. Wu, Opt. Express 22, 1906 (2014).
  * Romhányi  _et al._ (2011) J. Romhányi, K. Totsuka, and K. Penc, Phys. Rev. B 83, 024413 (2011).
  * Sparta  _et al._ (2001) K. Sparta, G. J. Redhammer, P. Roussel, G. Heger, G. Roth, P. Lemmens, A. Ionescu, M. Grove, G. Güntherodt, F. Hüning, H. Lueken, H. Kageyama, K. Onizuka, and Y. Ueda, European Physical Journal B 19, 507 (2001), arXiv:cond-mat/0012383 [cond-mat.str-el] .
  * Romhányi  _et al._ (2015b) J. Romhányi, K. Penc, and R. Ganesh, Nature Communications 6, 6805 (2015b), arXiv:1406.1163 [cond-mat.str-el] .
  * Owerre (2017b) S. A. Owerre, Journal of Physics Communications 1, 021002 (2017b).
  * Herring (1942) C. Herring, Journal of the Franklin Institute 233, 525 (1942).
  * BolinLiao (2012) BolinLiao, (2012).
  * (56) GAP,  _GAP – Groups, Algorithms, and Programming, Version 4.10.1_, The GAP Group (2019).
  * Matsumoto  _et al._ (2014) R. Matsumoto, R. Shindou, and S. Murakami, Phys. Rev. B 89, 054420 (2014).