# Spinon quantum spin Hall state in the kagome antiferromagnet with a
Dzyaloshinskii-Moriya interaction

Li-Wei He  National Laboratory of Solid State Microstructures and School of
Physics, Nanjing University, Nanjing 210093, China  Jian-Xin Li
<EMAIL_ADDRESS>National Laboratory of Solid State Microstructures and School
of Physics, Nanjing University, Nanjing 210093, China  Collaborative
Innovation Center of Advanced Microstructures, Nanjing University, Nanjing
210093, China

###### Abstract

We investigate the spin-$\frac{1}{2}$ antiferromagnetic Heisenberg model with
a Dzyaloshinskii-Moriya interaction on kagome lattice, making use of the
variational Monte Carlo technique. An exotic quantum spin state is found to
arise from a melting of the $\bm{Q}=0$ long-range magnetic order by a
topological transition, when a small anisotropic third nearest-neighbor
antiferromagnetic Heisenberg interaction is turned on. This novel state is a
gapped quantum spin liquid, characterized by a topological order with ground-
state degeneracy $n_{g}=4$ and topological entanglement entropy $\gamma=\ln
2$, suggesting it is an Abelian topological phase. Furthermore, the Chern
numbers of the spin-up (-down) spinon occupied bands of this state are
$C_{\uparrow\downarrow}=\pm 1$, respectively. From this perspective, this
state is also a time-reversal symmetric (total Chern number $C_{total}=0$)
topological insulator with spinons as the chiral edge states, which carry
opposite spin and move in the opposite direction. It is analogous to the
quantum spin Hall state but the spin current is carried by deconfined spinons
in a quantum spin liquid, so is dubbed as the spinon quantum spin Hall state.

††preprint: APS/123-QED

##  I introduction

One of the exotic and intriguing phase receiving extensive attention and
research in condensed matter physics is the quantum spin liquid (QSL) [1],
which is magnetic disordered and has the fractionalized elementary excitation
called spinon with spin-$\frac{1}{2}$. One class of this phase is the gapped
QSL, which is a crucial representative of the topological orders with the
long-range many-body quantum entanglement. Among the gapped QSLs, the chiral
spin liquid (CSL) [2] breaks the time-reversal (TR) symmetry and usually has
nontrivial Chern number to characterize itself. In a CSL, the spinons with up
and down spins usually couple the same gauge field and have the same Chern
number. For those gapped QSLs with TR symmetry, the ground-state degeneracy
(GSD) [3, 4] other than the Chern number is a good quantity to describe the
global topological property. Furthermore, another important quantity, the
topological entanglement entropy (TEE) [5, 6] related to the quantum dimension
of topological excitaitons, is very useful for both chiral and achiral
topological orders.

Another exotic and extensively studied phase is the quantum spin Hall (QSH)
state [7, 8], which is realized in the topological insulator [9] characterized
by an insulating bulk gap and gapless edge or surface states (the bulk-edge
correspondence of topological system) protected by TR symmetry. In the QSH
state, the electrons with opposite spins move along the opposite direction on
a given edge [7, 10, 11]. Consequently, the two states on the edge possess the
spin Chern number with opposite signs so that the QSH state is characterized
by a Z2 topological number.

Pictorially, both the gapped QSL with TR symmetry and the QSH state are
insulators and topological nontrivial phases. However, the conventional
topological insulator is a non (or weak)-interacting single particle state
whose quasiparticles are electrons in the framework of Landau Fermi liquid.
While, the QSL is a strong interacting Mott insulator without conventional
Landau quasiparticles but the fractionalized excitations such as spinons.
Hence, thus far, these two states are studied separatively from each other.
Thus, an issue arises as to if there exists a possible exotic quantum state
which is a gapped QSL with TR symmetry meanwhile the quasiparticles (spinons)
exhibit similar QSH effect.

To search for QSLs, the kagome antiferromagnet is an appealing platform [12,
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26], in view of its strong
geometric frustration and the resulting strong quantum fluctuations. Many
novel many-body quantum states have been explored in this lattice recently
[27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38]. It is generally believed
that the Heisenberg model with the nearest-neighbor (NN) antiferromagnetic
(AFM) interaction $J_{1}$ on a kagome lattice hosts a QSL ground state.
However, it is controversial that this QSL is gapped [12] or not [13, 14, 15].
Further, the long-range AFM Heisenberg spin interactions beyond the NN term
are always possible and will affect the properties of ground states. It is
reported that the second NN AFM interaction $J_{2}$ favors a $\bm{Q}=0$
magnetic order [16, 14, 17, 18]. However, the interplay between the $J_{2}$
term and one of the third NN AFM couplings which is across the diagonals of
hexagons $J_{d}$ can induce a CSL [19, 20, 17] with spontaneous TR symmetry
breaking. And even, the CSL can arise in the XXZ model with anisotropic
$J_{2}$ and $J_{d}$ terms [21, 22]. These results point to the dominant role
that a relatively large $J_{d}$ might play to stabilize the CSL though it
alone does not. Moreover, the scalar three-spin interaction
$\chi=\bm{S}_{i}\cdot\bm{S}_{j}\times\bm{S}_{k}$ (subscripts mean vertexes of
each triangle), which breaks TR symmetry, can naturally stabilize the CSL [20,
23]. Therefore, it suggests that the novel gapped QSL with TR symmetry we look
for does not exist in these models.

Another choice is to consider the Dzyaloshinskii-Moriya interaction (DM)
interaction, whose effects in a kagome lattice have recently been investigated
theoretically and experimentally [39, 40, 18, 41, 42, 43, 44, 45, 46].
However, it has been shown that this small interaction in fact favors the in-
plane $\bm{Q}=0$ magnetic order, when the direction of DM vector is
perpendicular to the $xy$ plane and leads to the local vector chirality
$\bm{\chi}_{\mathrm{v}}=\bm{S}_{1}\times\bm{S}_{2}+\bm{S}_{2}\times\bm{S}_{3}+\bm{S}_{3}\times\bm{S}_{1}$.
Hence, if proceeding along this direction, we need to consider additional
exchange interactions to melt this magnetic order. We notice that in this
magnetic order spins aline ferromagnetically along the diagonal direction of
the hexagon in the kagome lattice. So, it arises a possibility to consider the
interplay between the DM interaction and the third NN AFM interaction along
the diagonal direction.

Based on these considerations, in this work, we study the nature of the
quantum spin states in the $J_{1}$-$J_{2}$-$J_{3}$ kagome antiferromagnet with
additional DM interaction and the third NN AFM interaction $J_{d}$ along the
diagonal direction. We do not intend to obtain the global phase diagram with
so many spin interaction parameters, instead mainly focus on the study of
quantum spin states emerging out of the $\bm{Q}=0$ magnetic ordered state.
When only $J_{1}$ term exists, our numerical simulation suggests a Dirac spin
liquid (DSL) as the ground state, and the introduction of a weak DM
interaction transits the system into the long-range $\bm{Q}=0$ magnetic order.
These results are consistent with previous researches [13, 14, 15, 18, 41].
Therefore, without loss of generality, we will fix the magnitude of DM
interaction $D$ as $0\lesssim D\lesssim 0.2$ in our main calculations to
investigate the effects of other longer-range AFM Heisenberg interactions.
When a weak diagonal third NN interaction $J_{d}$ turns on, the magnitude $M$
of the magnetic moment for the $\bm{Q}=0$ state decreases and eventually drops
to zero when $J_{d}\simeq 0.21$ (setting $J_{1}=1$). Our calculation suggests
there is a continuous phase transition into a disordered phase, a gapped QSL
with TR symmetry. We also show that when a small additional $J_{2}$ is
present, this transition still occurs, for example a slightly larger
$J_{d}\simeq 0.27$ is required with $J_{2}=0.05$ because the $J_{2}$ term
favors the $\bm{Q}=0$ order, as mentioned above. Obviously, this state can
survive in a broad range of $J_{d}$. It is verified that the CSL and the so-
called $cuboc1$ magnetic state breaking TR symmetry are not found in the range
of $J_{d}$ we considered, while the $cuboc1$ order is indeed found in the
large $J_{d}$ range, such as $D=0.1,J_{2,3}=0,J_{d}=0.3$. Beyond the diagonal
third NN term $J_{d}$, we have also checked the effects of the usual NN term
$J_{3}$ and find that it further enhances the effect of suppressing the
$\bm{Q}=0$ order. Therefore, this phase transition is general in the sense of
the reasonable DM interaction and diagonal third NN term, and the resulting
gapped QSL with TR symmetry is robust.

To describe the intrinsic property of the gapped QSL, we calculate the TEE and
find $\gamma\simeq 0.748$ on a torus, which agrees numerically with the exact
value $\gamma_{\mathrm{ideal}}=\ln 2\simeq 0.693$. So, the quantum dimension
$D_{q}=2$ is obtained via $\gamma=\ln D_{q}$. In the meantime, we find that
its GSD is $n_{g}=4$. These results suggest that the gapped QSL holds an
Abelian topological order with $n_{g}=D_{q}^{2}$. It is expected that the
total Chern number of this QSL is zero because of TR symmetry. Interestingly,
we find that the spin-up and -down spinons see opposite gauge fluxes and there
is no coupling between them. As a result, we can independently define the
Chern numbers for the spin-up and -down spinons, respectively. It turns out
that the Chern number of the spin-up (-down) spinons is
$C_{\uparrow(\downarrow)}=+(-)1$, which is right the Z2 index [47]. So, the
spinons with different spins move along opposite directions on a given edge
with the opposite chiral central charges $c_{\pm}=\pm 1/2$. That shows that
this disordered state is an exotic gapped QSL with TR symmetry and shares the
same properties as those of a QSH state at the mean time. Thus, we name it the
spinon quantum spin Hall (SQSH) state. According to the topological long-range
entanglement, ground-state degeneracy and the Chern number of this SQSH state,
we suggest that it is a double-semion topological order (doubled Chern-Simons
state), which is described by a sum of two topological quantum field theories
with opposite chiralities and is suggested in the string-net model [48].
Generalizing to spin-$k$ ($k$ is integer) antiferromagnets, the spin liquids
with parity and time-reversal symmetry have been proposed, such as the doubled
CSL termed as CSL+CSL- with $k=1$ [49], which holds the same topological
properties with the SQSH state. In addition, this state has also been studied
by the wire deconstructionism of the long-range entanglement topological phase
[50].

##  II model and method

We start with the following model,

$\displaystyle H=$ $\displaystyle\
\sum_{\left(i,j\right)}J_{ij}\bm{S}_{i}\cdot\bm{S}_{j}+D\sum_{\langle
i,j\rangle}\bm{D}_{ij}\cdot\bm{S}_{i}\times\bm{S}_{j}.$ (1)

The first term in model (1) is the AFM Heisenberg spin interaction, including
the first, second and two anisotropic third NN spin couplings, expressed
respectively as the $J_{1}$, $J_{2}$, $J_{3}$ and $J_{d}$ terms. The second
term denotes the DM interaction connecting the first NN bond spins, and the
direction of the DM vector $\bm{D}_{ij}$ is perpendicular to the NN bond
$\langle i,j\rangle$ with $D$ its magnitude. All these terms are illustrated
in Fig. 1. In this paper, we only consider the case that the vector
$\bm{D}_{ij}$ is perpendicular to the lattice plane, namely only $D^{z}_{ij}$
is finite, so that total $S_{z}$ is conserved. We set $J_{1}=1$ as the energy
unit for convenience.

The model (1) will be investigated by the variational Monte Carlo (VMC)
method. First, we use the fermionic doublet representation to rewrite the spin
interactions as following,

$\displaystyle\bm{S}_{i}\cdot\bm{S}_{j}=-\frac{1}{4}(T_{ij}T_{ij}^{\dagger}+P_{ij}P_{ij}^{\dagger})+\mathrm{const},$
(2)
$\displaystyle\bm{S}_{i}\times\bm{S}_{j}=-\frac{i}{4}(T_{ij}\bm{T}_{ij}^{\dagger}+P_{ij}\bm{P}_{ij}^{\dagger}-\mathrm{h.c.})+\mathrm{const},$

where $T_{ij}=\psi^{\dagger}_{i}\psi_{j}$
($P_{ij}=\psi^{\dagger}_{i}\bar{\psi}_{j}$) is the singlet hopping (pairing)
term, while $\bm{T}_{ij}=\psi^{\dagger}_{i}\bm{\sigma}\psi_{j}$
($\bm{P}_{ij}=\psi^{\dagger}_{i}\bm{\sigma}\bar{\psi}_{j}$) the triplet
hopping (pairing) term. The fermionic doublet field and its particle-hole
partner are given by $\psi=(c_{\uparrow},c_{\downarrow}){{}^{\top}}$ and
$\bar{\psi}=(c^{\dagger}_{\downarrow},-c^{\dagger}_{\uparrow}){{}^{\top}}$,
respectvely. Considering the SU(2) gauge structure of this fermionic
representation [51], it is necessary to implement Lagrangian multipliers
$\bm{\lambda}$ to enforce the generators of the SU(2) gauge group
$\bm{\Lambda}_{i}=0$ to return the subspace of real physical states. Their
expressions with fermionic doublet representation are as following,

$\displaystyle\Lambda_{i}^{x}=-\frac{1}{4}(\psi_{i}^{\dagger}\bar{\psi}_{i}+\bar{\psi}_{i}^{\dagger}\psi_{i}),$
(3)
$\displaystyle\Lambda_{i}^{y}=-\frac{i}{4}(\psi_{i}^{\dagger}\bar{\psi}_{i}-\bar{\psi}_{i}^{\dagger}\psi_{i}),$
$\displaystyle\Lambda_{i}^{z}=\frac{1}{2}(1-\psi_{i}^{\dagger}\psi_{i}).$

Then, we use the fermionic parton approximation to decouple the spin
interactions into noninteracting quadratic structure, and obtain the mean-
filed Hamiltonian (irrelevant constants are omitted),

$\displaystyle H_{\mathrm{mf}}^{all}=$
$\displaystyle\sum_{i,j}(t_{ij}^{s}\psi_{i}^{\dagger}\psi_{j}+\bm{t}^{t}_{ij}\cdot\psi_{i}^{\dagger}\bm{\sigma}\psi_{j}+\Delta_{ij}^{s}\psi_{i}^{\dagger}\bar{\psi}_{j}$
(4)
$\displaystyle+\bm{\Delta}^{t}_{ij}\cdot\psi_{i}^{\dagger}\bm{\sigma}\bar{\psi}_{j}+\mathrm{H.c.})$
$\displaystyle+\sum_{i}\bm{\lambda}\cdot\bm{\Lambda}_{i}-\bm{M}_{i}\cdot\psi^{\dagger}_{i}\bm{\sigma}\psi_{i}/2,$

where $t_{ij}^{s}$, $\bm{t}^{t}_{ij}\parallel\bm{D}_{ij}$ are spinon hopping
parameters and $\Delta_{ij}^{s}$, $\bm{\Delta}^{t}_{ij}\parallel\bm{D}_{ij}$
are pairing ones, and a background field $\bm{M}_{i}$ is applied to induce a
static magnetic order. Therefore, all the variational parameters are
$\alpha^{all}=(t_{ij}^{s},\bm{t}^{t}_{ij},\Delta_{ij}^{s},\bm{\Delta}^{t}_{ij},\bm{\lambda},\bm{M}_{i})$.

Obviously, there must be various different ansatzes from the mean-field
Hamiltonian Eq. (4) with a plenty of variational parameters. We selectively
consider various singlet and triplet hopping terms and several pairing terms
combined with the projective symmetry group[52, 53, 54] (PSG).

Actually, we have found that the spinon-pairing instability (the allowed first
NN triplet pairing term with complex numbers) is vanishingly small
($|\bm{\Delta}_{\langle ij\rangle}^{t}|/t_{\langle ij\rangle}^{s}<10^{-3}$) in
our calculation. So, we will ignore all pairing terms reasonably. Moreover,
the $\bm{\lambda}$ term can also be ignored because of vanishing pairing terms
in the actual VMC procedure. In this way, we can obtain the reduced mean-field
Hamiltonian, which is written as,

$\displaystyle H_{\mathrm{mf}}=$
$\displaystyle\sum_{i,j}(t_{ij}^{s}\psi_{i}^{\dagger}\psi_{j}+\mathrm{H.c.})+\sum_{\langle
ij\rangle}(\bm{t}^{t}_{ij}\cdot\psi_{i}^{\dagger}\bm{\sigma}\psi_{j}+\mathrm{H.c.})$
(5)
$\displaystyle-\sum_{i}\bm{M}_{i}\cdot\psi^{\dagger}_{i}\bm{\sigma}\psi_{i}/2,$

where $t_{ij}^{s}$, $\bm{t}^{t}_{ij}\parallel\bm{D}_{ij}$ are the singlet and
triplet spinon hopping parameters, respectively. The former includes the first
and second NN hopping terms and the later only include the first NN one.

With the mean-field ground-state wave function
$|\mathrm{GS(\alpha)}\rangle_{\mathrm{mf}}$, we can obtain a trial wave
function
$|\Phi(\alpha)\rangle=P_{G}|\mathrm{GS(\alpha)}\rangle_{\mathrm{mf}}$, where
$P_{G}$ is the Gutzwiller projection to guarantee the single occupancy
condition, and $\alpha=(t_{ij}^{s},\bm{t}^{t}_{ij},\bm{M}_{i}$) are
variational parameters. We emphasize that these parameters could be complex
numbers in principle, so the whole parameter range is in fact large. In the
variational process, we employ the stochastic reconfiguration scheme [55] to
optimize so many parameters. We have considered various Z2 QSLs, U(1) QSLs,
$\bm{Q}=0$ and $cuboc1$ magnetic orders as initial trial states selectively
(see Appendix A for details). We adopt the torus geometry with
$L_{1}=L_{2}=12$ for main results, where $L_{1,2}$ are the lengths along the
two Bravais kagome-lattice vectors $\bm{a}_{1}=(1,0)$ and
$\bm{a}_{2}=(-1/2,\sqrt{3}/2)$.

##  III result

In this paper, we do not intend to obtain the comprehensive phase diagram of
the spin model in the kagome lattice defined by model (1). Instead, we will
focus on the study of the possible states stabilized or induced by the
additional DM interaction $D$ and the AFM interactions $J_{d}$ across the
diagonals of the hexagons of the kagome lattice, in the presence of the first,
second and usual third NN spin couplings.

(a)

(b)

(c)

(d)

Figure 1: (a) Illustration of the spin interactions included in model (1) on
kagome lattice. The dashed lines with different colors indicate different
Heisenberg terms. The arrows on the nearest-neighbor bonds indicate the
direction of the DM interaction and the DM vector $\bm{D}_{ij}$ is oriented
parallel (antiparallel) to $z$ axis in up (down) triangles indicated by
$\bigodot$ ($\bigotimes$). (b) denotes the ansatz for the exotic SQSH state.
The $+(-)\theta$ is the flux for spin-up (-down) spinon in each triangle,
respectively, and the $\pi-(+)2\theta$ is the same way in each hexagon. The
marked red and blue dashed bonds denote that the hopping terms along these
bonds in $H_{\mathrm{mf}}$ have the opposite sign compared with those unmarked
ones so that the unit cell is doubled. (c) is the schematic diagram indicating
the in-plane configuration of classical spins in the long-range magnetic order
with $\bm{Q}=0$, three spins in all triangles form a $120^{\circ}$
distribution, which preserves the original translational symmetry. (d)
illustrates the energy curves versus $D$ when $J_{2,d,3}=0$ (see Appendix A
for the details of the four states). We note that all of the standard errors
in this paper are considered as confidence intervals.

We start with the result in the presence of only the first NN $J_{1}$ AF
couplings. Our VMC results suggest that a Dirac spin liquid is the most
energetically favored state in this case, which is consistent with previous
results[13, 14, 15]. In this U(1) QSL, all triplet terms vanish and only the
singlet hopping terms survive. In the case of only the first NN hopping term
is included, there are zero ($\pi$) fluxes through triangles (hexagons) in the
kagome lattice, corresponding to the pattern for $\theta=0$ in Fig. 1. We also
consider the Jastrow-type wave functions for the magnetic ordered state, which
is widely used in VMC studies. This state is constructed based on the solution
with only a finite $\bm{M}_{i}$ term in $H_{\mathrm{mf}}$ (see Appendix A). It
has the magnetic order $\bm{Q}=0$ in the lattice (XY) plane as illustrated in
Fig. 1 , and is called as Jastrow + $(\bm{Q}=0)$ state. In the case of only
$J_{1}$ term, the energy of this magnetic order is obviously higher than that
of the Dirac spin liquid, as shown in Fig. 1.

When turning on the DM interaction $D$, we find that the energy of the Dirac
spin liquid does not depend on $D$, while the energy of the Jastrow +
$(\bm{Q}=0)$ state decreases nearly linearly with $D$ (see Fig. 1). However,
before the Jastrow + $(\bm{Q}=0)$ state surpasses energetically the Dirac spin
liquid, a novel competing state emerges. This state has the same forms of the
first and second NN singlet hopping terms as the Dirac spin liquid, but has a
finite and pure imaginary first NN triplet hopping term. The selection of this
state is done via a comprehensive comparison with other possible states. We
have considered the ansatzes with complex first and second NN singlet hopping
terms, but find that their imaginary parts are almost zero
($\mathrm{Im}(t_{\langle ij\rangle}^{s})/t_{\langle ij\rangle}^{s}<10^{-2}$,
$\mathrm{Im}(t_{\langle\langle ij\rangle\rangle}^{s})/t_{\langle
ij\rangle}^{s}<10^{-2}$) in our variational process. We have also checked the
existence of the first NN triplet pairing term (complex number), and it turns
out that this term vanishes ($|\bm{\Delta}_{\langle ij\rangle}^{t}|/t_{\langle
ij\rangle}^{s}<10^{-3}$) in our calculation. In particular, we have considered
two candidate states, which are regarded as the derivative states of the
uniform resonating valence bond state (see Appendix A). One candidate carries
the same fluxes in all triangles and the other carries the opposite fluxes in
the up and down triangles. Combining the triplet hopping terms (complex
number) and finite $\bm{M}_{i}$, we find that both these states are
energetically unfavored.

In this novel state, the complex triplet hopping term will lead to the
consequence that free spinons at the VMC mean-field level will carry nonzero
flux when hops along closed loops, as the direction of the DM vector
$\bm{D}_{ij}$ considered in this paper is perpendicular to the lattice plane.
Specifically, the spinon with up (down) spin carries a $+(-)\theta$ flux in
all triangles and $\pi-(+)2\theta$ flux in all hexagons. This shares the same
physics as a quantum spin Hall state or topological insulator. So, we dub it
as the spinon quantum spin Hall (SQSH) state and will discuss its properties
in detail later. When $J_{2}=J_{d}=J_{3}=0$, this state emerges firstly at
$D\sim 0.01$ in the sense that our numerical calculations can determine it.
However, the mix of this state with the $\bm{Q}=0$ magnetic order state, i.e.
SQSH + ($\bm{Q}=0$) state, is always lower in energy in the range
$0.01\lesssim D\lesssim 0.2$. Moreover, this magnetic order state is
energetically more favored than the Jastrow + ($\bm{Q}=0$) state, as shown in
Fig. 1.

The SQSH + ($\bm{Q}=0$) state is induced and stabilized by the DM interaction.
It is also analogous to the most possible ground state of the AFM $J_{1}$
Heisenberg model on the triangular lattice [56, 57, 58]. As the DM vector
$\bm{D}_{ij}$ is considered to be along the perpendicular direction of the
lattice plane, it confines the spins to lie in the lattice plane and
effectively plays a role of the easy-plane anisotropy. In the meantime, it is
able to induce the local vector chirality $\bm{\chi}_{\mathrm{v}}$. These two
properties make the DM interaction to favor the in-plane $\bm{Q}=0$ magnetic
order. In detail, this state requires a finite $\bm{M}_{i}$ and additionally
the same hopping term as the SQSH state. However, due to the existence of the
magnetic order, the elementary excitation of this SQSH + ($\bm{Q}=0$) state is
the magnon which is the confined state of two spinons. Therefore, there are no
free spinons existing in the sense of the VMC mean-field level and no SQSH
effects.

(a)

(b)

(c)

Figure 2: Magnitude of magnetic moment ($M$) in the SQSH + ($\bm{Q}=0$) state
versus $J_{d}$ term for different $J_{2},J_{3}$ interactions. (a) and (b) is
obtained with the DM interaction $D=0.1$, (c) with $D=0.2$. We note that the
legends of all these curves are labeled by $(J_{2},J_{3})$ and the lattice
size we employ in (b) and (c) is $12\times 6\times 3$.

It turns out that we need to melt the $\bm{Q}=0$ magnetic order in order to
realize SQSH state. We find that the AFM interaction $J_{d}$ across the
diagonals of hexagons can effectively weaken the magnitude of the magnetic
moment $M$. As can be seen in Fig. 2, $M$ decreases nearly linearly with
$J_{d}$ and vanishes at $J_{d}\simeq 0.21$, when only the DM interaction
besides the first NN term is considered, namely $D=0.1$ and $J_{2}=J_{3}=0$
[Fig. 2)]. The introduction of the $J_{2}$ term enhances effectively the
moment of the $\bm{Q}=0$ magnetic order, so that a larger $J_{d}$ is needed to
eliminate the magnetic order. As shown in Fig. 2, when $J_{2}$ starts from
zero to 0.05, the critical $J^{c}_{d}$ term leading to $M\simeq 0$ increases
to be 0.27. On the other hand, the effect of the $J_{3}$ term behaves
conversely with $J_{2}$, but cooperatively with $J_{d}$ [Fig. 2]. For example,
if the $J_{2}$ term is fixed to be $0.1$, we can find that the moment $M$ with
$J_{3}=0.1$ is significantly smaller than that with $J_{3}=0.05$ for the same
$J_{d}$, and $J^{c}_{d}$ drops to be near 0.16 for a fixed $J_{3}=J_{2}=0.1$.
To show more clearly the effect of the $J_{3}$ term, we present the results
for the moment $M$ in the case of $D=0.1$, $J_{2,d}=0$ in Table 1. To check
the effects of the DM interaction, we show results obtained with $D=0.2$ in
Fig. 2, one can see that $J^{c}_{d}$ increases noticeably compared to the case
of $D=0.1$ for various $J_{2}$ and $J_{3}$. We note that those spin exchange
parameters to melt the $\bm{Q}=0$ magnetic order can be nearly one magnitude
smaller than the dominant nearest-neighbor exchange interaction $J_{1}$, which
are considered to be acceptable physically and realizable experimentally.

$J_{3}$ |  0.05 |  0.08 |  0.1 |  0.12 |  0.15
---|---|---|---|---|---
$M$ |  0.2074 |  0.1239 |  0.0728 |  0.0279 |  0.0062
Table 1: Magnetic moment $M$ in the SQSH + ($\bm{Q}=0$) magnetic order state calculated with lattice size $12\times 12\times 3$ in the case of $D=0.1$, $J_{2,d}=0$. state |  $D$ |  $J_{2,3}$ |  $J_{d}$ |  $E$
---|---|---|---|---
SQSH | 0.1 | 0.0 | 0.29 | -0.43371
0.1 | 0.0 | 0.3 | -0.43345
$cuboc1$ | 0.1 | 0.0 | 0.29 | -0.43371
0.1 | 0.0 | 0.3 | -0.43379
Table 2: Variational energy (per site) for the SQSH state and $cuboc1$ magnetically ordered state with lattice size $12\times 12\times 3$. When $D=0.1$, $J_{2,3}=0.0$ and $J_{d}\gtrsim 0.3$, the $cuboc1$ order state with TR symmetry breaking will be dominant. All the errors of the energies are $\sim 10^{-5}$. state |  $D$ |  $J_{2}$ |  $J_{d}$ |  $J_{3}$ |  $E$
---|---|---|---|---|---
SQSH | 0.1 | 0.05 | 0.3 | 0 | -0.43444
0.1 | 0.05 | 0.32 | 0 | -0.43382
$cuboc1$ | 0.1 | 0.05 | 0.3 | 0 | -0.43443
0.1 | 0.05 | 0.32 | 0 | -0.43387
Table 3: Variational energy (per site) for the SQSH state and $cuboc1$
magnetically ordered state with lattice size $12\times 12\times 3$. When
$D=0.1$, $J_{2}=0.05$, $J_{d}\gtrsim 0.32$ and $J_{3}=0.0$, the $cuboc1$
ordered state with TR symmetry breaking will be dominant. All the errors of
the energies are $\sim 10^{-5}$.

After the $\bm{Q}=0$ magnetic order is melted, we further study the robustness
of this disordered SQSH state against the $cuboc1$ order and present the
results in Table 2 and 3 . It shows that the SQSH state can survive in a
relatively broad range of $J_{d}$ and eventually gives way to the $cuboc1$
order, such as for $J_{d}\gtrsim 0.3$. In addition, we find that the
transition from the finite $M$ magnetic order state to the $M=0$ spin liquid
state with SQSH is a continuous phase transition. Now, let us discuss the
detail properties of the pure SQSH state. As stated above, this state is a
quantum spin liquid and inherits the same singlet hopping pattern as the Dirac
quantum spin liquid [See Fig. 1]. The key difference between them is that the
SQSH state includes a pure imaginary $z$ component in the first NN triplet
hoppings induced by the DM interaction. As a result, the spin-up and -down
spinons see the opposite flux ($\pm\theta$) in all triangles. Because of the
absence of the coupling between spin-up and -down spinons, we can rewrite the
mean-field Hamiltonian,

$H_{\mathrm{mf}}=\begin{pmatrix}h&0\\\ 0&h^{*}\end{pmatrix},$ (6)

where $h$ denotes the Hamiltonian for the spin-up spinons while $h^{*}$ that
for the spin-down spinons. When we go into the $\bm{k}$ space, at an arbitrary
$\bm{k}$, there is always a couple of degenerate states, which are conjugate
to each other. So, the spin-up and -down spinons have opposite Berry phase in
any plaqutte of $\bm{k}$ space, and consequently the opposite Chern number. We
have calculated the Chern number of the filled three bands for the spin-up
(-down) spinons as shown in Fig. 3, and find that
$C_{\uparrow(\downarrow)}=+(-)1$. We then calculate the energy dispersions of
the spin-up (-down) spinons in the SQSH state with the period-open boundary
condition, and the results are shown in Fig. 4.

(a)

(b)

Figure 3: (a) Mean-field band structure of the SQSH state, where each band is
doubly degenerate so that there is a significant energy gap because of the one
spinon per site. (b) The chirality-chirality correlation of the SQSH state
$|\langle\chi_{i}\chi_{j}\rangle|$ ($i\neq j$) with respect to the distance
along $\bm{a}_{1}$ for different lattice sizes. And the
$|\bm{r}_{i}-\bm{r}_{j}|$ means the distance between two up-pointing
triangles. We only need to consider the distance up to 6, 9 and 12 because of
the period boundary condition for a system of $L_{1,2}=12$, $L_{1}=18,L_{2}=8$
and $L_{1}=24,L_{2}=6$, respectively. These results are obtained in the system
with $D=0.1$, $J_{d}=0.21$ and $J_{2}=J_{3}=0$.

It shows clearly that two edge states emerge in the gap of the energy bands
for each spin species spinon. In particular, the energy bands of the spin-up
and -down spinons are antisymmetric with respect to momentum $k$, indicating
that the spin-up state is just the time-reversal copy of the spin-down one.
So, the spin-up and -down spinons move along opposite direction on the edge.
We have also checked the chirality-chirality correlation in the SQSH state,
the results are shown in Fig. 3. It shows
$|\langle\chi_{i}\chi_{j}\rangle|\sim 0$ within numerical error, so the SQSH
indeed has TR symmetry. It suggests that the two edge states of this
topological SQSH are protected by TR symmetry. Thus, we believe that the SQSH
state is the spinon version of the quantum spin Hall state.

Figure 4: Left (right) panel is the spin-up (-down) energy dispersion of the
SQSH state with the same optimal parameters as that of Fig. 3, respectively.
In the calculation, the period (open) boundary condition is used in the
direction of basis vectors $\bm{n}_{1}=(2,0)$ ($\bm{n}_{2}=(-1/2,\sqrt{3}/2)$)
defined in the doubled unit cell and the reciprocal vector
$\mathbf{b}_{1}=(\pi,\pi/\sqrt{3})$.

On the other hand, the SQSH state is essentially the strongly correlated
state, which is different from the conventional quantum spin Hall state based
on the single-particle picture. Its elementary excitations are fractionalized
spinons out of the QSL ground state and could embrace more intrinsic
topological properties. To further explore its property, we calculate the
topological entanglement entropy and ground-state degeneracy. The TEE is
obtained numerically by partitioning the system into two subsystems and
calculating the second order Renyi entropy (see Appendix C for details). Then,
the Renyi entropy is expected to follow
$S(\mathcal{L})=\alpha\mathcal{L}-\gamma$, where $\alpha$ depends on the
details of the state, $\mathcal{L}$ represents the boundary length of a
contractible patch with codimension-1 boundary in the system and $\gamma$ is
the universal TEE. In order to eliminate the area-law contribution
$\alpha\mathcal{L}$, we calculate the entanglement entropy of plaquette
$\mathrm{P}_{1}$ (the shaded region in the inset of Fig. 5) with different
sizes, and the result is presented in Fig. 5. It shows that $S(\mathcal{L})$
increases linearly with $\mathcal{L}$. Then, we apply a linear extrapolation
to $\mathcal{L}\rightarrow 0$, and obtain the TEE for the SQSH state as
$\gamma\approx 0.748$, which is very close to $\ln 2=0.693$. It shows that the
SQSH state has intrinsic topological order. This is different from the quantum
spin Hall state which is a symmetry-protected topological state.

Figure 5: Entanglement entropy as a function of the area for the SQSH state with lattice size $12\times 12\times 3$. The corresponding optimal parameters are obtained when $D=0.1$, $J_{2}=0.05$, $J_{3}=0$ and $J_{d}=0.27$. The $x$-axis labeled as $\mathcal{L}$ means that the area is $\mathcal{L}^{2}$ times of the primitive cell, i.e., the shaded rhombus as shown in the insert, where different colored sites label the different sublattices of kagome lattice. state |  $\varepsilon_{1}$ |  $\varepsilon_{2}$ |  $\varepsilon_{3}$ |  $\varepsilon_{4}$ |  $n_{g}$
---|---|---|---|---|---
$\bm{Q}=0$ | 3.997 | 1.3$\times 10^{-3}$ | 9.72$\times 10^{-4}$ | 9.27$\times 10^{-4}$ | 1
SQSH | 1.985 | 0.6774 | 0.673 | 0.6646 | 4
Table 4: Eigenvalues and GSDs of the $\bm{Q}=0$ and SQSH states, calculated in
a kagome lattice of 12$\times$12$\times$3\. And $\varepsilon_{i}$ denotes the
eigenvalue of the overlap matrix while $n_{g}$ denotes GSD.

The nontrivial topological properties of this state can be further
characterized by its GSD. The calculation of GSD is carried out by
constructing four states $|\phi_{\pm,\pm}\rangle_{\mathrm{mf}}$ with a
Gutzwiller projection, where the + (-) subscript denotes the period (anti-
period) boundary condition along the directions of two lattice basis vectors,
respectively. After diagonalizing the overlap matrix on the basis vector of
these four states, we can obtain their eigenvalues, and the number of the
significant finite eigenvalues is just GSD. In Table. 4, the results of four
eigenvalues and the corresponding GSD are presented for the $\bm{Q}=0$
magnetic order state and SQSH state. It shows that there is only one
significant finite eigenvalue for the $\bm{Q}=0$ state, so its GSD is
$n_{g}=1$. In the category of topological order [59, 60], this magnetic order
only holds a trivial (identity) topological excitation $\mathbb{I}$ with its
quantum dimension $D_{q}=1$ due to $n_{g}=1$. According to the relation
$\gamma=\ln D_{q}$, we find its topological entanglement entropy is
$\gamma=0$. Besides, its Chern number is found to be zero. So, the $\bm{Q}=0$
magnetic order state is a topological trivial phase without long-range
entanglement with $\gamma=0$. For the SQSH state, there are four finite
eigenvalues and its GSD is $n_{g}=4$. From the above result of the TEE for the
SQSH state $\gamma\approx\ln 2$, we obtain its quantum dimension $D_{q}=2$.
For an Abelian topological phase, GSD = $D_{q}^{2}$. Therefore, we suggest
that this SQSH state is an Abelian topological phase.

##  IV conclusions

In summary, we have studied the interplay between the Dzyaloshinskii-Moriya
interaction and the long-range AFM Heisenberg interactions in the kagome
antiferromagnet by using the variational Monte Carlo method. We find that the
Dzyaloshinskii-Moriya interaction alone favors the $\bm{Q}=0$ long-range
magnetic order, and an additional antiferromagnetic interaction across the
diagonals of the hexagons of the kagome lattice can suppress and eliminate
eventually this order phase. This topological phase transition leads to an
exotic quantum spin state with fruitful topological properties. We elaborate
that it is a topological gapped quantum spin liquid with the ground-state
degeneracy $n_{g}=4$ and the topological entanglement entropy $\gamma=\ln 2$.
As the fractionalized excitations in a quantum spin liquid, the spinons
constitute the two chiral edge states protected by the time-reversal symmetry.
So, the spin-up and -down spinons move along opposite directions on a given
edge, and it gives rise to the quantum spin Hall effect existing in a
topological insulator.

We suggest that, by doping magnetic impurities into this spinon quantum spin
Hall state to suppress one of the so-call helical states and retain the
topological properties of the system as has been done for the celebrated
quantum anomalous Hall effect [61], it is possible to detect the
fractionalized spinons in this exotic quantum spin liquid.

###### Acknowledgements.

We would like to thank Q.-H. Wang, Z.-X. Liu, Z.-L. Gu, X.-M. Cui and J.-B.
Liao for many helpful and valuable discussions. This work was supported by
National Key Projects for Research and Development of China (Grant No.
2021YFA1400400) and the National Natural Science Foundation of China (No.
92165205).

##  Appendix A MEAN-FIELD ANSATZES

1\. $\bm{Q}=0$ magnetic order

Firstly, we consider a simple classical magnetic order $\bm{Q}=0$ in the
lattice (XY) plane. This state is constructed based on the solution with only
a finite $\bm{M}_{i}$ term in $H_{\mathrm{mf}}$, and is restricted to the
$S_{tot}^{z}=0$ subspace [14] by the application of the projector
$P_{S_{tot}^{z}=0}$. Then, quantum fluctuations are included by the long-range
Jastrow projector:

$P_{J}^{z}=\mathrm{exp}\left(1/2\sum_{ij}\mu_{ij}S_{i}^{z}S_{j}^{z}\right),$
(7)

where $\mu_{ij}$ is the pseudopotential that only depends on the absolute
distance $|\bm{R}_{i}-\bm{R}_{j}|$ of two sites. It decays exponentially with
the distance so that we can just consider the first three without loss of
generality. Finally, we get the state as
$|\mathrm{\Psi_{\mathrm{mag}}}\rangle=P_{S_{tot}^{z}=0}P_{J}^{z}|\mathrm{GS}\rangle_{\mathrm{mf}}$,
and call it the Jastrow + $(\bm{Q}=0)$ state.

2\. DSL

The Dirac quantum spin liquid (DSL) is one of the most competitive candidate
ground states in the Heisenberg model with only $J_{1}$ term on kagome
lattice. In this U(1) QSL, all triplet terms vanish and only singlet hopping
terms survive. If we only consider the first NN hopping ones, there are zero
($\pi$) fluxes through triangles (hexagons), respectively. Based on PSG, the
second singlet hopping terms can also exist. All bond patterns are shown in
Fig. 1 in the main text. We also emphasize that the third NN terms in this
state are forbidden by PSG, so that we rationally abandon them. In fact, the
second NN singlet pairing terms are allowed by PSG, namely, the so-called
Z${}_{2}[0,\pi]\beta$, a gapped QSL [53]. But, it’s not energy favorable in
our model by our calculation. In detail, our numerical results suggest that
the $J_{1}$, $J_{d}$ and $J_{3}$ interactions suppress this pairing term,
which is consist with Ref. [14]. For this reason, we also throw away the
second NN singlet pairing terms all the time.

3\. SQSH

The SQSH state is a novel state we found in this paper. In this state, the
first and second NN singlet hopping terms have the same forms as those of the
DSL discussed above, but the first NN triplet hopping term is finite and pure
imaginary number. According to the direction of the DM interaction vector
$\bm{D}_{ij}$ we chosen in this paper, this triplet hopping term will lead to
that the spin-up (-down) spinon sees $+(-)\theta$ flux in all triangles and
$\pi-(+)2\theta$ flux in all hexagons. Intrinsically, the spin-up and -down
spinons have opposite Chern number $C_{\uparrow,\downarrow}=\pm 1$.

As shown in Fig. 1, when $0.01\lesssim D<0.2$ and other interactions beyond
$J_{1}$ term are absent, another magnetic order state as a magnetic
instability of SQSH, we call it the SQSH + ($\bm{Q}=0$) state, is
energetically favored. For convenience, we stipulate the unique $\bm{Q}=0$
magnetic ordered state mentioned in the main text is just the SQSH +
($\bm{Q}=0$) state, because the Jastrow + ($\bm{Q}=0$) state is not energy
favorable. Henceforth, we also call it the $\bm{Q}=0$ state.

4\. CSL

We have considered the chiral quantum spin liquid (CSL), in which the first
and second NN singlet hopping terms are complex numbers as discussed in Ref.
[20]. But, we find that their imaginary parts are almost zero
($\mathrm{Im}(t_{\langle ij\rangle}^{s})/t_{\langle ij\rangle}^{s}<10^{-2}$,
$\mathrm{Im}(t_{\langle\langle ij\rangle\rangle}^{s})/t_{\langle
ij\rangle}^{s}<10^{-2}$) in our variational process. At least, in the range of
interactions we considered in this work, this result excludes the chiral QSL
stabilized by strong enough $J_{d}$ [20, 23].

5\. $cuboc1$ magnetic order

We also study the so-called $cuboc1$ magnetic order as the instability of the
TR symmetry breaking chiral spin liquid, i.e., CSL + $cuboc1$ [17]. Without
loss of generality, we additionally include the first NN triple hopping term
$t_{<ij>,z}^{t}$ into this ansatz. The detail of the classical $cuboc1$ order
is shown in Fig. 6. For the sake of simplicity, we call it $cuboc1$ state. We
indeed find that the variational energy of this magnetic order is lower than
that of the SQSH state, when the $J_{d}$ term is relatively large, as listed
in Table 2 and 3. Before the phase transition, we find that the $cuboc1$ order
almost reduces to SQSH state by VMC calculation because the optimal parameters
$t_{<ij>}^{s}$ and $t_{<<ij>>}^{s}$ are real number, the optimal parameter
$t_{<ij>,z}^{t}$ is pure imaginary and the magnetic moment $M$ is almost
vanishing. As a result, the energies of both $cuboc1$ order and SQSH state
look like degeneracy within numerical error. Therefore, we believe this exotic
SQSH state can survive in a relatively broad range of $J_{d}$.

Figure 6: 12 sublattices of the magnetic unit cell for classical $cuboc1$
order. $\bm{S}_{1}=(1,0,0)$, $\bm{S}_{2}=(-\frac{1}{2},\frac{\sqrt{3}}{2},0)$,
$\bm{S}_{3}=(-\frac{1}{2},\frac{\sqrt{3}}{6},\frac{\sqrt{6}}{3})$,
$\bm{S}_{4}=(0,\frac{\sqrt{3}}{3},-\frac{\sqrt{6}}{3})$,
$\bm{S}_{5}=(\frac{1}{2},-\frac{\sqrt{3}}{2},0)$,
$\bm{S}_{6}=(-\frac{1}{2},-\frac{\sqrt{3}}{2},0)$,
$\bm{S}_{7}=(1,-\frac{\sqrt{3}}{3},\frac{\sqrt{6}}{3})$,
$\bm{S}_{8}=(-\frac{1}{2},-\frac{\sqrt{3}}{6},-\frac{\sqrt{6}}{3})$,
$\bm{S}_{9}=(\frac{1}{2},-\frac{\sqrt{3}}{6},-\frac{\sqrt{6}}{3})$,
$\bm{S}_{10}=(-1,0,0)$,
$\bm{S}_{11}=(\frac{1}{2},\frac{\sqrt{3}}{6},\frac{\sqrt{6}}{3})$ and
$\bm{S}_{12}=(\frac{1}{2},\frac{\sqrt{3}}{2},0)$.

6\. uRVB state

In addition, we have considered another two candidate states. Both of them are
regarded as the derivative states of uRVB states. Without loss in generality,
we allow the first NN singlet $t_{ij}^{s}$ terms are complex number so that
there are gauge fluxes though triangular and hexagonal plaquette. And we
consider two flux patterns, one where fluxes through all triangles are the
same and another one where fluxes through up and down triangles are the
opposite. Combining the triplet hopping $t_{ij,z}^{t}$ terms (complex number)
and finite $\bm{M}_{i}$ of the $\bm{Q}=0$ magnetic order, we find that the
trial energies of both states are much higher than aforementioned states. So,
we discard them and do not draw up their energy curves in Fig. 1 in main text.

##  Appendix B  CHERN NUMBER

Nonzero Chern number is one of the fundamental topological numbers to
characterize topological phase of matter. Here we won’t go into details about
the concepts of Berry connection, Berry phase and Chern number with formal
analytical expression. We just introduce the numerical calculation of Chern
number in the filled bands [62].

Firstly, we obtain the corresponding mean-field Hamiltonian with the optimized
variational parameters by VMC and then transform it into $\bm{k}$-space,
$H(\bm{k})$. And we emphasize that the $H(\bm{k})$ is periodic along the
directions of reciprocal basis vectors $\bm{b}_{1,2}$,
$H(\bm{k})=H(\bm{k}+n_{1}\bm{b}_{1}+n_{2}\bm{b}_{2})$, where $n_{1,2}$ are any
integer numbers. In other word, $H(\bm{k})$ is in Bloch form. Therefore, this
Fourier transformation must be handled with care. To be specific, it usually
needs a gauge transformation, $c_{\bm{k}}\longrightarrow
c_{\bm{k}}e^{i\bm{k}\cdot\bm{\delta}}$. In general, the $\bm{\delta}$ is
different for different $c_{\bm{k}}$ and not unique.

Then, for a lattice with finite size, the Brillouin zone is filled with
discrete $\bm{k}$ points. And we define intervals of $\bm{k}$ points in two
directions of reciprocal basis vectors,

$\bm{u}_{i}=\frac{l_{i}\bm{b}_{i}}{2N_{i}\pi},\quad(i=1,2;\quad
N_{i}/l_{i}\in\bm{\mathrm{N}}^{*}).$ (8)

In our calculation, we take $l_{i}=1$ to guarantee the highest numerical
precision. We also note larger intervals are also allowed as long as the
result is convergent. And then, we require that the eigenstate
$|n(\bm{k})\rangle$ of $H(\bm{k})$ is also periodic in the Brillouin zone to
eliminate the effect of any U(1) gauge of eigenstate. Now we can define the
U(1) quantity for a certain $\bm{k}$ as following,

$\eta(\bm{k})_{\bm{u}_{i}}\equiv\frac{\langle
n(\bm{k})|n(\bm{k}+\bm{u}_{i})\rangle}{|\langle
n(\bm{k})|n(\bm{k}+\bm{u}_{i})\rangle|}.$ (9)

$\eta(\bm{k})_{\bm{u}_{i}}$ is well defined as long as the denominator of Eq.
9 is nonzero. And then, we can define another variable about phase in a loop
with $\eta(\bm{k})_{\bm{u}_{i}}$,

$\displaystyle\theta(\bm{k})=\frac{1}{i}\ln\left(\eta(\bm{k})_{\bm{u}_{1}}\eta(\bm{k}+\bm{u}_{1})_{\bm{u}_{2}}{\eta(\bm{k}+\bm{u}_{2})^{\dagger}_{\bm{u}_{1}}}{\eta(\bm{k})^{\dagger}_{\bm{u}_{2}}}\right),$
(10) $\displaystyle-\pi<\theta(\bm{k})\leq\pi.$

Finally, the Chern number of the $n$th filled band is obtained by,

$C_{n}\equiv\frac{1}{2\pi}\sum_{\bm{k}\in\mathrm{BZ}}\theta(\bm{k}).$ (11)

##  Appendix C  GROUND-STATE DEGENERACY AND TOPOLOGICAL ENTANGLEMENT ENTROPY

The low-energy gauge fluctuations of gapped QSLs (topological orders) are
characterized by ground-state degeneracy (GSD). When one compacts the lattice
to a torus, in the thermodynamics limit, there is no energy cost when a Z2
$\pi$ flux is inserted in any hole of the torus. In the mean-field theory,
this procedure is equivalent to changing the boundary condition of
$H_{\mathrm{mf}}$ from period to anti-period. For a two-dimensional system, in
general, we can always construct four states
$|\phi_{\pm,\pm}\rangle_{\mathrm{mf}}$, where the + (-) subscript denotes the
periodic (anti-periodic) boundary condition along the directions of two
lattice basis vectors, respectively. Then, we enforce a Gutzwiller projection
to these four ground states of $H_{\mathrm{mf}}$ to recover physical Hilbert
space. Thus, we rewrite the symbols of the four states for convenience,

$\displaystyle|1\rangle=P_{G}|\phi_{+,+}\rangle_{\mathrm{mf}},\quad|2\rangle=P_{G}|\phi_{+,-}\rangle_{\mathrm{mf}},$
(12)
$\displaystyle|3\rangle=P_{G}|\phi_{-,+}\rangle_{\mathrm{mf}},\quad|4\rangle=P_{G}|\phi_{-,-}\rangle_{\mathrm{mf}}.$

We can calculate the 4 by 4 overlap matrix $\mathcal{O}$ based on these four
states. In detail, the matrix element $\mathcal{O}_{ij}=\langle
i|j\rangle/\sqrt{\langle i|i\rangle\langle j|j\rangle}$, where $i,j=1,2,3,4$.
After diagonalizing this overlap matrix, we can obtain its eigenvalues. The
number of the significantly finite eigenvalues is just GSD. We calculate the
overlap matrices of the $\bm{Q}=0$ magnetically ordered state and SQSH states
as following,

$\displaystyle\mathcal{O}_{\bm{Q}=0}\simeq\begin{pmatrix}1&e^{-i0.29}&e^{-i0.79}&e^{-i0.76}\\\
e^{i0.29}&1&e^{-i0.5}&e^{-i0.46}\\\ e^{i0.79}&e^{i0.5}&1&e^{i0.03}\\\
e^{i0.76}&e^{i0.46}&e^{-i0.03}&1\end{pmatrix},$ (13)
$\displaystyle\mathcal{O}_{\mathrm{SQSH}}\simeq$
$\displaystyle\begin{pmatrix}1&0.32e^{i2.64}&0.33e^{-i2.13}&0.33e^{i3.05}\\\
0.32e^{-i2.64}&1&0.32e^{i1.51}&0.33e^{i0.42}\\\
0.33e^{i2.13}&0.32e^{-i1.51}&1&0.33e^{-i1.1}\\\
0.33e^{-i3.05}&0.33e^{-i0.42}&0.33e^{i1.1}&1\end{pmatrix}.$

And now, we can obtain their eigenvalues and GSDs, as shown in Table 4.

Another important quantity to characterize topological order is topological
entanglement entropy (TEE) [5, 6]. Firstly, we divide a system into two parts,
as shown in Fig. 7.

(a)

(b)

Figure 7: (a) shows a schematic diagrm of bipartition on an arbitrary system.
(b) indicates the system is compacted to a torus. And this bipartition is
trivial so that the boundaries are contractile.

And then, the von Neumann entanglement entropy of P1 for a system can be
represented as follows,

$S({\mathrm{P_{1}}})=\alpha\mathcal{L}-\gamma,$ (14)

where the coefficient $\alpha$ is not universal and depends on the details of
the state, $\mathcal{L}$ is the codimensional-1 boundary of P1, and $\gamma$
is just the universal TEE. As we know, $S({\mathrm{P_{1}}})\geq 0$ for an
arbitrary gapped system [63]. For those systems with $\gamma=0$, it is
possible to deform the ground state to obtain $\alpha=0$, namely one can
remove the leading area-low contribution ($\alpha\mathcal{L}=0$). Thus, the
state with $\alpha=0$ is the pure direct product one without long-range
entanglement, i.e. non-topological gapped phase. On the contrary, the
deformation could never occur for the state with $\gamma>0$, which is the case
for a topological gapped phase. Therefore, the nonzero $\gamma$ is a direct
sign of long-range entanglement.

Numerical calculation of TEE from the von Neumann entropy is difficult, so we
focus on the Renyi entropy here. The Renyi entropy for the gapped state with
bipartition is defined as [64],

$S_{n}=\frac{1}{1-n}\ln\left[\mathrm{Tr}\left(\rho_{1}^{n}\right)\right],$
(15)

where $\rho_{1}$ is the reduced density matrix obtained by tracing out the
subsystem P2, $\rho_{1}=\mathrm{Tr}_{2}|\Psi\rangle\langle\Psi|$, where
$|\Psi\rangle$ is a normalized wave function of the system. In this paper, we
just focus on the Renyi entropy with index $n=2$,
$S_{2}=-\ln\left[\mathrm{Tr}\left(\rho_{1}^{2}\right)\right]$. We define a
swap operator $X$[65] with the purpose as follow,

$X|\alpha_{1}\rangle\otimes|\alpha_{2}\rangle=|\beta_{1}\rangle\otimes|\beta_{2}\rangle,$
(16)

where $|\alpha_{1}\rangle=|a\rangle|b\rangle$ and
$|\alpha_{2}\rangle=|m\rangle|n\rangle$ are two configurations, the
$|a\rangle$ and $|m\rangle$ are in P1 while $|b\rangle$ and $|n\rangle$ are in
P2, and then $|\beta_{1}\rangle=|m\rangle|b\rangle$ and
$|\beta_{2}\rangle=|a\rangle|n\rangle$.

Finally, we can rewrite $S_{2}$ in terms of the expectation of $X$ with
respect to the wave function $|\Psi\rangle\otimes|\Psi\rangle$,
$S_{2}=-\ln\langle X\rangle$. Empirically, $\langle X\rangle$ is predicted to
be a complex number in actual calculation if $|\Psi\rangle$ is complex. So we
can divide this expectation into two parts, $\langle X\rangle=\langle
X_{mod}\rangle\langle X_{phase}\rangle$, which can be individually calculated
by Monte Carlo (MC) method as shown in Eq. 17. It is worth mentioning that
$\tilde{\rho}_{\alpha_{1},\alpha_{2}}$ is a joint probability distribution.
Besides, we note that for large size $\mathcal{L}$, the computational cost is
relatively high because we have to taken more samples in the MC process to
reduce numerical error. Therefore, we calculate the $S(\mathcal{L})$ with
$\mathcal{L}=1$ up to $4$ to eliminate the area contribution and obtain the
TEE.

$\displaystyle\langle
X_{mod}\rangle=\sum_{\alpha_{1},\alpha_{2}}\rho_{\alpha_{1}}\rho_{\alpha_{2}}\left|f(\alpha_{1},\alpha_{2})\right|,$
(17) $\displaystyle\langle
X_{phase}\rangle=\sum_{\alpha_{1},\alpha_{2}}\tilde{\rho}_{\alpha_{1},\alpha_{2}}e^{i\theta(\alpha_{1},\alpha_{2})},$
$\displaystyle\rho_{\alpha_{i}}=\frac{\left|\langle\alpha_{i}|\Psi\rangle\right|^{2}}{\langle\Psi|\Psi\rangle},f(\alpha_{1},\alpha_{2})=\frac{\langle\beta_{1}|\Psi\rangle\langle\beta_{2}|\Psi\rangle}{\langle\alpha_{1}|\Psi\rangle\langle\alpha_{2}|\Psi\rangle},$
$\displaystyle\tilde{\rho}_{\alpha_{1},\alpha_{2}}=\frac{\left|\langle\alpha_{1}|\Psi\rangle\langle\alpha_{2}|\Psi\rangle\right|^{2}\left|f(\alpha_{1},\alpha_{2})\right|}{\sum_{\alpha_{1},\alpha_{2}}\left|\langle\alpha_{1}|\Psi\rangle\langle\alpha_{2}|\Psi\rangle\right|^{2}\left|f(\alpha_{1},\alpha_{2})\right|},$
$\displaystyle
e^{i\theta(\alpha_{1},\alpha_{2})}=\frac{f(\alpha_{1},\alpha_{2})}{\left|f(\alpha_{1},\alpha_{2})\right|}.$

## References

  * Anderson [1973] P. W. Anderson, Mater. Res. Bull. 8, 153 (1973).
  * Wen  _et al._ [1989] X.-G. Wen, F. Wilczek, and A. Zee, Phys. Rev. B 39, 11413 (1989).
  * Wen [1989] X.-G. Wen, Phys. Rev. B 40, 7387 (1989).
  * Wen and Niu [1990] X.-G. Wen and Q. Niu, Phys. Rev. B 41, 9377 (1990).
  * Levin and Wen [2006] M. Levin and X.-G. Wen, Phys. Rev. Lett. 96, 110405 (2006).
  * Kitaev and Preskill [2006] A. Kitaev and J. Preskill, Phys. Rev. Lett. 96, 110404 (2006).
  * Kane and Mele [2005] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005).
  * Bernevig and Zhang [2006] B. A. Bernevig and S. C. Zhang, Phys. Rev. Lett. 96, 106802 (2006).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Wu  _et al._ [2006] C. Wu, B. A. Bernevig, and S. C. Zhang, Phys. Rev. Lett. 96, 106401 (2006).
  * Xu and Moore [2006] C. Xu and J. E. Moore, Phys. Rev. B 73, 045322 (2006).
  * Depenbrock  _et al._ [2012] S. Depenbrock, I. P. McCulloch, and U. Schollwock, Phys. Rev. Lett. 109, 067201 (2012).
  * Ran  _et al._ [2007] Y. Ran, M. Hermele, P. A. Lee, and X.-G. Wen, Phys. Rev. Lett. 98, 117205 (2007).
  * Iqbal  _et al._ [2015] Y. Iqbal, D. Poilblanc, and F. Becca, Phys. Rev. B 91, 020402(R) (2015).
  * He  _et al._ [2017] Y.-C. He, M. P. Zaletel, M. Oshikawa, and F. Pollmann, Phys. Rev. X 7, 031020 (2017).
  * Tay and Motrunich [2011] T. Tay and O. I. Motrunich, Phys. Rev. B 84, 020404(R) (2011).
  * Gong  _et al._ [2015] S.-S. Gong, W. Zhu, L. Balents, and D. N. Sheng, Phys. Rev. B 91, 075112 (2015).
  * Hering and Reuther [2017] M. Hering and J. Reuther, Phys. Rev. B 95, 054418 (2017).
  * Bieri  _et al._ [2015] S. Bieri, L. Messio, B. Bernu, and C. Lhuillier, Phys. Rev. B 92, 060407(R) (2015).
  * Hu  _et al._ [2015] W.-J. Hu, W. Zhu, Y. Zhang, S. Gong, F. Becca, and D. N. Sheng, Phys. Rev. B 91, 041124(R) (2015).
  * He  _et al._ [2014] Y.-C. He, D. N. Sheng, and Y. Chen, Phys. Rev. Lett. 112, 137202 (2014).
  * He and Chen [2015] Y.-C. He and Y. Chen, Phys. Rev. Lett. 114, 037201 (2015).
  * Wietek  _et al._ [2015] A. Wietek, A. Sterdyniak, and A. M. Läuchli, Phys. Rev. B 92, 125122 (2015).
  * Jiang  _et al._ [2008] H.-C. Jiang, Z.-Y. Weng, and D. N. Sheng, Phys. Rev. Lett. 101, 117203 (2008).
  * Yan  _et al._ [2011] S. Yan, D. A. Huse, and S. R. White, Science 332, 1173 (2011).
  * Liao  _et al._ [2017] H. J. Liao, Z. Y. Xie, J. Chen, Z. Y. Liu, H. D. Xie, R. Z. Huang, B. Normand, and T. Xiang, Phys. Rev. Lett. 118, 137202 (2017).
  * Guo and Franz [2009] H.-M. Guo and M. Franz, Phys. Rev. B 80, 113102 (2009).
  * Wen  _et al._ [2010] J. Wen, A. Rüegg, C.-C. J. Wang, and G. A. Fiete, Phys. Rev. B 82, 075125 (2010).
  * Tang  _et al._ [2011] E. Tang, J.-W. Mei, and X.-G. Wen, Phys. Rev. Lett. 106, 236802 (2011).
  * Yu and Li [2012] S.-L. Yu and J.-X. Li, Phys. Rev. B 85, 144402 (2012).
  * Yin  _et al._ [2018] J.-X. Yin, S.-S. Zhang, H. Li, K. Jiang, G. Chang, B. Zhang, B. Lian, C. Xiang, I. Belopolski, H. Zheng, T. A. Cochran, S.-Y. Xu, G. Bian, K. Liu, T.-R. Chang, H. Lin, Z.-Y. Lu, Z. Wang, S. Jia, W. Wang, and M. Z. Hasan, Nature 562, 91 (2018).
  * Ye  _et al._ [2018] L. Ye, M. Kang, J. Liu, F. von Cube, C. R. Wicker, T. Suzuki, C. Jozwiak, A. Bostwick, E. Rotenberg, D. C. Bell, L. Fu, R. Comin, and J. G. Checkelsky, Nature 555, 638 (2018).
  * Liu  _et al._ [2019] D.-F. Liu, A.-J. Liang, E.-K. Liu, Q.-N. Xu, Y.-W. Li, C. Chen, D. Pei, W.-J. Shi, S.-K. Mo, P. Dudin, T. Kim, C. Cacho, G. Li, Y. Sun, L.-X. Yang, Z.-K. Liu, S. S. P. Parkin, C. Felser, and Y.-L. Chen, Science 365, 1282 (2019).
  * Ortiz  _et al._ [2020] B. R. Ortiz, S. M. L. Teicher, Y. Hu, J. L. Zuo, P. M. Sarte, E. C. Schueller, A. M. M. Abeykoon, M. J. Krogstad, S. Rosenkranz, R. Osborn, R. Seshadri, L. Balents, J. He, and S. D. Wilson, Phys. Rev. Lett. 125, 247002 (2020).
  * Chen  _et al._ [2021] H. Chen, H. Yang, B. Hu, Z. Zhao, J. Yuan, Y. Xing, G. Qian, Z. Huang, G. Li, Y. Ye, S. Ma, S. Ni, H. Zhang, Q. Yin, C. Gong, Z. Tu, H. Lei, H. Tan, S. Zhou, C. Shen, X. Dong, B. Yan, Z. Wang, and H. J. Gao, Nature 599, 222 (2021).
  * Nie  _et al._ [2022] L. Nie, K. Sun, W. Ma, D. Song, L. Zheng, Z. Liang, P. Wu, F. Yu, J. Li, M. Shan, D. Zhao, S. Li, B. Kang, Z. Wu, Y. Zhou, K. Liu, Z. Xiang, J. Ying, Z. Wang, T. Wu, and X. Chen, Nature 604, 59 (2022).
  * Zheng  _et al._ [2022] L. Zheng, Z. Wu, Y. Yang, L. Nie, M. Shan, K. Sun, D. Song, F. Yu, J. Li, D. Zhao, S. Li, B. Kang, Y. Zhou, K. Liu, Z. Xiang, J. Ying, Z. Wang, T. Wu, and X. Chen, Nature 611, 682 (2022).
  * Yin  _et al._ [2022] J.-X. Yin, B. Lian, and M. Z. Hasan, Nature 612, 647 (2022).
  * Elhajal  _et al._ [2002] M. Elhajal, B. Canals, and C. Lacroix, Phys. Rev. B 66, 014422 (2002).
  * Kumar  _et al._ [2015] K. Kumar, K. Sun, and E. Fradkin, Phys. Rev. B 92, 094433 (2015).
  * Messio  _et al._ [2017] L. Messio, S. Bieri, C. Lhuillier, and B. Bernu, Phys. Rev. Lett. 118, 267201 (2017).
  * Owerre [2017] S. A. Owerre, Phys. Rev. B 95, 014422 (2017).
  * Zorko  _et al._ [2008] A. Zorko, S. Nellutla, J. van Tol, L. C. Brunel, F. Bert, F. Duc, J. C. Trombe, M. A. de Vries, A. Harrison, and P. Mendels, Phys. Rev. Lett. 101, 026405 (2008).
  * Matan  _et al._ [2011] K. Matan, B. M. Bartlett, J. S. Helton, V. Sikolenko, S. Mat’aš, K. Prokeš, Y. Chen, J. W. Lynn, D. Grohol, T. J. Sato, M. Tokunaga, D. G. Nocera, and Y. S. Lee, Phys. Rev. B 83, 214406 (2011).
  * Lee and Nagaosa [2013] P. A. Lee and N. Nagaosa, Phys. Rev. B 87, 064423 (2013).
  * Jiang  _et al._ [2012] H.-C. Jiang, Z. Wang, and L. Balents, Nat. Phys. 8, 902 (2012).
  * Sato and Ando [2017] M. Sato and Y. Ando, Rep. Prog. Phys. 80, 076501 (2017).
  * Levin and Wen [2005] M. A. Levin and X.-G. Wen, Phys. Rev. B 71, 045110 (2005).
  * Scharfenberger  _et al._ [2011] B. Scharfenberger, R. Thomale, and M. Greiter, Phys. Rev. B 84, 140404 (2011).
  * Neupert  _et al._ [2014] T. Neupert, C. Chamon, C. Mudry, and R. Thomale, Phys. Rev. B 90, 205101 (2014).
  * Affleck  _et al._ [1988] I. Affleck, Z. Zou, T. Hsu, and P. W. Anderson, Phys. Rev. B 38, 745 (1988).
  * Wen [2002] X.-G. Wen, Phys. Rev. B 65, 165113 (2002).
  * Lu  _et al._ [2011] Y.-M. Lu, Y. Ran, and P. A. Lee, Phys. Rev. B 83, 224413 (2011).
  * Bieri  _et al._ [2016] S. Bieri, C. Lhuillier, and L. Messio, Phys. Rev. B 93, 094437 (2016).
  * Sorella and Capriotti [2000] S. Sorella and L. Capriotti, Phys. Rev. B 61, 2599 (2000).
  * Iqbal  _et al._ [2016] Y. Iqbal, W.-J. Hu, R. Thomale, D. Poilblanc, and F. Becca, Phys. Rev. B 93, 144411 (2016).
  * Zhao and Liu [2021] Q.-R. Zhao and Z.-X. Liu, Phys. Rev. Lett. 127, 127205 (2021).
  * He and Li [2023] L.-W. He and J.-X. Li, Phys. Rev. B 108, 245102 (2023).
  * Wen [1990] X. G. Wen, Int. J. Mod. Phys. B 04, 239 (1990).
  * Wen [2015] X.-G. Wen, Natl. Sci. Rev. 3, 68 (2015).
  * 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, Science 340, 167 (2013).
  * Fukui  _et al._ [2005] T. Fukui, Y. Hatsugai, and H. Suzuki, J. Phys. Soc. Jpn. 74, 1674 (2005).
  * Savary and Balents [2017] L. Savary and L. Balents, Rep. Prog. Phys. 80, 016502 (2017).
  * Zhang  _et al._ [2011] Y. Zhang, T. Grover, and A. Vishwanath, Phys. Rev. B 84, 075128 (2011).
  * Hastings  _et al._ [2010] M. B. Hastings, I. Gonzalez, A. B. Kallin, and R. G. Melko, Phys. Rev. Lett. 104, 157201 (2010).


# Policy Gradient RL Algorithms as Directed Acyclic Graphs

Juan Jose Garau Luis
Massachusetts Institute of Technology
Department of Aeronautics and Astronautics
<EMAIL_ADDRESS>

###### Abstract

Meta Reinforcement Learning (RL) methods focus on automating the design of RL
algorithms that generalize to a wide range of environments. The framework
introduced in [1] addresses the problem by representing different RL
algorithms as Directed Acyclic Graphs (DAGs), and using an evolutionary meta
learner to modify these graphs and find good agent update rules. While the
search language used to generate graphs in the paper serves to represent
numerous already-existing RL algorithms (e.g., DQN, DDQN), it has limitations
when it comes to representing Policy Gradient algorithms. In this work we try
to close this gap by extending the original search language and proposing
graphs for five different Policy Gradient algorithms: VPG, PPO, DDPG, TD3, and
SAC.

##  1 Introduction

Novel Reinforcement Learning (RL) algorithms are recurrently being proposed to
address a wide range of performance gaps. Each of these algorithms involves a
human-driven design process that does not scale well. In order to reduce the
design cost, the research community is focusing on meta RL algorithms, which
automate the process of developing new policy update rules that generalize to
multiple unseen environments.

Broadly, meta RL methods consider a double-loop algorithm in which an agent
attempts to learn good policies for arbitrary environments in an inner loop
and a meta learner attempts to produce good update rules in an outer loop. The
inner agent uses a loss function provided by the meta learner, whereas the
meta learner uses the agent’s return to guide its search.

Among all the meta RL methods that are currently being studied, in this work
we focus on the framework presented in [1], in which authors introduce a
search language with 26 operators used by the meta learner to form different
update rules by connecting multiple operators and constructing a computation
graph in the form of a Directed Acyclic Graph (DAG). The meta learner explores
the space of graphs by means of evolutionary strategies and keeps track of a
population with the best loss functions. The goal of the meta learner is to
find the best RL algorithm given by loss function $L^{*}$, defined as

$L^{*}=\underset{L}{\text{arg
max}}\left[\sum_{\varepsilon}\text{Eval}(L,\varepsilon)\right]$ (1)

where each $\varepsilon$ belongs to a set of different training environments,
all of them used to evaluate the agent’s performance, defined in turn as

$\text{Eval}(L,\varepsilon)=\frac{1}{M_{\varepsilon}}\sum_{m=1}^{M_{\varepsilon}}\frac{R_{m}-R_{min}}{R_{max}-R_{min}}$
(2)

where $R_{m}$ corresponds to the return during episode $m$, $M_{\varepsilon}$
is the total number of episodes using environment $\varepsilon$, and $R_{min}$
and $R_{max}$ are the minimum and maximum possible returns for the
environment, respectively. This method can start a search from scratch or use
existing algorithms as warm-starts.

The authors test their approach using a set of training environments
consisting of 4 classical control tasks from OpenAI Gym [2] and 12 multitask
gridworld styile environments from MiniGrid [3]. They use 4 of these
environments to meta train and 12 to meta test, and the DQN loss function [4]
as an input algorithm. They find two improved versions of the DQN loss
function; one of them, defined as DQNReg, outperforms both DQN and DDQN [5] in
15 out of the 16 environments. DQNReg is then tested on 4 Atari environments
[6], where it is also shown to outperform PPO [7].

One of the salient conclusions of the paper is that representing RL algorithms
as DAGs is a good way to obtain interpretable solutions. The authors reinforce
this idea by including a comprehensive analysis of two of the novel update
rules found by their framework based on their graph structure. While it is a
promising research direction, the search language proposed in the original
paper is not complete enough to build on other types of RL algorithms, such as
Policy Gradient algorithms [8]. In this work we try to close this research gap
by proposing extensions to their framework and representing five different
Policy Gradient algorithms as DAGs: VPG, PPO, DDPG, TD3, and SAC.

##  2 Related work

Representing algorithms as computational graphs is also proposed in other
studies. In [9] the same concept is used to find appropriate image
classifiers. This idea is also present in [10], where a meta learner produces
curiosity graphs that drive the agent’s exploration by shaping the
environment’s reward. The agent then uses the modified reward with a fixed
update rule. Conversely, the evolutionary strategy in [1] allows flexible
update rules and is meta trained on a diverse set of environments.

As oposed to changing the form of RL algorithms, other studies consider a
fixed-form update rule which takes in parameters provided by the meta learner.
It is the case of [11], where authors propose an evolutionary strategy to
learn good policy update parameters $\hat{\pi}$ that guide the policy
optimization. They train the meta RL framework on several MuJoCo tasks [12]
and test on similar versions of those tasks. The same train and test procedure
is followed by the authors in [13], although a meta-gradient descent strategy
is adopted instead of evolutionary methods. This approach is also followed by
the authors in [14], which propose MetaGenRL to attain higher generalization
capabilities. Those are proven in test environments based on unseen MuJoCo
tasks.

In addition to learning appropriate policy update rules $\hat{\pi}$, the
method presented in [15], defined as Learned Policy Gradient, also learns a
prediction update rule $\hat{y}$, i.e., the semantics of the agent’s
prediction. In some experiments, the discovered semantics converge towards a
notion of value function. The authors claim LPG is able to attain better
generalizability this way, as proven on the performance on unseen Atari test
environments after meta training on toy environments.

##  3 Methods

The focus of this work is to extend the applicability of the framework
presented in [1] by adapting the search language to Policy Gradient algorithms
[8]. First, we present a high-level description of the method from the
original paper, summarized in Figure 1. Then, we introduce our extensions to
the search language.

Figure 1: Evolutionary Reinforcement Learning method [1] overview. The meta
learning framework constructs a population of RL algorithms or loss functions
from a predefined search language, randomly creating each one or using an
already-existing algorithm as a warm-start. For a fixed number of iterations,
the best RL algorithm from a random set of $T$ algorithms undergoes mutation
and a child algorithm is generated. This new algorithm is evaluated, first on
a hurdle environment, and then on a set of training environments. The meta
learner collects the returns from each environment and computes the score of
the new algorithm. The oldest algorithm from the population is then removed.

###  3.1 General Description

The input to the algorithm consists of a set of training environments, one
additional “hurdle environment”, and a search language that conforms the nodes
of the computational graphs. The algorithm starts by creating a population of
$N$ 111This symbol is defined by us, it is not in the original paper RL
algorithms or loss functions, which can be constructed from scratch using the
search language or from an input already-existing RL algorithm. The starting
population is scored in each of the training environments.

For a defined number of $C$ iterations, the meta learner first samples a set
of $T$ algorithms from the population and selects the one with the highest
score as a parent algorithm. Then, a child algorithm is created by mutating
the parent algorithm. The hurdle environment is used to assess a minimum
performance threshold $\alpha$, which determines whether the child algorithm
is considered for the population or not. If it passes the threshold, the child
algorithm is evaluated on each of the training environments by independent
agents, and an aggregated score is generated. The new algorithm is added to
the population, which is kept to a constant size of $N$ individuals by
removing the oldest one.

###  3.2 Search Language Extensions

The framework expresses RL algorithms as DAGs of nodes with typed inputs and
outputs. The input nodes consist of the elements in the tuple
$(s_{t},a_{t},r_{t},s_{t+1})$, the network parameters $\theta$, the target
network parameters $\theta^{\prime}$, the discount factor $\gamma$, and other
numerical constants (e.g., 1, 0.1). Other types of nodes include parameter
nodes, which encode neural network operations, and operation nodes, which
perform scalar and list operation over the inputs (e.g., Add, MaxList). A full
list of the available nodes can be found in the appendix of the original paper
[1]. For reference, Figure 2 shows the DDQN algorithm [5] encoded as a
computational graph using the _original_ set of nodes from the paper.

Figure 2: Computational graph representing DDQN loss function
$L_{\theta}=\left(Q_{\theta}\left(s_{t},a_{t}\right)-\left(r_{t}+\gamma
Q_{\theta^{\prime}}\left(s_{t+1},\underset{a}{\text{arg
max}}Q_{\theta}\left(s_{t+1},a\right)\right)\right)\right)^{2}$.

The output of the computational graph, given by the node Output, is then used
as the loss function $L_{\theta}$ the agent must minimize, updating its
parameters by gradient descent
$\theta\longleftarrow\theta-\alpha_{lr}\nabla_{\theta}L_{\theta}$. We follow
this gradient descent convention throughout the rest of this work.

We now introduce each of the extensions to the original language, divided into
new input nodes, new operation nodes, and new constant nodes. We follow the
data type notation convention proposed by the authors, which includes state
$\mathbb{S}$, action $\mathbb{Z}$, float $\mathbb{R}$, and list
$List[\mathbb{X}]$, where $\mathbb{X}$ indicates it can be of $\mathbb{S}$ or
$\mathbb{R}$222Notation extracted from the original paper [1].

####  3.2.1 Input nodes

In the paper it is assumed each independent agent initializes Q-value
parameters $\theta$, target network parameters $\theta^{\prime}$, and an empty
replay buffer $\mathcal{D}$, where the agent stores all the transition tuples
encountered in order to do off-policy learning. To adapt the framework to
Policy Gradient algorithms, we propose the following extensions:

  * •

We register the episode-termination signal $d_{t}$ in the replay buffer
$\mathcal{D}$.

  * •

We add the parameter $\lambda$ [8] to the language.

  * •

We allow the tuple $(s_{t},a_{t},r_{t},d_{t},s_{t+1})$ to represent a sequence
of consecutive timesteps, instead of a single timestep. According to the
original paper, operations broadcast in the cases in which inputs have a
different number of dimensions.

  * •

We use additional network parameters $\phi$ (e.g., for the value network when
needed) and their target counterparts $\phi^{\prime}$.

  * •

In addition to the network parameter nodes from the original paper, which
allow mappings from states $\mathbb{S}$ to floats $\mathbb{R}$ and lists
$List[\mathbb{R}]$, we also allow to map both a state $\mathbb{S}$ and an
action $\mathbb{Z}$ to a float $\mathbb{R}$, i.e.,
$\theta:\mathbb{S},\mathbb{Z}\longrightarrow\mathbb{R}$. Likewise, we
introduce the parameter node
$NN:\mathbb{S},\mathbb{Z}\longrightarrow\mathbb{R}$ to carry out operations of
such kind. We do this in order to address _continuous_ state and action
spaces, as required by certain algorithms, such as DDPG [16].

####  3.2.2 Operation nodes

Apart from the operation nodes proposed in the original paper, we introduce
the additional nodes:

  * •

SumAndDiscount: It takes in a vector $v$ of length $l_{v}$ and a numerical
value $b$, and outputs a vector $v^{\prime}$ of length $l_{v}$ where each
element is computed as $v^{\prime}_{i}=\sum_{k=i}^{l_{v}}v_{k}\cdot b^{k-i}$.
This is necessary to compute certain elements of different Policy Gradient
methods, such as computing the rewards-to-go (in that case, $v=r$ and $b=1$).

  * •

Clip: It clips a scalar or vector component-wise given minimum and maximum
float values.

  * •

Squashing: Given a state $s_{t}$ from data type $\mathbb{S}$ and network
parameters $\theta:\mathbb{S}\longrightarrow\mathbb{R}$, it computes actions
$\tilde{a}$ following a squashed Gaussian policy, i.e.,
$\tilde{a}=\tanh(\mu_{\theta}(s_{t})+\sigma_{\theta}(s_{t})\odot\xi)$, with
$\xi\sim\mathcal{N}(0,I)$. This node is specifically introduced to represent
the Soft Actor-Critic algorithm [17].

  * •

Prob: In cases in which the action space is continuous, and therefore the node
Softmax can not be used, this node computes the probability of taking a
certain action in a given state with given network parameters $\theta$. It
assumes the probabilities can be computed, such as in the case of Gaussian or
categorical policies.

####  3.2.3 Constant nodes

Finally, we encode useful constant values as new constant nodes for the
language. These are general values, such as $-1$, or algorithm-specific
constants, such as $1+\epsilon$ and $1-\epsilon$, for the example case of PPO.
We introduce all algorithm-specific constant nodes in their related sections.

##  4 Example Graphs

Based on the language extensions introduced in the previous section, we
propose computational graphs in the form of DAGs for the following Policy
Gradient algorithms: Vanilla Policy Gradient (VPG), Proximal Policy
Optimization (PPO) [7], Deep Deterministic Policy Gradient (DDPG) [16], Twin
Delayed DDPG (TD3) [18], and Soft Actor-Critic (SAC) [17]. We base our graph
constructions on the descriptions and implementations provided in [19]. For
all algorithms, we leave the target parameters update considerations (i.e.,
$\theta^{\prime}\longleftarrow\theta$, $\phi^{\prime}\longleftarrow\phi$) for
the framework user to decide.

###  4.1 Vanilla Policy Gradient (VPG)

The work in [1] only considers off-policy algorithms that sample batches of
tuples $(s_{t},a_{t},r_{t},s_{t+1})$. In this work we also allow $s_{t}$ to
represent a sequence of consecutive states. We apply this idea to construct a
computational graph for the VPG algorithm that can be used as an input to the
meta learning framework. We need to consider two different loss functions: one
for the policy $\pi_{\theta}$ and another for the value function $V_{\phi}$.

First, Figure 3 shows the computational graph for the policy loss, in which we
assume the use of GAE-Lambda advantage estimation method [8] to compute the
discounted advantage for each timestep. In this example, we introduce the
novel node SumAndDiscount and a constant node with the value $-1$, in order to
follow the loss minimization convention from the original paper.

Figure 3: Computational graph representing Vanilla Policy Gradient (VPG)
policy $\pi_{\theta}$ loss
$L_{\theta}=-\frac{1}{T}\sum_{t=0}^{T}\log\pi_{\theta}\left(a_{t}|s_{t}\right)\hat{A}_{t}$.

Then, Figure 4 shows the computational graph for the value function loss,
which computes the MSE between the value function and the rewards-to-go.

Figure 4: Computational graph representing Vanilla Policy Gradient (VPG) value
function $V_{\phi}$ loss
$L_{\phi}=\frac{1}{T}\sum_{t=0}^{T}\left(V_{\phi}(s_{t})-\hat{R}_{t}\right)^{2}$.

###  4.2 Proximal Policy Optimization (PPO)

In the case of the PPO algorithm we consider the policy parameters $\theta$ as
well as an early version $\theta_{k}$ to compute the policy change ratio.
Figure 5 shows the computational graph for this algorithm, which includes the
novel node Clip, as well as the constant nodes $1+\epsilon$ and $1-\epsilon$,
input to the clipping function. We use the same advantage estimation method
used in the VPG case, and also consider the same value function update method,
already presented in Figure 4.

Figure 5: Computational graph representing Proximal Policy Optimization (PPO)
policy $\pi_{\theta}$ loss
$L_{\theta}=-\frac{1}{T}\sum_{t=0}^{T}\min\left(\frac{\pi_{\theta}(a_{t}|s_{t})}{\pi_{\theta_{k}}(a_{t}|s_{t})}\hat{A}_{t},\text{clip}\left(\frac{\pi_{\theta}(a_{t}|s_{t})}{\pi_{\theta_{k}}(a_{t}|s_{t})},1-\epsilon,1+\epsilon\right)\hat{A}_{t}\right)$.

###  4.3 Deep Deterministic Policy Gradient (DDPG)

Next, we examine the DDPG algorithm, which in the original paper [16] was
proposed to address continuous action spaces. As explained in the previous
section, we adapt the framework to satisfy continuity requirements by
introducing the parameter node
$NN:\mathbb{S},\mathbb{Z}\longrightarrow\mathbb{R}$, which accounts for a
Q-function taking a state and a continuous action as inputs. DDPG algorithm
considers both a policy $\mu_{\theta}$ and a Q-function $Q_{\phi}$, as well as
their target counterparts, $\mu_{\theta^{\prime}}$ and $Q_{\phi^{\prime}}$,
respectively. Figure 6 shows the computational graph to update the Q-function
parameters $\phi$, which makes use of the episode-termination signal $d_{t}$,
also introduced in this work.

Figure 6: Computational graph representing Deep Deterministic Policy Gradient
(DDPG) Q-function $Q_{\phi}$ loss
$L_{\phi}=\left(Q_{\phi}\left(s_{t},a_{t}\right)-\left(r_{t}+\gamma\left(1-d_{t}\right)Q_{\phi^{\prime}}\left(s_{t+1},\mu_{\theta^{\prime}}\left(s_{t+1}\right)\right)\right)\right)^{2}$.

Then, Figure 7 shows the computational graph to update policy parameters
$\theta$.

Figure 7: Computational graph representing Deep Deterministic Policy Gradient
(DDPG) policy $\mu_{\theta}$ loss
$L_{\theta}=-Q_{\phi}\left(s_{t},\mu_{\theta}\left(s_{t}\right)\right)$.

###  4.4 Twin Delayed DDPG (TD3)

The TD3 loss functions can be constructed with all input, parameter, and
operation nodes presented so far. However, new constant nodes are needed for
these functions. Figure 8 shows the computational graph to update the
Q-function parameters $\phi_{i}$. In the TD3 algorithm, clipped noise is added
to the actions as a regularizer. In the graph, this is achieved by using input
nodes with constant values $-c$ and $c$, as well as a constant $\epsilon$
sampled from the normal distribution $\mathcal{N}(0,\sigma)$. To make sure
actions are kept within bounds, additional constant nodes with action bounds
$a_{low}$ and $a_{high}$ are used.

Figure 8: Computational graph representing Twin Delayed DDPG (TD3) Q-function
$Q_{\phi_{i}}$ loss, for $i=1,2$, given by
$L_{\phi_{i}}=\left(Q_{\phi_{i}}\left(s_{t},a_{t}\right)-\left(r_{t}+\gamma\left(1-d_{t}\right)\underset{i=1,2}{\min}Q_{\phi^{\prime}_{i}}\left(s_{t+1},a\left(s_{t+1}\right)\right)\right)\right)^{2}$,
where
$a(s_{t+1})=\text{clip}(\mu_{\theta^{\prime}}(s_{t+1})+\text{clip}(\epsilon,-c,c),a_{\text{low}},a_{\text{high}})$,
with $\epsilon\sim\mathcal{N}(0,\sigma)$.

Note this graph corresponds to generic Q-function parameters $\phi_{i}$; TD3
algorithm considers two different Q-functions, $Q_{\phi_{1}}$ and
$Q_{\phi_{2}}$, the framework user can decide whether to use the same graph to
update both or use two different DAGs. In contrast, only Q-function
$Q_{\phi_{1}}$ is used to update policy parameters $\theta$, as shown in
Figure 9.

Figure 9: Computational graph representing Twin Delayed DDPG (TD3) policy
$\pi_{\theta}$ loss $L_{\theta}=-Q_{\phi_{1}}(s_{t},\mu_{\theta}(s_{t}))$.

###  4.5 Soft Actor Critic (SAC)

Finally, we address SAC loss functions. Following the implementation in [19],
we assume the use of environments with continuous action spaces. Figure 10
shows the computational graph corresponding to the generic Q-function
$Q_{\phi_{i}}$ – SAC uses two different Q-functions following the same
rationale as TD3. In this case we use the two last operation nodes introduced
in this work, Squashing and Prob. We follow the convention from the original
paper of considering single-output nodes. In case the framework user wants to
consider multiple-output nodes, the Squashing and Prob nodes could be
integrated in a single node. To construct this graph, we also use an
additional constant node with hyperparameter $\alpha$ [17].

Figure 10: Computational graph representing Soft Actor-Critic (SAC) Q-function
$Q_{\phi}$ loss
$L_{\phi}=\left(Q_{\phi_{i}}\left(s_{t},a_{t}\right)-\left(r_{t}+\gamma\left(1-d_{t}\right)\underset{i=1,2}{\min}Q_{\phi^{\prime}_{i}}\left(s_{t+1},\tilde{a}\right)-\alpha\log_{\pi_{\theta}}\left(\tilde{a}|s_{t+1}\right)\right)\right)^{2}$,
where
$\tilde{a}=\tanh(\mu_{\theta}(s_{t+1})+\sigma_{\theta}(s_{t+1})\odot\xi)$,
with $\xi\sim\mathcal{N}(0,I)$.

Figure 11 shows the computational graph to update policy parameters $\theta$,
which makes use of the same new nodes.

Figure 11: Computational graph representing Soft Actor-Critic (SAC) policy
$\pi_{\theta}$ loss
$L_{\theta}=-\left(\underset{i=1,2}{\min}Q_{\phi_{i}}(s_{t},\tilde{a})-\alpha\log_{\pi_{\theta}}(\tilde{a}|s_{t})\right)$,
where $\tilde{a}=\tanh(\mu_{\theta}(s_{t})+\sigma_{\theta}(s_{t})\odot\xi)$,
with $\xi\sim\mathcal{N}(0,I)$.

##  5 Conclusions

In this work we have presented an extension of the search language introduced
in [1], which serves as a set of building blocks to generate multiple RL
algorithms in the form of Directed Acyclic Graphs. Specifically, we have
focused on adapting the language to Policy Gradient algorithms and, to that
end, proposed four additional operation nodes as well as five modifications to
the input structure of the graphs. By means of these extensions, we have shown
computational graphs for five different Policy Gradient algorithms: VPG, PPO,
DDPG, TD3, and SAC.

Overall, representing RL algorithms as graphs offers high interpretability to
the framework users, especially when RL algorithms are novel. As part of the
future work, these graphs should be tested as warm starts to the original
evolutionary framework and the resulting child graphs validated in different
unseen environments. In addition, other RL algorithms not addressed in this
work could be represented in the same format.

##  6 Github Repository

A Github repository with all of the figures shown in this work can be found at
https://github.com/jjgarau/DAGPolicyGradient.

## References

  * [1] John D. Co-Reyes, Yingjie Miao, Daiyi Peng, Esteban Real, Sergey Levine, Quoc V. Le, Honglak Lee, and Aleksandra Faust.  Evolving reinforcement learning algorithms, 2021.
  * [2] Greg Brockman, Vicki Cheung, Ludwig Pettersson, Jonas Schneider, John Schulman, Jie Tang, and Wojciech Zaremba.  Openai gym, 2016.
  * [3] Maxime Chevalier-Boisvert, Lucas Willems, and Suman Pal.  Minimalistic gridworld environment for openai gym.  https://github.com/maximecb/gym-minigrid, 2018.
  * [4] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A Rusu, Joel Veness, Marc G Bellemare, Alex Graves, Martin Riedmiller, Andreas K Fidjeland, Georg Ostrovski, et al.  Human-level control through deep reinforcement learning.  nature, 518(7540):529–533, 2015.
  * [5] Hado van Hasselt, Arthur Guez, and David Silver.  Deep reinforcement learning with double q-learning, 2015.
  * [6] M. G. Bellemare, Y. Naddaf, J. Veness, and M. Bowling.  The arcade learning environment: An evaluation platform for general agents.  Journal of Artificial Intelligence Research, 47:253–279, 2013\.
  * [7] John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov.  Proximal policy optimization algorithms, 2017.
  * [8] Richard S Sutton and Andrew G Barto.  Reinforcement learning: An introduction.  MIT press, 2018.
  * [9] Esteban Real, Alok Aggarwal, Yanping Huang, and Quoc V Le.  Regularized evolution for image classifier architecture search, 2019.
  * [10] Ferran Alet, Martin F. Schneider, Tomas Lozano-Perez, and Leslie Pack Kaelbling.  Meta-learning curiosity algorithms, 2020.
  * [11] Rein Houthooft, Richard Y. Chen, Phillip Isola, Bradly C. Stadie, Filip Wolski, Jonathan Ho, and Pieter Abbeel.  Evolved policy gradients, 2018.
  * [12] E. Todorov, T. Erez, and Y. Tassa.  Mujoco: A physics engine for model-based control.  2012 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 5026–5033, 2012.
  * [13] Sarah Bechtle, Artem Molchanov, Yevgen Chebotar, Edward Grefenstette, Ludovic Righetti, Gaurav Sukhatme, and Franziska Meier.  Meta-learning via learned loss, 2020.
  * [14] Louis Kirsch, Sjoerd van Steenkiste, and Jürgen Schmidhuber.  Improving generalization in meta reinforcement learning using learned objectives, 2020.
  * [15] Junhyuk Oh, Matteo Hessel, Wojciech M. Czarnecki, Zhongwen Xu, Hado van Hasselt, Satinder Singh, and David Silver.  Discovering reinforcement learning algorithms, 2020.
  * [16] Timothy P Lillicrap, Jonathan J Hunt, Alexander Pritzel, Nicolas Heess, Tom Erez, Yuval Tassa, David Silver, and Daan Wierstra.  Continuous control with deep reinforcement learning.  arXiv preprint arXiv:1509.02971, 2015.
  * [17] Tuomas Haarnoja, Aurick Zhou, Pieter Abbeel, and Sergey Levine.  Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor, 2018.
  * [18] Scott Fujimoto, Herke van Hoof, and David Meger.  Addressing function approximation error in actor-critic methods, 2018\.
  * [19] Joshua Achiam.  Spinning Up in Deep Reinforcement Learning.  2018\.