# Homotopy of periodic $2\times 2$ matrices

Joseph E. Avron  Department of Physics, Technion, Haifa, Israel  Ari M. Turner
Department of Physics, Technion, Haifa, Israel

###### Abstract

We describe the homotopy classes of loops in the space of $2\times 2$ simple
(=non-degenerate) matrices with various symmetries. This turns out to be an
elementary exercise in the homotopy of closed curves in
$\mathbbm{R}^{3}/\\{0\\}$. Since closed curves in $\mathbbm{R}^{3}/\\{0\\}$
can be readily visualized, no advanced tools of algebraic topology are needed.
The matrices represent gapped Bloch Hamiltonians in 1D with a two dimensional
Hilbert space per unit cell.

##  1 Introduction

We study the homotopy classes of loops in the space of simple $2\times 2$
matrices that satisfy various unitary or anti-unitary symmetries. This is
motivated by the program of classifying the phases of insulating quantum
matter [1, 2, 3, 4]. The classification is enriched by allowing for phases
with symmetry [5, 6, 7, 8].

An insulator is characterized by a (many body) ground state separated from the
excited states by a spectral gap which is bounded away from zero in the limit
of large systems. The ground states of two insulators belong to the same phase
if it is possible to interpolate between the two ground states without closing
the gap. One can picture a phase diagram in a space of Hamiltonians with all
possible parameters; the problem of finding phases can be described
topologically as finding the connected components of the region corresponding
to Hamiltonians with a gap111 Further topological properties of this space are
also interesting; for example, its fundamental group plays a role in quantized
pumping [9] and topologically protected quantum gates [10].. This is a
challenging mathematical program, especially for interacting systems, as it
concerns many-body quantum mechanics in the thermodynamic limit [11, 12, 13,
14, 15, 16, 17, 18].

The phases of topological insulators are of interest for the field of
condensed matter because of their mobile edge states that are impervious to
disorder [19, 20] and the search for quasi-particles with exotic statistics
[10, 21, 22]. They are also of interest in quantum information theory because
of their connections with computational complexity of many-body systems [23].

Arguably the simplest setting for studying phases of insulators is in the
context of Bloch Hamiltonians222For the study of models of free Fermions in
disordered systems see e.g. [24, 25]. for non-interacting Fermions [1, 2].
Because Bloch Hamiltonians have discrete translation symmetry, they can be
studied directly in the limit of infinite systems. If the $m$ lowest bands
(for some $m$) are separated from the rest by a gap, an insulator with $m$
Fermions per period can be formed by filling these states. We focus on the
simplest case where $m=1$ (with spinless Fermions) and where the dimension of
space is one.

In one dimension, a Bloch band gives rise to a loop $P(k)$ of orthogonal
projections where $k$ takes values in the Brillouin zone (BZ). The loop of
band projections corresponds to a notion of a quantum state of a system with
infinitely many particles.

Suppose that one restricts oneself to Bloch bands with a given symmetry, e.g.
time-reversal. Two Bloch bands with projections $P_{0}(k)$ and $P_{1}(k)$ are
said to be homotopic if there is a family of projections $P_{t}(k)$, jointly
continuous333The continuity follows from the locality in real space of the
Hamiltonian and the gap condition. in $(t\in[0,1],k\in\text{BZ})$, respecting
the symmetry, that interpolates444An alternative definition of a homotopy uses
the space of loops which has the topology known as the compact-open topology:
a homotopy is a continuous path in this space.[26] between $P_{0}(k)$ and
$P_{1}(k)$. This corresponds to continuous changes between many-body quantum
states of free Fermions. Describing the equivalence classes of Bloch bands
with various symmetries is an interesting and non-trivial theoretical problem.

Our modest aim is to give an elementary invitation to this field by describing
the homotopic classification of non-degenerate two-band Bloch Hamiltonians in
one dimension with various symmetries. The results, summarised in Table 1, go
somewhat beyond the first column in the periodic table of topological
insulators555A comparison is made in Appendix B.[3, 7, 8] for the Altland-
Zirnbauer symmetry classes [5] as we allow for certain space symmetries and
some results are specific to two bands. Our methods are elementary, as we
reduce the problem to the study of the homotopy of curves in three dimensions.
The elementary nature of the problem is a consequence of:

  * •

The space of one dimensional projections in $\mathbbm{C}^{2}$ is easy to
visualize. It is simply the two sphere $\mathbbm{S}^{2}$. In contrast, the
space of $n\times n$ rank $m$ projections,666This space is the Grassmannian,
$Gr_{m}(\mathbbm{C}^{n})$, a basic object in algebraic geometry and topology.
which describes insulators of $m$ occupied bands (see sec. 6), is difficult to
visualize.

  * •

One-dimensional Bloch Hamiltonians with 2 gapped bands can be represented by
closed curves in $\mathbbm{R}^{3}/\\{0\\}$ hence questions regarding the
homotopy of Bloch Hamiltonians reduce to questions about the homotopy of
curves in 3 dimensional Euclidean space with the origin removed. The study of
the homotopy classes of curves in two and three dimensional spaces is the most
intuitive part of homotopy theory. This is in contrast with the study of the
homotopy of Bloch Hamiltonians in higher dimensions which translates to the
study of the homotopy of higher dimensional closed surfaces which are not
easily visualized and require more advanced tools777This is true even in the
case of two gapped bands, since the high homotopy groups of the 2-sphere are
non-trivial..

Symmetry | Constraint | Homotopy classes
---|---|---
None | None | $\\{\sigma_{z}\\}$
$\Theta_{+}$ | $H(k)=H^{*}(-k)$ | $\\{\sigma_{z}\\}$
$\Theta_{-}$ | $H(k)=\sigma_{y}H^{*}(-k)\sigma_{y}$ | —
$C_{+}$ | $H(k)=-{\sigma_{z}}H^{*}(-k){\sigma_{z}}$ | $\left\\{\pm{\sigma_{x}},\ \pm R_{1}\right\\}$
$C_{-}$ | $H(k)=-{\sigma_{y}}H^{*}(-k){\sigma_{y}}$ | $\\{\sigma_{z}\\}$
B | $H(k)={\sigma_{x}}H(-k){\sigma_{x}}$ | $\\{\pm{\sigma_{x}},\pm R_{1}\\}$
S | $H(k)=G_{k}H(-k)G^{*}_{k}$ | $\\{\pm\sigma_{z}\\}$
$C_{+}\circ\Theta_{+}$ | $H(k)=-{\sigma_{z}}H(k){\sigma_{z}}$ | $\displaystyle{\\{R_{n}|n\in\mathbbm{Z}\\}}$
$B\circ\,\Theta_{+}$ | $H(k)={\sigma_{x}}H^{*}(k){\sigma_{x}}$ | $\\{R_{n}|n\in\mathbbm{Z}\\}$
$S\circ\,\Theta_{+}$ | $H(k)=G_{k}H^{*}(k)G^{*}_{k}$ | $\\{\pm\sigma_{z}\\}$
S $\land\,\Theta_{+}$ | $H(k)=G_{k}H(-k)G^{*}_{k}=H^{*}(-k)$ | $\\{\pm\sigma_{z}\\}$
B $\land\,\Theta_{+}$ | $H(k)={\sigma_{x}}H(-k){\sigma_{x}}=H^{*}(-k)$ | $\displaystyle{\\{\pm R_{n}|n\in\mathbbm{Z}\\}}$
$C_{+}\land\Theta_{+}$ | $H(k)=-{\sigma_{z}}H^{*}(-k){\sigma_{z}}=H^{*}(-k)$ | $\displaystyle{\\{\pm R_{n}|n\in\mathbbm{Z}\\}}$
$C_{-}\land\Theta_{+}$ | $H(k)=-{\sigma_{y}}H^{*}(-k){\sigma_{y}}=H^{*}(-k)$ | $\\{\sigma_{y}\\}$
Table 1: In the first column $\Theta_{\pm}$ stands for time-reversal with
$\Theta_{\pm}^{2}=\pm 1$, similarly for Particle-Hole symmetry $C_{\pm}$. $B$
stands for bond-reflection and $S$ for site reflection. $\land$ means that
both symmetries are present and $\circ$ that only the product is a symmetry.
$C_{+}\circ\Theta_{+}$ is chiral symmetry. In the second column,
${\sigma_{x}},{\sigma_{y}}$ and ${\sigma_{z}}$ denote the Pauli matrices, $*$
is complex conjugation, and $G_{k}$ is the diagonal unitary of Eq. 2.10. The
third column gives the equivalence classes $2\times 2$ periodic gapped
matrices. $\\{\sigma_{x}\\},\\{\sigma_{y}\\},$ and $\\{\sigma_{z}\\}$
represent the trivial (= null-homotopic) Bloch matrix associated with a
lattice of disconnected unit cells. The choice of Pauli matrix is not always
canonical (see the relevant sections). $R_{n}$ is the unitary given in Eq.
4.4. The entries above the top dashed line correspond to generating
symmetries. The entries below show few composite symmetries. The list of
composites is incomplete: Only products with $\Theta_{+}$ are listed.

##  2 Bloch Hamiltonians with two bands

Consider a one-dimensional, single particle Hamiltonian acting on the Hilbert
space

$\ell^{2}(\mathbbm{Z})\otimes\mathbbm{C}^{2}$ (2.1)

This represents a chain of unit cells, each hosting a two dimensional Hilbert
space, e.g. a spin or a pair of atoms; see Fig. 1. The basis states are
$\left|j,a\right\rangle$ labeled by $j\in\mathbbm{Z}$ and $a\in\\{0,1\\}$.

Figure 1: A periodic array with two atoms per unit cell. The red ellipse shows
a unit cell.

A one-particle lattice Hamiltonian is a self-adjoint infinite dimensional
matrix with elements

$\left\langle j,b\right|H\left|i,a\right\rangle=\left\langle
i,a\right|H\left|j,b\right\rangle^{*}$ (2.2)

A Bloch Hamiltonian is, by definition, translation invariant, i.e.

$\left\langle j,b\right|H\left|i,a\right\rangle=\left\langle
0,b\right|H\left|i-j,a\right\rangle.$ (2.3)

Translation invariance implies that $H$ can be partially diagonalized by
Fourier transform (going to momentum space). For the Fourier transform to be
well defined one needs to assume locality, i.e.

$\sum_{j}\Big{|}\left\langle 0,b\right|H\left|j,a\right\rangle\Big{|}<\infty$
(2.4)

Locality, Eq. 2.4, implies that the Bloch Hamiltonian $H(k)$ defined by

$H_{ba}(k)=\sum_{j}\left\langle j,b\right|H\left|0,a\right\rangle e^{-ikj}$
(2.5)

is a continuous periodic matrix-valued function of $k$. $H(k)$ is Hermitian
for each $k$ since

$H_{ba}(k)=\sum_{j}\left\langle j,b\right|H\left|0,a\right\rangle e^{-ikj}\\\
=H_{ab}^{*}(k)$ (2.6)

$H(k)$ is not uniquely defined: There is a $k$-dependent gauge freedom
associated with the choice of unit cell as shown in the next section.

###  2.1 Gauge ambiguity of Bloch Hamiltonians

Changing the unit cell permutes the basis vectors, i.e.

$G^{\ell}\left|j,b\right\rangle=\left|j-b\ell,b\right\rangle,\quad
b\in\\{0,1\\},\quad\ell\in\mathbbm{Z}$ (2.7)

See Fig. 1 for an illustration. Suppose $\ell=1$.

Figure 2: The two red ellipses correspond to two different choices of unit
cell. Changing the unit cell is a gauge transformation (=diagonal unitary)
acting on $H(k)$.

In coordinate space $H$ transforms by

$\left\langle j,b\right|GHG^{\dagger}\left|i,a\right\rangle=\left\langle
j+b,b\right|H\left|i+a,a\right\rangle$ (2.8)

and $H(k)$ transforms by

$\big{(}G_{k}H(k)G_{k}^{\dagger}\big{)}_{b,a}=H_{ba}(k)e^{ik(b-a)}.$ (2.9)

The action of $G_{k}$ in $k$-space is implemented by a diagonal, k-dependent,
unitary $2\times 2$ matrix:

$G_{k}={\begin{pmatrix}1&0\\\ 0&e^{ik}\end{pmatrix}}$ (2.10)

###  2.2 Periodic matrix valued functions and loops

A Hermitian $2\times 2$ matrix can be written as

$H=\mathbf{x}\cdot\boldsymbol{\sigma}+\mathds{1}\,t,\quad\mathbf{x}\in\mathbbm{R}^{3},\
t\in\mathbbm{R}$ (2.11)

where $\boldsymbol{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ is a vector of
Pauli matrices. The matrix is simple (non-degenerate) if

$\|\mathbf{x}\|>0$ (2.12)

The projection on the ground state of $H$ is given by

$P=\frac{\mathds{1}-\mathbf{\hat{x}}\cdot\boldsymbol{\sigma}}{2},\quad\mathbf{\hat{x}}=\frac{\mathbf{x}}{\|\mathbf{x}\|}$
(2.13)

In particular, if $H$ is a Bloch Hamiltonian, the projection on its bottom
band (=valence band) only depends on the traceless part of $H$. We can
therefore restrict ourselves to studying traceless matrices.

A loop of Hermitian gapped $2\times 2$ matrices is in $1-1$ correspondence
with a closed curve $\gamma$ in Euclidean 3-space with the origin removed:

$\gamma=\\{\mathbf{x}_{k}|-\pi\leq k\leq\pi\\}\in\mathbbm{R}^{3}/\\{0\\}.$
(2.14)

Similarly, a loop of rank 1 projections is in $1-1$ correspondence with a
closed curve $\hat{\gamma}$ on the 2-sphere

$\hat{\gamma}=\\{\mathbf{\hat{x}}_{k}|-\pi\leq k\leq\pi\\}\in\mathbbm{S}^{2}.$
(2.15)

Time-reversal, Particle-Hole, and reflection symmetries all relate $H(k)$ with
$H(-k)$ and therefore $\mathbf{x}_{k}$ with $\mathbf{x}_{-k}$. In order to
study deformation of $\gamma$ that respects the symmetry it is useful to think
of the loop $\gamma$ as the concatenation of two curves

$\gamma=\gamma_{-}\circ\gamma_{+}$ (2.16)

where

$\gamma_{-}=\\{\mathbf{x}_{k}|-\pi\leq k\leq
0\\},\quad\gamma_{+}=\\{\mathbf{x}_{k}|0\leq k\leq\pi\\}.$ (2.17)

with $\gamma_{\pm}$ related by the symmetry.

$k$$0$$\pi$ Figure 3: The one dimensional Brillouin zone is a circle
parameterized by the angular variable $k$. The time-reversal, Particle-Hole
and reflection symmetry, relate $k$ with $-k$ and leave invariant the two
points $k=0$ and $k=\pi$.

##  3 Symmetries

This section is a multilingual dictionary that translates several symmetries
of $H$ in coordinate space to symmetries of $H(k)$ and then to symmetries of
the loops $\gamma$. This will allow us to determine the homotopy classes for
each type of symmetry in the next section888Composed symmetries, such as
chiral symmetry, are relegated to Appendix A.1.. (One can alternate between
the two sections in order to understand each symmetry all at once.)

###  3.1 Time reversal

Time-reversal $\boldsymbol{\Theta}$ is an anti-unitary map. This is a
consequence of the presence of $i$ in the Heisenberg equation

$\dot{\mathbf{A}}=i[\mathbf{H},\mathbf{A}]$ (3.1)

where $\mathbf{A,H}$ are operators. Since
$\\{\boldsymbol{\Theta},\frac{d}{dt}\\}=0$ it follows that
$\\{\boldsymbol{\Theta},i\\}=0$ and so is anti-linear.

Now for the problem studied in this article, where the operators act in Fock
space, it is useful to represent time-reversal as an operation on single-
particle states. By definition, time reversal maps particle creation to
particle creation999The operation that maps creation to annihilation is
Particle-Hole transformation., although a creation operator for a given state
may map to a creation operator for some other state:

$\boldsymbol{\Theta}a^{\dagger}(f)\boldsymbol{\Theta}^{\dagger}=a^{\dagger}\big{(}\Theta(f)\big{)}.$
(3.2)

The bold and ordinary $\Theta$’s refer respectively to the transformation of
the creation operator and the transformation of the wave function associated
to it. As $a^{\dagger}(f)$ is linear in $f$ and $a(g)$ is anti-linear in $g$,
it follows from the commutation relation that ${\Theta}(f)$ is anti-unitary

$\boldsymbol{\Theta}\langle
g|f\rangle\boldsymbol{\Theta}^{\dagger}=\\{\boldsymbol{\Theta}a(g)\boldsymbol{\Theta}^{\dagger},\boldsymbol{\Theta}a^{\dagger}(f)\boldsymbol{\Theta}^{\dagger}\\}=\\{a({\Theta}g),a^{\dagger}({\Theta}f)\\}=\langle\Theta
g|\Theta f\rangle$ (3.3)

Since $\boldsymbol{\Theta}$ is antiunitary, this expression is also equal to
$\langle f|g\rangle$ which forces $\Theta$ to be antiunitary. Hence

$\Theta=U\circ*$ (3.4)

with $U$ unitary and $*$ complex conjugation.

One reasonably assumes that reversing time twice produces no observable
effect, i.e., it is a pure phase. Then it follows that, since $\Theta$
commutes with $\Theta^{2}$, the only phase consistent with the anti-unitarity
of $\Theta$ is $\pm 1$, i.e.:

$\Theta^{2}=UU^{*}=\pm\mathds{1}\Longrightarrow U=\pm U^{t}$ (3.5)

When $\Theta^{2}=\pm\mathds{1}$, the symmetry will be denoted
${\Theta}_{\pm}$. Under unitary change of bases $W$

$\Theta\mapsto W\Theta W^{\dagger}=WUW^{t}\circ*$ (3.6)

$U$ then undergoes a congruence transformation.

We assume that $U$ acts on the internal coordinates $\left|a\right\rangle$ in
real space while leaving the spatial coordinate $j$ unchanged. It follows from
$[H,\Theta]=0$ in real space that

$\Theta H(k)\Theta^{-1}=UH^{*}(-k)U^{\dagger}$ (3.7)

Hence, $H(k)$ is time-reversal invariant if

$H(k)=UH^{*}(-k)U^{\dagger}$ (3.8)

for an appropriate unitary $U=\pm U^{t}$.

####  3.1.1 $\Theta_{+}$

For the symmetry $\Theta_{+}$, $U$ is a symmetric matrix. Any symmetric
unitary $2\times 2$ matrix $U$ is congruent to the identity101010We thank Oded
Kenneth for the proof below.. This follows from the fact that $U$ can be
written as

$U(\theta,\phi,\alpha)=e^{i\alpha}[\mathds{1}\cos\theta+i\sin\theta({\sigma_{x}}\cos\phi+{\sigma_{z}}\sin\phi)]$
(3.9)

$U$ has a symmetric square root $\sqrt{U}=U(\theta/2,\phi,\alpha/2)$. It
follows that $U$ is congruent to the identity:

$U=\sqrt{U}\mathds{1}(\sqrt{U})^{t},\quad$ (3.10)

Pick $U=\mathds{1}$, so

$\Theta_{+}=*$ (3.11)

Then

$H(k)=H^{*}(-k)$ (3.12)

translated to a loop $\gamma$ with

$z_{k}=z_{-k},\quad x_{k}=x_{-k},\quad y_{k}=-y_{-k}$ (3.13)

$\gamma$ is symmetric under reflection in $y$ and is anchored at $k=0,\pi$ to
the punctured $x-z$ plane, see Fig. 4.

####  3.1.2 $\Theta_{-}$

There is a unique $2\times 2$ anti-symmetric unitary matrix $U$ up to
multiplication by a scalar. Since an overall phase is for free in an anti-
unitary transformation, we may chose $U={\sigma_{y}}$

$\Theta={\sigma_{y}}\circ*,\quad{\sigma_{y}}=\begin{pmatrix}0&-i\\\
i&0\end{pmatrix}$ (3.14)

whose action on $H(k)$ is

$\Theta H(k)\Theta^{-1}={\sigma_{y}}H^{*}(-k){\sigma_{y}}$ (3.15)

A $2\times 2$ matrix $H(k)$ is (Fermionic) time reversal symmetric if

$\mathbf{x}_{k}=-\mathbf{x}_{-k}$ (3.16)

Since

$\mathbf{x}_{0}=\mathbf{x}_{\pi}=0$ (3.17)

the Bloch Hamiltonian is gapless and the projection $P$, Eq. 2.13, is ill-
defined at $k=0,\pi$. The family of Fermionic time reversal lies outside the
gapped 2 bands framework and shall not be considered further.

Figure 4: The curve $\gamma_{-}$ associated with time-reversal symmetric
Hamiltonians is anchored on the punctured $x-z$ plane. The rest of the loop is
fixed by symmetry of reflection in the plane.

###  3.2 Particle-Hole transformation

Particle-Hole transformation $\mathbf{C}$ maps particle creation to particle
annihilation:

$\mathbf{C}a^{\dagger}(f)\mathbf{C}^{\dagger}=a\big{(}Cf\big{)}$ (3.18)

The particle creation operator $a^{\dagger}(f)$ depends linearly on $f$ while
the annihilation operator is anti-linear. $\mathbf{C}$ (the many-body
operator) must commute with $i$ in order to preserve the direction of time.
Thus, in order for both sides to be linearly dependent on $f$, $C$ (the one-
particle operator) must be anti-linear. Being anti-unitary, $C$ comes in two
varieties

$C^{2}=\pm\mathds{1}$ (3.19)

just like time-reversal. ${C}_{\pm}$ denotes the symmetry corresponding to
$C^{2}=\pm\mathds{1}$.

####  3.2.1 Particle-Hole symmetry

Particle-Hole symmetry says that for any eigenstate of the single particle
Hamiltonian $H$ there is a mirror state with an opposite energy and momentum.
$C$ maps from the single particle state to its mirror image. For Bloch
Hamiltonians this means

${(H)(k)+(CHC^{-1})(-k)=0}$ (3.20)

In contrast, the many-particle Hamiltonian $\mathbf{H}$ commutes with
$\mathbf{C}$. To see this write

$\mathbf{H}=\sum\epsilon_{n,k}^{\phantom{\dagger}}a^{\dagger}_{n,k}a_{n,k}^{\phantom{\dagger}}$
(3.21)

and then, formally denoting the mirror of band $n$ by $-n$,

$\displaystyle\mathbf{C}\mathbf{HC}^{\dagger}$
$\displaystyle=\sum\epsilon_{n,k}^{\phantom{\dagger}}a_{-n,-k}^{\phantom{\dagger}}a_{-n,-k}^{\dagger}=\sum{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\epsilon}_{-n,-k}^{\phantom{\dagger}}a_{n,k}^{\phantom{\dagger}}a_{n,k}^{\dagger}$
$\displaystyle=-\sum\epsilon_{n,k}^{\phantom{\dagger}}\big{(}\\{a_{n,k}^{\phantom{\dagger}},a_{n,k}^{\dagger}\\}-a^{\dagger}_{n,k}a_{n,k}^{\phantom{\dagger}}\big{)}=\mathbf{H}-\underbrace{\sum\epsilon_{n,k}}_{=0}$
$\displaystyle=\mathbf{H}$ (3.22)

It follows from Eq. 2.6, the anti-unitarity of $C$ and Eq. 3.20 that a Bloch
Hamiltonian has Particle-Hole symmetry if

$H(k)=-(CHC^{-1})({-}k)=-UH^{*}(-k)U^{\dagger};\quad C=U\circ*.$ (3.23)

####  3.2.2 $\mathbf{C}_{+}=1$

By section 3.1.1 any anti-unitary transformation whose square is $\mathds{1}$
is equivalent to any another. Pick111111Had we picked the same $U$ for time-
reversal and Particle-Hole, the only Hamiltonian which is symmetric under both
would have been $H=0$. $U={\sigma_{z}}$, i.e.

$C=\sigma_{z}\circ*$ (3.24)

The symmetry of $\gamma$ corresponding to

$H(k)=-\sigma_{z}H^{*}(-k)\sigma_{z}$ (3.25)

is

$z_{k}=-z_{-k},\quad x_{k}=x_{-k},\quad y_{k}=-y_{-k}$ (3.26)

$\gamma$ is symmetric under rotation by $\pi$ around the $x-$axis and is
anchored at $k=0,\pi$ to the punctured $x-$axis, see Fig. 5.

${x}$ Figure 5: The curve $\gamma_{-}$ associated with Particle-Hole symmetric
Hamiltonians with $C^{2}=\mathds{1}$ is anchored on the punctured $x$-axis.
$\gamma_{+}$ is fixed by the symmetry as the $\pi$ rotation of $\gamma_{-}$
around the $x$-axis.

####  3.2.3 $\mathbf{C}_{-}$

By section 3.1.2 we have:

$C={\sigma_{y}}\circ*$ (3.27)

whose action on $H(k)$ is

$CH(k)C^{-1}={\sigma_{y}}H^{*}(-k){\sigma_{y}}.$ (3.28)

So a $2\times 2$ matrix $H(k)$ is Particle-Hole symmetric if

$H(k)=-{\sigma_{y}}H^{*}(-k){\sigma_{y}}=H(-k)$ (3.29)

The loop $\gamma$ corresponding to Eq. 3.29 is self-retracing

$\mathbf{x}_{k}=\mathbf{x}_{-k}$ (3.30)

$\mathbf{x}_{0}$ and $\mathbf{x}_{\pi}$ can then be anywhere in
$\mathbbm{R}^{3}/\\{0\\}$.

###  3.3 Bond reflection symmetry

Figure 6: Bond reflection symmetry is a reflection about the green dashed
line.

Bond reflection (see Fig. 6) permutes the basis vectors

$B\left|j,a\right\rangle=\left|-j,a\oplus 1\right\rangle,\quad B^{\dagger}=B$
(3.31)

where $\oplus$ is addition modulo 2. The transformation reflects unit cells,
and exchanges the “atoms” in a unit cell. $B$ is an involution,
$B^{2}=\mathds{1}$.

The action of bond-reflection on the Hamiltonian in the coordinate basis is

$\left\langle
j,b\right|BHB^{\dagger}\left|i,a\right\rangle=\left\langle-j,b\oplus
1\right|H\left|-i,a\oplus 1\right\rangle$ (3.32)

It follows from Eq. 2.5 that the action in $k$-space is

$(BHB^{\dagger})_{b,a}(k)=H_{b\oplus 1,a\oplus 1}(-k)$ (3.33)

In matrix form:

$(BHB^{\dagger})(k)=\begin{pmatrix}-z_{-k}&x_{-k}{+}iy_{-k}\\\
x_{-k}{-}iy_{-k}&z_{-k}\end{pmatrix}={\sigma_{x}}H(-k){\sigma_{x}}$ (3.34)

$H(k)$ is bond-reflection symmetric if

$H(k)={\sigma_{x}}H(-k){\sigma_{x}}$ (3.35)

The loop $\gamma$ corresponding to Eq. 3.35 is symmetric under a rotation by
$\pi$ about the $x$-axis:

$z_{k}=-z_{-k},\quad x_{k}=x_{-k},\quad y_{k}=-y_{-k}$ (3.36)

and is anchored to the punctured $x-$axis, see Fig. 5.

###  3.4 Site reflection symmetry

Figure 7: Site reflection symmetry is a reflection about the green dashed line

Reflection about a site (see Fig. 7) permutes the basis vectors by

$S\left|j,a\right\rangle=\left|-j-a,a\right\rangle,\quad S=S^{\dagger}$ (3.37)

$S$ is an involution that mixes different unit cells. Its action on $H$ in
real space is given by

$\left\langle j,b\right|SHS^{\dagger}\left|i,a\right\rangle=\left\langle-
j-b,b\right|H\left|-i-a,a\right\rangle$ (3.38)

By Eq. 2.5 site reflection acts in $k$-space by

$(SHS^{\dagger})_{b,a}(k)=H_{b,a}(-k)e^{ik(b-a)}$ (3.39)

and in matrix form

$(SHS^{\dagger})(k)=G_{k}H(-k)G^{\dagger}_{k},$ (3.40)

where $G_{k}$ is given in Eq. 2.10. Hence, $H(k)$ is site-reflection symmetric
if

$H(k)=G_{k}H(-k)G^{\dagger}_{k}$ (3.41)

The symmetry of the loop corresponding to Eq. 3.41 is a relation between
$\gamma_{-}$ and $\gamma_{+}$:

$z_{k}=z_{-k},\quad\begin{pmatrix}x_{k}\\\
y_{k}\end{pmatrix}=\begin{pmatrix}\cos k&-\sin k\\\ \sin k&\cos
k\end{pmatrix}\begin{pmatrix}x_{-k}\\\ y_{-k}\end{pmatrix}$ (3.42)

$\gamma$ is anchored at one end, $k=\pm\pi$, to the punctured $z$-axis

${x_{\pm\pi}=y_{\pm\pi}=0,\quad|z_{\pm\pi}|>0}$ (3.43)

###### Example 3.1

A gapped Bloch Hamiltonian with site reflection symmetry is the SSH model [27]
with equal hopping and on-site potential:

$H(k)=\begin{pmatrix}v&1+e^{-ik}\\\ 1+e^{ik}&-v\end{pmatrix},\quad v>0$ (3.44)

##  4 Homotopy of loops of gapped $2\times 2$ matrices

Two Bloch Hamiltonians are homotopic if one can be deformed to the other,
within a given symmetry class, without closing the gap. By section 2.2 the
question reduces to the homotopy of continuous closed loops
$\gamma\in\mathbbm{R}^{3}/\\{0\\}$ constrained by the symmetry. The symmetry
constraint can sometimes complicate the question of homotopy equivalence, see
section 4.6. Fortunately, for most of the constraints we consider, this is not
the case.

When the symmetries exchange $k$ and $-k$, the problem is simplified by
writing $\gamma$ as the concatenation $\gamma=\gamma_{-}\circ\gamma_{+}$ with
$\gamma_{\pm}$ given by 2.17. The symmetry says that $\gamma_{-}$ determines
$\gamma_{+}$ (and vice versa), and imposes constraints on the end-
points121212In some cases $\gamma$ is also constrained to be planar. of
$\gamma_{\pm}$. Then any homotopy of the full curve $\gamma$ is determined by
the homotopy of the arc $\gamma_{-}$, with the end-points required to satisfy
the constraint. The homotopy equivalence of curves $\gamma_{-}$ is readily
visualized and the determination of the equivalence classes reduces to an
elementary exercise. Homotopy classifications for chiral symmetry are
relegated to an appendix.

###  4.1 No symmetry

$\mathbbm{R}^{3}/\\{0\\}$ is simply connected. This means that every
(continuous) loop can be contracted to a point. Since
$\mathbbm{R}^{3}/\\{0\\}$ is connected, every point is homotopic to every
other point. Hence any loop of $2\times 2$ gapped Bloch matrices $H(k)$ is
null-homotopic (and similarly for any loop of projections). We identify the
equivalence class with the Bloch matrix $H(k)=\sigma_{z}$ representing a
periodic array of isolated cells. This choice is not canonical. Any $2\times
2$ simple (=non-degenerate) Bloch matrix is homotopic to any other non-
degenerate Bloch matrix and the choice $\\{\sigma_{z}\\}$ is just that, a
choice.

The corresponding many-body ground state is a (formal) pure product state

$\prod_{j}\left|0\right\rangle_{j}$ (4.1)

Being homotopic to a pure product states is taken as the defining property of
a trivial phase. This gives the first line in Table 1.

###  4.2  $\Theta_{+}$

By section 3.1.1 $\gamma$ is symmetric under reflection in $y-$axis, and is
anchored to the $x-z$ plane for $k=0,\pi$. Consider the corresponding
$\hat{\gamma}$ on the unit sphere whose end-points are anchored to the equator
in the $x-z$ plane, see Fig. 4.

$\hat{\gamma}_{-}$ can be deformed provided it respects the anchoring to the
equator. $\hat{\gamma}_{+}$ follows by the symmetry constraint. Suppose first
that the anchoring points coincide. Since $\mathbbm{S}^{2}$ is simply
connected, $\hat{\gamma}_{-}$ can be contracted to a point. The mirror image
$\hat{\gamma}_{+}$ will contract simultaneously.

In the case that the anchoring points do not coincide the end-point at $k=\pi$
can be brought to the end-point at $k=0$ by applying a rotation around the
$y-$axis $R_{y}(\phi)$ (since both lie on the unit circle in the $x-z$ plane).
The angle $\phi$ is imagined to increase continuously from 0 as a function of
$k$ so that the points coincide. The $k$ dependent rotation

$R_{y}(k\phi/\pi)\hat{\mathbf{x}}_{k},\quad 0\leq k\leq\pi$ (4.2)

is a homotopy of $\hat{\gamma}_{-}$ which respects the gap condition, since a
rotation preserves the length of vectors. Hence $\hat{\gamma}_{-}$ is
homotopic to a point. Finally, since the equator in the $x-z$ plane is
connected, all points on the equator are homotopic.

It follows that the space of time reversal 1D Bloch Hamiltonians with a gap
condition has a single component under homotopy, which we label by a single
Pauli matrix $\\{\sigma_{z}\\}$. This accounts for the second line in table 1.

###  4.3 Particle-Hole symmetry

####  4.3.1 $C_{+}$

By section 3.2.2 the closed loop $\gamma$ is symmetric under rotation by $\pi$
about the $x-$axis, and is anchored to the (punctured) $x-$axis for $k=0,\pi$.
$\gamma$ is homotopic to $\hat{\gamma}$ on the unit sphere, anchored at the
poles:

$\hat{\mathbf{x}}_{0},\hat{\mathbf{x}}_{\pi}=\pm(1,0,0)$ (4.3)

When the anchoring points coincide $\hat{\gamma}_{-}$ is homotopic to a point.
In the case that the anchoring points are antipodal, $\hat{\gamma}_{-}$ can be
continuously shortened to a geodesic. This means that $\hat{\gamma}$ is
homotopic to a circle.

${\sigma_{x}}$ Figure 8: $\gamma_{-}$ with Particle-Hole symmetry. The blue
spherical curve is anchored at $x=\pm 1$. The red one is anchored at $x=-1$.

In conclusion, Particle-Hole symmetric Hamiltonians have four distinct
equivalence classes represented by the four matrices:

$\\{\pm{\sigma_{x}},\pm R_{1}\\},{\mathrm{where\
}}R_{n}=\begin{pmatrix}0&e^{{-}ink}\\\ e^{ink}&0\end{pmatrix}$ (4.4)

($R\cong R^{\dagger}$ by a $\pi$ rotation about the $x$ axis.) The four
elements are distinguished by the combinations of signs of the anchoring
points on $x$.

$\pm{\sigma_{x}}$ is a Bloch Hamiltonian with disconnected cells, so the
ground states are pure product states

$\prod_{j}\left|+\right\rangle_{j},\quad\prod_{j}\left|-\right\rangle_{j},\qquad\left|\pm\right\rangle=\frac{\left|0\right\rangle\pm\left|1\right\rangle}{\sqrt{2}}$
(4.5)

representing trivial phases, while the Hamiltonians $\pm R_{1}$ when solved on
a finite chain turn out to have edge states at the ends.

####  4.3.2 $C_{-}$

By section 3.2.3 the loop $\gamma$ is self-retracing and is not anchored. Thus
it can be contracted to a point along itself, so it is null-homotopic. This
gives the fifth line in table 1.

###  4.4 Bond reflection symmetry

By section 3.3 the curve $\gamma$ is anchored on the punctured $x$-axis. This
is the same scenario as in section 4.3.1 and hence there are four equivalence
classes labeled by the four matrices

$\\{\pm{\sigma_{x}},\ \pm R_{1}\\}$ (4.6)

This gives the sixth line in table 1.

###  4.5 Site reflection symmetry

By section 3.4 the loop $\gamma$ of Eq. 3.42 is anchored at one point,
$k=\pm\pi$, to the punctured $z$-axis:

$x_{\pm\pi}=y_{\pm\pi}=0,\quad|z_{\pm\pi}|>0$ (4.7)

Since $\gamma_{-}$ is anchored at the single point $k=-\pi$ it can be
contracted to the anchoring point by the homotopy

$\mathbf{x}_{k}(t)=\mathbf{x}_{(1-t)k-t\pi},\quad t\in[0,1]$ (4.8)

It follows that $\hat{\gamma}$ is homotopic to either the north or south pole,
and hence site symmetric Bloch Hamiltonians are homotopic to either one of two
matrices

$\\{\pm{\sigma_{z}}\\}$ (4.9)

The corresponding (many-body) ground state represents a trivial phase as it is
(formally) a pure product state. This gives the seventh line in table 1.

Figure 9: Curves anchored at the north pole are homotopic to the north pole.
And curves anchored at the south pole are homotopic to the south pole. This
shows that the space of loops is made of two disconnected sets.

###  4.6 $S\circ\Theta_{+}$: The world-line of the loop.

It is sometimes useful to consider the loop $\gamma$ as a world-line in the
four dimensional space $\mathbbm{S}^{1}\times\mathbbm{R}^{3}$ where the first
coordinate is $k$. A case in point is the composed symmetry
$S\circ\Theta_{+}$, where the Bloch Hamiltonian satisfies the constraint

$H(k)=G_{k}H^{*}(k)G_{k}^{\dagger},\quad G_{k}={\begin{pmatrix}1&0\\\
0&e^{ik}\end{pmatrix}}$ (4.10)

This constrains the coordinates $(x_{k},y_{k})$:

$x_{k}-iy_{k}=(x_{k}+iy_{k})e^{-ik}$ (4.11)

andd gives no constraint on $z_{k}$. The constraint implies

$(x_{k},y_{k})=\lambda_{k}\left(\cos\frac{k}{2},\sin\frac{k}{2}\right),$
(4.12)

where $\lambda_{k}$ is a real-valued function of $k$.

The loop $\gamma$ may be viewed as the world line in the 3-dimensional
manifold $\mathbbm{M}_{S\circ\Theta_{+}}$

$\mathbbm{M}_{S\circ\Theta_{+}}=\left\\{\left(k,\lambda\cos\frac{k}{2},\lambda\sin\frac{k}{2},z\right)\Big{|}\lambda,z\in\mathbbm{R},-\pi\leq
k\leq\pi\right\\}$ (4.13)

embedded in the 4-dimensional space $\mathbbm{S}^{1}\times\mathbbm{R}^{3}$.
$\gamma$ is not a closed curve in the $(\lambda,z)$ coordinates because the
$\lambda-z$ plane rotates relative to the $x-y$ plane by $\pi$ as $k$ goes
from $-\pi$ to $\pi$. The end-point $\mathbf{x}_{-\pi}=\mathbf{x}_{\pi}$ has
different coordinates in the $(\lambda,z)$ plane

$\lambda_{-\pi}=-\lambda_{\pi},\quad z_{-\pi}=z_{\pi}$ (4.14)

The homotopy classes of Hamiltonians with $\mathrm{S}\circ\Theta$ symmetry are
thus homotopy classes of curves in a punctured $\lambda-z$ plane whose end-
points are constrained to be mirror-images in the $z$-axis, see Fig. 10.

${\sigma_{z}}$$-\sigma_{z}$$\lambda$$z$ Figure 10: $\gamma$ in the $\lambda-z$
plane. The start and end points $(\lambda_{\pm\pi},z_{\pm\pi})$ are the two
cyan circles that are mirror imagees. The red circle in the puncture in the
plane. The curve can be retracted to ${\sigma_{z}}$ by retracting the end
points symmetrically relative to the $z$ axis. The curve can not be retracted
to $-{\sigma_{z}}$.

The end-points can be deformed to be on the $+z$-axis (since any point and its
mirror image can be moved continuously to $+z$). This curve is closed so has a
winding number. If one of the end-points is taken around the origin $n$ times
and brought back (while the other end-point mirrors the first), the winding
number changes by $2n$. Thus, there are only two classes, characterized by
even and odd winding number around the origin (after the end-points are moved
to the $+z$-axis). The curves with even winding number can be contracted to
constants $+\sigma_{z}$ and the curves with odd winding number can be
contracted to $-\sigma_{z}$ (by moving the two end-points half-way around a
circle in opposite directions) so all Hamiltonians are homotopic to either
$\pm\sigma_{z}$131313To distinguish between these two cases without first
moving the end-points together, count the number of times the interior of the
curve crosses the $-z$ axis. The parity is a homotopy invariant, and is
equivalent to the previous classification..

##  5 Changing unit cell permutes homotopy classes

The labeling of the homotopy classes depends, in general, on the choice of
unit cell [2, 28]. Different choices of unit cell are related by a gauge
transformation, section 2.1,

$H(k)\mapsto G_{k}^{\ell}H(k)G_{k}^{-\ell},\quad\ell\in\mathbbm{Z}$ (5.1)

Under such a transformation the equivalence classes transform as

${\sigma_{z}}\mapsto{\sigma_{z}},\quad{\sigma_{x}}\mapsto R_{\ell},\quad
R_{n}\mapsto R_{n+\ell}.$ (5.2)

This does not mean that the phases are equivalent, just that the labeling of
the phase depends on how the unit cell is chosen.

For example, in the case of bond-reflection symmetry, the classes
$\\{\pm\sigma_{x}\\}$ are exchanged with $\\{\pm R_{1}\\}$ under a gauge
transformation141414The topological class of a perfectly periodic system has
consequences for the behavior when the periodicity is broken by ending the
chain, as long as far away from the edge the system is described by the
periodic Hamiltonian. If the chain is stopped after a unit cell (whichever
unit cell has been used to express the Bloch Hamiltonian), there will be an
edge state for the $R_{1}$ e̊quivalence class but not for the $\sigma_{x}$
class.. This is the case in the Su-Schrieffer-Heeger model [27].

##  6 Fragile and New phases

The classification of Bloch Hamiltonians with two bands in 1D can be extended
to Bloch Hamiltonians with more bands. Fixing the number of bands to be $n$
and the number of occupied bands to be $m<n$, we are interested in the
homotopy of loops of $n\times n$ invertible Hermitian matrices with $m$
negative and $n-m$ positive eigenvalues.

The homotopic classification in the general case may differ from the
classification of two bands: First, distinct phases in the two bands case may
dissolve into one equivalence class upon the addition of bands. Second, new
phases may appear.

###  6.1 A fragile phase

A phase (with a given symmetry) is fragile [29] if its homotopy class changes
when the dimension of the matrices of band projection $P(k)$ is increased:

$P(k)\mapsto\begin{pmatrix}P(k)&0\\\ 0&0\end{pmatrix}$ (6.1)

As an example, of a fragile topological phase consider an $n\times n$ real,
periodic, projection $P(k)$ onto a single band:

$P(k)=P^{*}(k)=P(k+2\pi)=P^{2}(k),\quad Tr\,P(k)=1$ (6.2)

The reality condition may be viewed to be a consequence of the symmetry
constraint $B\circ\Theta_{+}$ upon choosing $\Theta_{+}=\sigma_{x}\circ*$.

Since $P(k)$ is a real symmetric matrix it can be diagonalized by an
orthogonal transformation and its non-trivial eigenvector

$P(k)\left|k\right\rangle=\left|k\right\rangle$ (6.3)

can be chosen to be real vectors of unit length in $\mathbbm{R}^{n}$. The
eigenvalue equation then determines $\left|k\right\rangle$ uniquely, up to a
$\pm 1$ phase ambiguity. Since $P(k)$ is continuous in $k$, we may choose
$\left|k\right\rangle$ so that it is a continuous family for $k\in[-\pi,\pi]$.
Then $\left|k\right\rangle$ traces a continuous curve on the unit sphere in
$\mathbbm{R}^{n}$. Since $P(-\pi)=P(\pi)$ and $\left|k\right\rangle$ has $\pm
1$ phase freedom, the curve is either a closed loop or one with antipodal end-
points:

$\left|k=-\pi\right\rangle=\pm\left|k=\pi\right\rangle$ (6.4)

The $\pm$ alternative divides the space of curves $\left|k\right\rangle$ on
the sphere into two classes: A symmetric class where
$\left|k=\pi\right\rangle=\left|k=-\pi\right\rangle$ and an anti-symmetric
class where $\left|k=\pi\right\rangle=-\left|k=-\pi\right\rangle$. This holds
for any $n\geq 2$.

In the case $n=2$ the curves can be classified further according to the
winding number of $\left|k\right\rangle$ on the circle, which is a whole
number (for the symmetric class) or a half-integer for the anti-symmetric
class, see Fig. 11. This can be related to the winding number of the curve
$\gamma$ representing the Hamiltonian: Since $H$ is real, $\gamma$ is in the
$x-z$ plane. $\gamma$ winds twice as fast as the eigenvector, so its winding
number is an integer (since the Hamiltonian is periodic). It is odd for an
anti-symmetric eigenvector and even for a symmetric one.

For $n\geq 3$ the vector $\left|k\right\rangle$ lies on $\mathbbm{S}^{n-1}$
which is simply connected, so all integer winding can be undone, reducing to a
single class for the symmetric case and a single class for the anti-symmetric
case.

Thus, we find that the periodic real projections $P(k)$ are labeled by:

$\mathbbm{Z}\text{ for }n=2,\quad\mathbbm{Z}_{2}\text{ for }n\geq 3$ (6.5)

Hence all real Bloch bands with an even winding number of $\gamma$ collapse to
a single equivalence class, while those with an odd winding number collapse to
a second equivalence class, once the dimension of the Hamiltonian is increased
above 2.

Figure 11: The blue solid semi-circle and the dashed semi-circle are
homotopically distinct when the embedding space is the plane but are
equivalent in three dimensions.

###  6.2 New phases

As an example of new phases that emerge when the number of bands is increased
consider the reflection symmetry ${R}$. The points $k=0,\pi$ are special as

$[{R},H(k)]=0,\quad k\in\\{0,\pi\\}$ (6.6)

Hence, eigenstates of $H(k),$ for $k=0,\pi$ are also eigenvectors of ${R}$.
Define the index $N_{-}(k)\in\\{0,\dots,m\\}$ as the sum of negative states
with negative parity:

$N_{-}(k)=\sum_{\epsilon_{n}<0}\frac{1-\left\langle\psi_{n}\right|{R}\left|\psi_{n}\right\rangle}{2},\quad
H(k)\left|\psi_{n}\right\rangle=\epsilon_{n}\left|\psi_{n}\right\rangle,\quad
k\in\\{0,\pi\\}$ (6.7)

$N_{-}$ is an invariant under deformations of $H(k)$ that respect the parity
and do not close the gap151515The number of negative energy states with
eigenvalues $+1$ is not a separate invariant since the total number of
negative eigenvalues is fixed, see [30].

In the case of $m$ occupied bands one expects to find, in general, $(m+1)^{2}$
equivalence classes labeled by $(N_{-}(0),N_{-}(\pi))$ and the number of
invariants to increase with the number of the bands.

It is instructive to examine the results in Table 1 where $m=1$, from this
perspective:

  * •

In the case of bond reflection $R_{k}=\sigma_{x}$ and $N_{-}(k)\in\\{0,1\\}$.
Table 1 indeed gives four equivalence classes, as it should.

  * •

In the case of site reflection $R_{k}=G_{k}$ is the identity for $k=0$ so
automatically $N_{-}(0)=0$ and the two equivalence classes are labelled by
$N_{-}(\pi)\in\\{0,1\\}$ in agreement with Table 1.

###  6.3 Stable phases

Since the number of invariants may either shrink or grow when the number of
bands increases, it is natural to introduce a notion of “stable invariants”.
For sufficiently large (invertible) matrices where both the numbers of
positive and negative eigenvalues are large, the topological classification
often stops changing, and the “stable invariants” describe these classes. We
shall not discuss further this issue.

##  7 Concluding remarks

Our aim here was to give an elementary invitation to the field of topological
insulators by studying the homotopic classification of $2\times 2$ simple
(=non-degenerate) periodic matrices in one dimension with various symmetries.
Using simple intuitive facts about the homotopy of curves, we have shown how
various symmetries are expressed in the homotopic classification of Bloch
Hamiltonians. We have also presented an introduction to the notions of fragile
and stable topological phases.

## Acknowledgement

JA thanks Johannes Kellendonk, Oded Kenneth, Roy Meshulam, Daniel Podolsky,
Immanuel Bloch and Raffaele Resta for useful discussions.

##  Appendix A Various symmetry compositions

###  A.1 Chiral symmetry: $C_{+}\circ\Theta_{+}$

Chiral symmetry (aka sublattice symmetry) is the _composition_ of Particle-
Hole and time reversal symmetry. As such it is a unitary map that preserves
$k$,

$H(k)\mapsto-UH(k)U^{\dagger},\quad U=U_{c}U_{t}^{*}$ (A.1)

where $\Theta_{+}=U_{t}\circ*$ and $C_{+}=U_{c}\circ*$.  Note that, in order
for chiral symmetry to lead to an interesting classification, one must choose
$U_{t}\neq U_{c}$. If $U_{t}=U_{c}$, this would have the unfortunate
consequence that $U=U_{c}U^{*}_{c}=U_{c}U^{\dagger}_{c}=\mathds{1}$ and chiral
symmetry would imply $H(k)=0$.

Choosing $\Theta_{+}$ and $C_{+}$ as in Eqs. 3.11,3.24, $H(k)$ has chiral
symmetry if

$H(k)=-\sigma_{z}H(k)\sigma_{z}$ (A.2)

This says that the curve $\gamma$ lies in a punctured $x-y$ plane161616Had we
chosen $U=\sigma_{y}$ the curve $\gamma$ would lie in the $x-z$ plane and
$H(k)$ would then be a real matrix., which is not simply connected and has its
winding number as a topological invariant

$\text{Winding}(\gamma)\in\mathbbm{Z}$ (A.3)

$H(k)$ with winding $n$ is homotopic to $R_{n}$ of Eq. 4.4. It follows that
chiral Bloch Hamiltonians are homotopic to

$\\{R_{n}|n\in\mathbbm{Z}\\}$ (A.4)

###  A.2
$C_{+}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\land}\Theta_{+}$

Particle-Hole and chiral symmetry together also imply symmetry under their
composition, chiral symmetry. Chiral symmetry restricts $\gamma$ to the $x-y$
plane and time reversal forces the anchoring points to lie in the $x-z$ plane.
It follows that the anchoring points at $0,\pi$ must lie on a punctured
$x$-axis.  Although $\gamma$ and $-\gamma$ have the same winding numbers and
can be rotated into one another by applying a $\pi$ rotation about the
$z$-axis, they are homotopically distinct since such a rotation does not
respect the anchoring condition. This doubles the homotopy classes

$R_{n}\mapsto\pm R_{n}.$ (A.5)

To summarize, the homotopy equivalence classes are given by

$\\{\pm R_{n}|n\in\mathbbm{Z}\\}$ (A.6)

###  A.3 $B\circ\Theta_{+}$

For Hamiltonians with only one symmetry, the _composition_ of time reversal
and bond reflection, one has

$H(k)=(B\Theta H\Theta^{-1}B)(k)=\sigma_{x}H^{*}(k)\sigma_{x}$ (A.7)

which implies

$z_{k}=0$ (A.8)

so $\gamma$ lies in the punctured $x-y$ plane, and there is no constraint on
the anchoring points at $k=0,\pi$. The equivalence classes are labelled by
their winding number and all Bloch Hamiltonians with winding $n$ are homotopic
to the matrix $R_{n}$. This gives the ninth line in table 1:

$\\{R_{n}|n\in\mathbbm{Z}\\}.$ (A.9)

###  A.4 $B\land\Theta_{+}$

Now, in addition to Eqs. A.7,A.8 bond reflection gives the additional
constraint that $\gamma$ is anchored on the punctured $x-$axis. It follows
that $\hat{\gamma}$ is distinguished by its winding, and the anchoring point
$\mathbf{x}_{0}=(\pm 1,0,0)$. The equivalence classes are represented by the
matrices

$\\{\pm R_{n}|n\in\mathbbm{Z}\\}$ (A.10)

This gives the twelfth line in the table 1.

###  A.5 $S\land\Theta_{+}$

When $S$ and $\Theta_{+}$ symmetries are combined, then (as in section 4.6)
the curve $\gamma$ is not in a fixed space, independent of $k$. This is not
anymore the case for Bloch Hamiltonians with time reversal and site reflection
symmetry.

The conditions are that

$(x_{k},y_{k})=\lambda_{k}\left(\cos\frac{k}{2},\sin\frac{k}{2}\right),\quad(\lambda_{k},z_{k})=(\lambda_{-k},z_{-k}),\quad{\lambda_{\pm\pi}=0}$
(A.11)

The first condition follows from the site-reflection condition Eq. 3.42, as in
section 4.6. Time-reversal implies that $x_{k}=x_{-k}$ and $y_{k}=-y_{-k}$;
i.e., $\lambda_{k}=\lambda_{-k}$. Also, $\lambda_{\pm\pi}=0$ by Eq. 4.7.
$\gamma_{-}$ is anchored at $-\pi$ on the punctured $z$-axis. The arc can be
contracted to the point it is anchored to on the north or to the south pole
while respecting the constraint, by contracting it along itself in the
$z-\lambda$ space. I.e. for $k<0$,

$z_{k}(t)=z_{(1-t)k-t\pi},\quad(x_{k},y_{k})=\lambda_{(1-t)k-t\pi}\left(\cos\frac{k}{2},\sin\frac{k}{2}\right)$
(A.12)

There are two homotopy classes determined by the anchoring points:

$\\{\pm{\sigma_{z}}\\}$ (A.13)

This gives the fourth line from the bottom of table 1.

##  Appendix B Comparison with the standard classification

The standard classification of topological insulators is concerned with the 10
types of symmetry that do not have spatial symmetries other than translations.
It is made from all the possible combinations of time reversal and Particle-
Hole symmetry. The standard classification and the two-band classification are
compared in the table below.

  * •

The standard classification labels the homotopy equivalence classes by Abelian
groups where $0$ says that there is one equivalence class. Since Bloch
Hamiltonians do not come with a natural group structure, we found it more
intuitive to mark the equivalence classes by a particularly simple Bloch
matrix in the set. This choice, is, of course, arbitrary. For example,
$\\{\sigma_{x}\\},\\{\sigma_{y}\\}$ and $\\{\sigma_{z}\\}$ in the third
column, all represent a single equivalence class. The number of elements in
the third column gives the number of homotopically distinct sets.

  * •

In the two cases marked by $*$ in first column the two band classification
gives twice as many equivalence classes as in the standard classification. The
extra $\mathbbm{Z}_{2}$ can be traced to the sign of the Pfaffian of $-iH$
where $H$ is Particle-Hole symmetric matrix (and so imaginary and anti-
symmetric). The standard classification disregards this “weak”
$\mathbbm{Z}_{2}$ index because it is a zero-dimensional property.

  * •

In three cases, marked by $-$ in the second column, no comparison is made as
there are no gapped 2-band Hamiltonian with $\Theta_{-}$ symmetry.

For these symmetries, there are no fragile phases, which would have also
caused there to be more phases in our classification than in the standard
classification. Fragile phases do occur when there is reflection symmetry, see
Sec. 6.1.

| Symmetry | Two-Band Classification | Standard Classification
---|---|---|---
A | None | $\\{\sigma_{z}\\}$ | 0
AIII | Chiral | $\\{R_{n},n\in\mathbbm{Z}\\}$ | $\mathbbm{Z}$
AI | $\Theta_{+}$ | $\\{\sigma_{z}\\}$ | 0
BDI∗ |  $\Theta_{+},$C+ | $\\{\pm R_{n}\\}$ | $\mathbbm{Z}$
D∗ | C+ | $\\{\pm\sigma_{x},\pm R_{1}\\}$ | $\mathbbm{Z}_{2}$
DIII† |  $\Theta_{-}$,C+ | – | $\mathbbm{Z}_{2}$
AII† | $\Theta_{-}$ | – | 0
CII† |  $\Theta_{-}$,C- | – | $\mathbbm{Z}$
C | C- | {$\sigma_{z}$} | 0
CI |  $\Theta_{+}$,C- | $\\{\sigma_{y}\\}$ | 0
Table 2: Comparison between standard classification of topological insulator
classification in one dimension and classification for $2\times 2$
Hamiltonians. The classifications agree except that in the standard
classification a weak $\mathbbm{Z}_{2}$ index is ignored in two cases (marked
∗), corresponding to the choice of $\pm$ in the Hamiltonians, and whenever
there is $\Theta_{-}$ symmetry (marked †), there are no $2\times 2$
Hamiltonians at all that are gapped.

## References

  * [1] F. Duncan M. Haldane.  Nobel lecture: Topological quantum matter.  Reviews of Modern Physics, 89(4):040502, 2017.
  * [2] Liang Fu, Charles L. Kane, and Eugene J. Mele.  Topological insulators in three dimensions.  Physical Review Letters, 98(10):106803, 2007.
  * [3] Alexei Kitaev.  Periodic table for topological insulators and superconductors.  In AIP Conference Proceedings, volume 1134, pages 22–30. American Institute of Physics, 2009.
  * [4] Zheng-Cheng Gu and Xiao-Gang Wen.  Tensor-entanglement-filtering renormalization approach and symmetry-protected topological order.  Physical Review B, 80(15):155131, 2009.
  * [5] Alexander Altland and Martin R. Zirnbauer.  Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures.  Physical Review B, 55(2):1142, 1997.
  * [6] Ching-Kai Chiu, Jeffrey C.Y. Teo, Andreas P. Schnyder, and Shinsei Ryu.  Classification of topological quantum matter with symmetries.  Reviews of Modern Physics, 88(3):035005, 2016.
  * [7] Andreas P. Schnyder, Shinsei Ryu, Akira Furusaki, and Andreas W.W. Ludwig.  Classification of topological insulators and superconductors in three spatial dimensions.  Physical Review B, 78(19):195125, 2008.
  * [8] Andreas P Schnyder, Shinsei Ryu, Akira Furusaki, and Andreas WW Ludwig.  Classification of topological insulators and superconductors.  In AIP Conference Proceedings, volume 1134, pages 10–21. American Institute of Physics, 2009.
  * [9] Erez Berg, Michael Levin, and Ehud Altman.  Quantized pumping and topology of the phase diagram for a system of interacting bosons.  Physical Review Letters, 106(11):110405, 2011.
  * [10] Alexei Kitaev.  Anyons in an exactly solved model and beyond.  Annals of Physics, 321(1):2–111, 2006.
  * [11] Xie Chen, Zheng-Cheng Gu, and Xiao-Gang Wen.  Local unitary transformation, long-range quantum entanglement, wave function renormalization, and topological order.  Physical Review B, 82(15):155138, 2010.
  * [12] Michael Levin and Xiao-Gang Wen.  Detecting topological order in a ground state wave function.  Physical Review Letters, 96(11):110405, 2006.
  * [13] Hal Tasaki.  Physics and Mathematics of Quantum Many-Body Systems.  Springer International Publishing, 2020.
  * [14] Yoshiko Ogata.  A $Z_{2}$-index of symmetry protected topological phases with time reversal symmetry for quantum spin chains.  Communications in Mathematical Physics, 374(2):705–734, 2020.
  * [15] Lukasz Fidkowski and Alexei Kitaev.  Topological phases of Fermions in one dimension.  Physical Review B, 83(7):075103, 2011.
  * [16] Frank Pollmann, Ari M. Turner, Erez Berg, and Masaki Oshikawa.  Entanglement spectrum of a topological phase in one dimension.  Physical Review B, 81(6):064439, 2010.
  * [17] Sven Bachmann, Spyridon Michalakis, Bruno Nachtergaele, and Robert Sims.  Automorphic equivalence within gapped phases of quantum lattice systems.  Communications in Mathematical Physics, 309(3):835–871, 2012.
  * [18] Sven Bachmann, Alex Bols, Wojciech De Roeck, and Martin Fraas.  A many-body index for quantum charge transport.  Communications in Mathematical Physics, 375(2):1249–1272, 2019\.
  * [19] Bertrand I. Halperin.  Quantized hall conductance, current-carrying edge states, and the existence of extended states in a two-dimensional disordered potential.  Physical Review B, 25(4):2185, 1982.
  * [20] Gian Michele Graf and Marcello Porta.  Bulk-edge correspondence for two-dimensional topological insulators.  Communications in Mathematical Physics, 324(3):851–895, 2013.
  * [21] Jason Alicea, Yuval Oreg, Gil Refael, Felix von Oppen, and Matthew P. A. Fisher.  Non-abelian statistics and topological quantum information processing in 1d wire networks.  Nature Physics, 7(5):412–417, 2011.
  * [22] Chetan Nayak, Steven H. Simon, Ady Stern, Michael Freedman, and Sankar Das Sarma.  Non-abelian anyons and topological quantum computation.  Reviews of Modern Physics, 80(3):1083, 2008.
  * [23] Julia Kempe, Alexei Kitaev, and Oded Regev.  The complexity of the local Hamiltonian problem.  Siam Journal on Computing, 35(5):1070–1097, 2006.
  * [24] Emil Prodan and Hermann Schulz-Baldes.  Bulk and boundary invariants for complex topological insulators: From K-theory to physics.  Springer, 2016.
  * [25] Jacob Shapiro.  Topology and Localization: Mathematical Aspects of Electrons in Strongly-Disordered Media.  PhD thesis, ETH Zurich, 2018.
  * [26] Ralph H Fox.  On topologies for function spaces.  Bull. Amer. Math. Soc., 51:429–432, 1945.
  * [27] WP Su, JR Schrieffer, and AJ Heeger.  Solitons in polyacetylene.  Physical review letters, 42(25):1698, 1979.
  * [28] Shinsei Ryu and Yasuhiro Hatsugai.  Entanglement entropy and the berry phase in the solid state.  Physical Review B, 73(24):245115, 2006.
  * [29] Hoi Chun Po, Haruki Watanabe, and Ashvin Vishwanath.  Fragile topology and Wannier obstructions.  Physical Review Letters, 121(12):126402, Sep 2018.
  * [30] Ari M. Turner, Yi Zhang, Roger S. K. Mong, and Ashvin Vishwanath.  Quantized response and topology of magnetic insulators with inversion symmetry.  Physical Review B, 85(16):165120, 2012.


# Locomotion Policy Guided Traversability Learning using Volumetric
Representations of Complex Environments

Jonas Frey, David Hoeller, Shehryar Khattak, Marco Hutter J. Frey, D. Hoeller,
S. Khattak and M. Hutter are with the Department of Mechanical and Process
Engineering, ETH Zürich, 8092 Zürich, Switzerland {jonfrey, skhattak,
dhoeller<EMAIL_ADDRESS>D. Hoeller is also associated with NVIDIA.This
work was supported by the Swiss National Science Foundation (SNSF) through
project 166232, 188596, the National Centre of Competence in Research Robotics
(NCCR Robotics), and the European Union’s Horizon 2020 research and innovation
program under grant agreement No.780883.

###### Abstract

Despite the progress in legged robotic locomotion, autonomous navigation in
unknown environments remains an open problem. Ideally, the navigation system
utilizes the full potential of the robots’ locomotion capabilities while
operating within safety limits under uncertainty. The robot must sense and
analyze the traversability of the surrounding terrain, which depends on the
hardware, locomotion control, and terrain properties. It may contain
information about the risk, energy, or time consumption needed to traverse the
terrain. To avoid hand-crafted traversability cost functions we propose to
collect traversability information about the robot and locomotion policy by
simulating the traversal over randomly generated terrains using a physics
simulator. Thousand of robots are simulated in parallel controlled by the same
locomotion policy used in reality to acquire 57 years of real-world locomotion
experience equivalent. For deployment on the real robot, a sparse
convolutional network is trained to predict the simulated traversability cost,
which is tailored to the deployed locomotion policy, from an entirely
geometric representation of the environment in the form of a 3D voxel-
occupancy map. This representation avoids the need for commonly used elevation
maps, which are error-prone in the presence of overhanging obstacles and
multi-floor or low-ceiling scenarios. The effectiveness of the proposed
traversability prediction network is demonstrated for path planning for the
legged robot ANYmal in various indoor and natural environments.

##  I INTRODUCTION

Figure 1: The legged robot ANYmal is deployed in a forest. A robot-centric
occupancy map is generated using two onboard LiDAR sensors. Our real-time
capable proposed sparse 3D neural network predicts the traversability of the
surrounding terrain in 3D, tailored to the deployed locomotion policy. Blue
indicates risk-free traversable, while red indicates traversable under high-
risk.

Humans and animals have the remarkable ability to estimate where they can go
safely and efficiently. The capability to quickly assess one’s surrounding is
an essential selection criterion for survival in a variety of scenarios e.g.,
when a mouse has to escape from a cat into a narrow hole in the wall. Similar
constraint-aware terrain assessment capabilities are desirable for autonomous
mobile robot deployment in complex unstructured environments. This assessment
of the surrounding terrain can be used to navigate safely from a start to a
goal position, during exploration, or to build a better semantic understanding
of the scene. A specific terrain is associated with a cost, which depends on
the capability of the robotic system, locomotion policy, and the structure of
the deployment environment. The traversability assessment is performed based
on multiple application-dependent objectives, which may include: Desired goal
position, time and risk minimization, collision avoidance, safety margins, or
energy efficiency [22].

To acquire information about its surrounding, a robot is equipped with a suite
of exteroceptive sensors. Typically, geometrical information about the
environment is accumulated in a map representation by fusing depth
measurements from LiDARs or depth cameras. For this, ground robots typically
utilize a 2.5D elevation map representation due to its computational
efficiency, in which each cell within a regularly sampled grid at a fixed
resolution stores the estimated height value [5]. During path planning for
ground robots the elevation map can be used to restrict the search space,
given that the base height of the robot is constrained by the underlying
terrain geometry. This reduces the planning time given that not all poses in
3D space have to be considered for path planning.

However, this 2.5D elevation map representation cannot represent multi-floor
environments adequately [14], as without additional semantic understanding,
integration of depth measurements coming from over-hang obstacles or the
ceiling corrupts the elevation map and leads to sub-optimal path planning.
Transitioning to a more expressive 3D voxel occupancy representation [20] of
the environment resolves the aforementioned problems, but requires new methods
for efficient traversability assessment. To this end, this work proposes a
robot-centric traversability estimation module, which predicts the
traversability of the terrain, while being aware of the robots’ specific
locomotion capabilities. As an input, the volumetric occupancy representation
is provided, which captures geometric information sufficient to avoid positive
obstacles, negative obstacles, and steep slopes. To exploit the sparsity of
the underlying problem, we use a sparse 3D convolutional neural network (CNN)
to estimate the terrain traversability in real-time. The network is trained on
data gathered within a highly-efficient physics simulator, where the success
rate of executing different local motion commands on randomly generated
terrain is evaluated inspired by the work of [2, 8, 30]. We deploy the same
locomotion policy in simulation as on the real robots to achieve successful
sim-to-real transfer. The main contributions of the presented work are:

  * •

Predicting terrain traversability using volumetric environment representations
to allow for motion planning in complex environments with over-hanging
obstacles.

  * •

Application of a sparse 3D CNN trained exclusively on simulated data for
traversability estimation.

  * •

Creating a framework for large-scale traversability experience collection
within a physics simulator.

  * •

Deployment in real complex and natural environments and showcasing path
planning application for ANYmal.

##  II RELATED WORK

Figure 2: Collecting traversability information in simulation (left to right):
Birds-eye-view of randomly generated terrain mesh; Zoomed-in view of
collision-free starting points in yellow; Multiple simulated ANYmal robots
executing motion commands; Simulated traversability estimate color-coded.

In this section, we provide an overview of traversability estimation methods
deployed for ground robots. We categorize the approaches based on their input
representation into visual and geometrical methods. Additionally, we provide
an overview of the most relevant work for deep learning on volumetric voxel
data.

###  II-A Visual Traversability Estimation

Images contain rich geometrical and semantic information about the
environment. Teams within the Darpa LAGR project [12] deployed self-supervised
learning approaches to classify from images the affordance necessary to
traverse with wheeled robots outdoor. To overcome the issue of limited
supervision signal Otsu et al. [21] use a support vector machine trained in a
self-supervised fashion to classify terrain for a mars rover. The supervision
signal is extracted from vibration data experienced by the robot and
associated in hindsight with the previously captured images of the terrain.
Wellhausen et al. [27] extend this concept for legged robots and train a CNN
to estimate terrain properties with a supervision signal acquired by force-
torque sensors in the robot’s feet. In conclusion visual data contains rich
information to estimate terrain friction, terrain roughness, foothold quality,
vegetation, support surfaces, and obstacles. Extraction of this information
under sensor noise, motion blur, bad lighting conditions remains a challenging
problem and collecting labeled high-quality data for supervised learning
methods is expensive. Further, using synthetic camera data is often prohibited
by the sim-to-real gap induced by the mismatch between the simulated camera
data and data encountered by a robot during its mission [6].

###  II-B Geometrical Traversability Estimation

Apart from cameras, most mobile robots are equipped with LiDAR and depth
cameras to acquire geometrical information. This data is less rich in semantic
information but easier to process and simulate. Individual depth measurements
can be used for precise localization and fused into different representations
including pointclouds, 2.5D elevation maps [5, 17] and 3D voxel-based
occupancy maps. The most common frameworks to create the latter are Voxblox
[20], Octomap [11], Fiesta [10] or implementations build around the VDB
representation [19].

Wermerlinger et al. [29] use an elevation map and identify the traversability
for a legged robot based on hand-crafted heuristics such as slopes, terrain
roughness, and steps. Different risk sources are assessed from elevation data
[1, 4]. The individual risks are stored in a Multi-Layer Traversability Map
(MLT), which by superposition of all layers, is used to estimate the total
risk for planning. Wellhausen et al. [28] train a CNN to estimate good
footholds for a legged robot from an elevation map. The training dataset of
good footholds is hand-labeled.

In contrast to previous work, Kruesi et al. [14] use pointclouds to estimate
the static terrain traversability for a wheeled robot to achieve safe motion
planning. Methods such as GBPlanner2 [15] illustrate the benefit of a
volumetric environment representation for exploration and global path
planning. The planner is primarily developed for aerial robots, but when
adapted to legged systems, the ground surface estimation is solved based on
heuristics and traversability not considered for planning. All of the previous
geometrical methods rely on hand-crafted heuristics and human expertise to
evaluate traversability.

In contrast, Chavez-Garcia et al. [2] simulate a wheeled robot executing local
motion commands and gather traversability information in simulation. They
train a CNN to classify if a given elevation map is traversable. Guzzi et al.
[8] extend this concept for legged robots and, additionally to the robot-
centric elevation map provide a relative target pose to estimate the time and
probability of reaching the target. Lastly, Yang et al. [30] follow a similar
procedure to simulate the traversability of a given terrain for a legged
robot. Information about the time, risk, and energy needed for a robot to
execute a local motion command is collected in simulation. A CNN network is
trained using this data and predicts the collected quantities similar to [8].
Our approach follows the idea of [8] to simulate the robot to acquire
locomotion policy-guided traversability information. All previous work rely on
elevation maps, which are error-prone in multi-floor or overhanging obstacles
scenario, to estimate traversability during the mission. We instead use a more
expressive occupancy input representation resolving this failure mode.

###  II-C Voxel-based Deep Learning

Predicting traversability from 3D voxel-based occupancy data requires
computational efficient methods mitigating the underlying cubic complexity by
exploiting sparsity. Wang et al. [26] develop a neural network to operate on
an octree volumetric data structure for shape analysis, restraining
computations on occupied voxel octants. Sparse 3D CNNs restrain computation to
occupied voxels [7]. These networks have been generalized to higher dimensions
by Choy et al. [3]. We adopt the sparse encoder-decoder network of Gwak et al.
[9] for semantic segmentation to predict traversability. The network utilizes
sparse generative convolutions and pruning with sequential up-sampling and
pruning voxels in 3D.

##  III APPROACH

###  III-A Problem Definition:

Given an environment representation and the internal state of the robot we aim
to assess if the robot is capable to reside over a given terrain. The
environment is represented in the form of an occupancy voxel map
$\mathbf{O}(x,y,z)\in\\{0,1\\}$ with a resolution of $d_{voxel}$, where each
voxel is either labeled as free(0) or occupied(1). The state of the robot
$\mathbf{s}$ is composed of its position $\mathbf{p}~{}=~{}(x,y,z)^{T}$,
heading $\alpha$, and locomotion policy $\pi_{\theta}$,

$\mathbf{s}=(\mathbf{p},\alpha,\pi_{\theta})^{T}$ (1)

To store the traversability information associated with a terrain we introduce
the notion of a traversability cost tensor $\mathbf{T}$. The tensor stores the
traversability cost at a spatial location, for a set of robot configurations
and motion commands $a$. We simplify the robot configuration to solely include
the heading parameter $\alpha$ and discard velocities, individual joint
positions and base pitch and roll, given that the traversability assessment is
less dependent on these fast varying variables under nominal operating
conditions of the studied legged robot ANYmal. Including the robot
configuration and motion command is motivated by the fact that traversing up
or down a slope imposes varying risk and e.g. within tight corridors only a
subset of actions can be executed collision-free. The traversability cost
tensor allows to capture this information adequately and measures
traversability in a range of 0 to 1, where 0 is untraversable and 1 risk-free
traversable.

$\mathbf{T}_{\pi_{\theta}}(x,y,z,\alpha_{j},a_{j})\in[0,1]$ (2)

Our objective is to find the function $f_{trav}$ which predicts each element
within the traversability cost tensor $\mathbf{T}$, given the occupancy map
representation $\mathbf{O}\in\mathcal{O}$ and state of the robot
$\mathbf{s}\in\mathcal{S}$.

$\begin{split}\mathbf{T}_{\pi_{\theta}}(x,y,z,\alpha_{i},a_{j})=f_{trav}(\mathbf{s},\mathbf{O})\\\
f_{trav}:\quad\mathcal{S}\times\mathcal{O}\rightarrow\mathbb{R}_{[0,1]}\end{split}$
(3)

###  III-B Simulating Traversability:

The traversability cost is evaluated by simulating locomotion over randomly
generated terrain in NVIDIA’s IsaacGym [16]. To achieve accurate sim-to-real
transfer the robot is controlled by the same locomotion policy later deployed
on the real robot. In this work, the traversability cost is defined as the
success rate of reaching the desired goal pose, which is of main interest for
autonomous navigation. The goal pose is determined by the motion command with
respect to the robots starting pose. Evaluation of multiple trials per motion
command under domain randomization results in a robust geometric
traversability risk assessment. We introduce the terrain generation procedure
in (Sec. III-B1), followed by the traversability data collection in (Sec.
III-B2).

####  III-B1 Terrain Generation

We aim to generate a large and diverse terrain dataset covering the domain of
terrains that the robot may encounter during deployment. Acquisition of a
large amount of experience has shown to increase generalization performance
for various tasks, including locomotion [18] and traversability estimation
[30]. The ground structure of the terrain map is generated using Perlin noise
with either a smooth surface or height-varying steps, following the commonly
used terrain generation for training locomotion policies [18, 23]. To generate
a complex environment with a variety of static obstacles, multiple objects
with a diameter between $0.1\text{\,}\mathrm{m}$ and $6\text{\,}\mathrm{m}$
are randomly spawned in the map. These objects were originally generated by
human artists to create 3D worlds for computer games. In total between 300 and
1000 objects are spawned with a random location and orientation. The spawning
probability of each object is inversely scaled with the diameter of its axis-
aligned bounding box leading to more smaller objects. Each spawned object’s
size is randomly scaled with a factor of $0.5$ and $1.5$. With a probability
of $90\text{\,}\mathrm{\char 37\relax}$ the spawned object is aligned with the
ground surface and otherwise spawned uniformly sampled between
$0\text{\,}\mathrm{m}$ and $3\text{\,}\mathrm{m}$ above the ground resulting
in overhanging obstacles. One example of a $32\text{\,}\mathrm{m}$ $\times$
$32\text{\,}\mathrm{m}$ terrain patch is shown in Figure 2.

####  III-B2 Traversability Data Collection

An initial set of robot starting poses is evaluated by uniformly sampling
points on the generated terrain mesh. We evaluate if the robot can be spawned
collision-free with nominal joint configuration above the sampled point. This
is conducted for all headings $\alpha_{k}$ sampled with $10^{\circ}$
resolution, resulting in a total of 36 possible start headings per point. The
robots’ pitch and roll base orientation are aligned with the underlying mesh.
The base height is adjusted such that all feet are in proximity to the mesh
surface. Only a single admissible robot starting pose for each heading is
further evaluated per voxel with the same spatial resolution as the
traversability cost tensor $\mathbf{T}$. Finally, to verify that the start
position is safe for the robot to be initialized, $2\text{\,}\mathrm{s}$ of
standing is simulated within the physics simulator with the current locomotion
policy. If the robot is capable to remain at the same position successfully
without collisions, the pose is labeled as a valid starting pose. Determining
good starting poses is essential to reduce failure during traversability data
collection due to a bad initialization of the robot. In addition, data
collection is only performed for a reduced subset of initial starting poses in
order to reduce the compute requirements.

For each of the valid spawn pose, $n_{total}=10$ trials of each motion command
$a_{i}$ are executed. In particular, 4 translation and 2 rotation motions are
performed. The translation distance is set to
$40\text{\,}\mathrm{c}\mathrm{m}$ along four cardinal directions (forward,
backward, left, right) and the rotation is given by a
$\pm$$45\text{\,}\mathrm{\SIUnitSymbolDegree}$ yaw. The selected motion
commands determine the relative pose offset with respect to the starting pose.
During execution the current pose offset of the robot with respect to the goal
is used to compute the velocity command tracked by the locomotion policy. Four
different trajectories with start (transparent) and goal poses (dense) with
purple arrows indicating the motion command are shown in Figure 2. Each
attempt results in a robot trajectory, which is evaluated for success based on
foot contacts, base orientation, environment collision, time to reach the goal
and distance to the goal. A maximum execution time of $4\text{\,}\mathrm{s}$
is set for each attempt. In Figure 2 the forward and side-ward moving robots
reach their goal successfully (green check mark), whereas the other robots
fail to execute the desired motion command given a collision with the
environment (red cross). The overall traversability cost is given by the
number of successful attempts $n_{suc}$ divided by total attempts $n_{total}$.

We choose to use relatively small pose offsets, such that a motion command can
be correlated with the terrain surrounding the starting pose. Therefore the
outcome of each trajectory is solely associated with the single starting entry
within the traversability cost tensor, which is indexed by the starting robot
position. Evaluating longer trajectories leads to the problem of determining
which terrain part is responsible for the failed execution. Additionally, it
is possible to increase the safety limits by adjusting the size of the used
collision bodies. For each attempt the robot base mass and friction
coefficient are randomized. The procedure to evaluate a trajectory on success
can be further adapted depending on the application, which may include time,
number of steps, or the configuration of the robot. The collected
traversability cost for an example terrain patch is color-coded in Figure 2,
where blue indicates risk-free traversable and red traversable under high
risk.

###  III-C Sparse 3D Traversability Prediction Network

####  III-C1 Input

We train a sparse 3D CNN to predict traversability information in real-time.
Offline, for each terrain patch, the mesh is converted to a voxel occupancy
map with a resolution of $10\text{\,}\mathrm{c}\mathrm{m}$ using ray tracing.
Every voxel containing a triangle is labeled as occupied. The occupancy map
generated during the mission provides fused information of all past
observations while functioning as an abstraction layer. The previously
computed robot starting poses are used to crop a robot-centric training
example from the larger terrain patch. The robot’s base position is at the
center of the extracted volume. The volume is rotated to align with robot’s
base yaw, while the z-axis is gravity aligned. This allows to implicitly
provide the network the robot’s pose, while the inclination of the terrain
remains unchanged. We choose to use a window size of $8\text{\,}\mathrm{m}$
$\times$ $8\text{\,}\mathrm{m}$ $\times$ $4\text{\,}\mathrm{m}$ or
respectively a volume size of $80\times 80\times 40$ voxels at a
$10\text{\,}\mathrm{c}\mathrm{m}$ resolution. The sparse 3D CNN is implemented
using the Minkowski Engine [3].

####  III-C2 Network

We adapt the generative sparse detection network developed for semantic
segmentation [9]. The network consists of an encoder-decoder structure with
skip connections between the respective feature maps. Each encoder block
applies a sparse 3D convolution (stride=2, kernel=4), batch normalization and
exponential linear unit as an activation function, which downsamples the
spatial resolution by a factor of 2. The encoder generates feature maps with
[8, 16, 32, 64, 128] channels. In contrast to the original network
architecture skip connections are only applied in the two lower layers of the
architecture. The intuition behind the choice of reducing the number of skip
connections is discussed in Section IV-A. The decoder spatially up-samples the
feature maps to the original resolution using generative transposed
convolutional layers. At each decoder level a pruning mask is predicted by a
sparse convolutional layer (stride=1, kernel=1) to increase the sparsity and
reduce the computation needed. Only non-zero elements are propagated to the
following layers. Furthermore, we tested applying our network in three
different output configurations, providing different levels of information
that can be used for path planning purposes later. In the simplest case only
the total risk is assessed, where the commands and robot configuration are
marginalized within the traversability cost tensor resulting in a $c=1$
dimensional target. The orientation risk when marginalizing over the motion
commands results in a $c=18$ dimensional traversability score with a
$10\text{\,}\mathrm{\SIUnitSymbolDegree}$ resolution when accounting for the
symmetry of the robot. Lastly, the c-bin motion direction risk for which the
relative motion direction of the applied action-heading combination, is used
to bin them in $c$-equally sized motion direction bins.

####  III-C3 Training

The network is trained to minimize the reconstruction and traversability
estimation loss. The reconstruction loss is evaluated per decoder layer and is
given by the weighted binary cross-entropy between predicted pruning mask and
the respectively down-sampled target traversability cost tensor. The
traversability loss is given by the mean squared error between the predicted
and target traversability score. Here, the traversability loss is only
assessed for correctly instantiated (true positive) points. ADAM with weight
decay is used to optimize the network for 100.000 steps with a batch size of
8. The learning rate is scheduled according to the 1cycle learning rate with a
target learning rate of $0.001$ [24]. During training at each network layer
the correct-pruning mask is enforced for training stability.

The occupancy map created based on the full available mesh data does not
resemble the data sensed by a robot during its mission. Surfaces are mostly
observed from limited viewpoints, under sensor noise and odometry errors,
which degrade the quality of the occupancy map generated online during a
mission. We increase the robustness of the network using data augmentation.
Firstly, only voxels are marked as traversable if a continuous traversable
path exists, connecting the voxel to the starting robot position at the
center. This is done by applying a flood-fill algorithm. The designed data
augmentation exploits two underlying assumptions. The traversability score is
independent of the occupancy state of a voxel further away than
$1\text{\,}\mathrm{m}$ from the next traversable cell and the traversability
cost prediction is supposed to be invariant to minor surface imperfections.
First, single voxels are dropped randomly with an increasing probability the
further away from the center. Secondly, the surface imperfections are
generated by randomly spawning voxels within the proximity of occupied voxels.

##  IV EXPERIMENTS

###  IV-A Training Results

Figure 3: ETH HPH building test dataset (left to right): Local Voxblox
occupancy map, height is color-coded; Predicted total risk traversability cost
for local occupancy map (blue risk-free traversable, red traversable under
high risk); Simulated total risk traversability cost for full map-patch, (blue
risk-free traversable, red traversable under high risk);

To train our proposed network, we created 20 training and 10 validation
terrain patches with a size of $32\text{\,}\mathrm{m}$ x
$32\text{\,}\mathrm{m}$ following the procedure explained in Section III-A.
Separate object asset datasets are used for generating the training and
validation terrain. These objects were generated by professional artists and
are freely available on the Unity Asset Store. We simulate 4096 robots to
collect traversability information in parallel facilitating an accumulated
real-time speed up factor of 400. Per patch 2500 minutes wall time are spent
collecting traversability using an RTX2080TI and 18 CPU cores of a AMD EPYC
7763. In total 57 years of real-time equivalent traversability information was
gathered.

The test dataset is generated by applying the same traversability simulation
procedure for two real-world meshes consisting of an underground mine in
Switzerland (1 patch) and a multi-floor university building at the ETH Zurich
(4 patches). The mine mesh was collected during a preparatory run of team
CERBERUS [25] for the DARPA Subterranean Challenge. All meshes were captured
using a Leica RTC360 laser scanner. While the training and validation datasets
are drawn from the same distribution, the test dataset allows assessing
performance on real-world data. In the absence of ground truth, the generated
test dataset allows for assessing the network’s capability to reproduce the
simulated traversability cost from the limited geometrical information
provided by the occupancy representation. In Section IV-B we experimentally
validate how this traversability estimate transfers to the real-world and can
be used for path planning. In Table I we provide an overview of the generated
datasets.

TABLE I: Datasets Overview: Assets reefers to the number of assets gathered. Instances is the range of spawned assets per patch.  | Patches | Size [m] | Assets | Instances | Real
---|---|---|---|---|---
Train | 20 | 32x32x16 | 300 | 300-1000 | ✗
Val | 10 | 32x32x16 | 200 | 300-1000 | ✗
Test | 5 | 32x32x16 | - | - | ✓

The root-mean-squared error (RMSE) between the predicted traversability
$\hat{y}$ and simulated traversability score $y$ is used to evaluate the
performance. Not instantiated voxels are assumed to have a traversability of
zero, and hence true negative (TN) voxels are not included when calculating
the RMSE metric.

$\begin{split}L_{RMSE}&=\sqrt{\frac{1}{N}\sum_{TP,FN,FP}(y-\hat{y})^{2}}\\\
N&=|TP|+|FN|+|FP|\end{split}$ (4)

TABLE II: RMSE reported for different network architectures M, output channels c, and datasets (Train, Validation, Test) | $c=1$ |  | $c=4$ |  | $c=18$
---|---|---|---|---|---
| Tr | Val | Test |  | Tr | Val | Test |  | Tr | Val | Test
M1 | 0.124 | 0.131 | 0.101 |  | 0.232 | 0.236 | 0.162 |  | 0.450 | 0.465 | 0.385
M2 | 0.119 | 0.126 | 0.097 |  | 0.221 | 0.226 | 0.155 |  | 0.438 | 0.450 | 0.379

During the training of the network, multiple hyper-parameters, including the
learning rate, model capacity (number of channels), loss function hyper-
parameters, data augmentations, and model architecture, were validated. We
performed experiments for altering the number of skip-connections, which allow
information forwarding between the encoder and decoder layers. Full skip-
connectivity results in a large number of false positively classified voxels
and overall worse performance (M1, Table II). Intuitively this can be
explained by the fact that the decoder has to correctly predict the pruning
mask for all forwarded features by the encoder. Only using skip connections in
the two lower layers resulted in the lowest RMSE (M2, Table II). Furtermore,
with the increase of task complexity (the number of output channels), the RMSE
increases for all network architectures. This increase can be attributed to
two reasons. First, the supervision signal includes more noise for a more
complex task given that fewer robot trajectories are accumulated within a
traversability estimate. Second, the model complexity is kept constant and a
single pruning mask is predicted for all channels within the decoder.

Additionally, it can be noted that the performance on the test data is better
than training data. We suspect that the real world data overall is an easier
objective to estimate than the randomly generated terrain. This observation
also demonstrates that the network is not overfitting to the training dataset.
In Figure 3 we deploy ANYmal, simulated in Gazebo, on the mesh of the
university building at the ETH Zurich. We simulate the LiDAR sensors and use
their measurements to build a robot-centric occupancy map (left). The network
prediction (middle) and simulated traversability cost for the full map are
color-coded. The trained network is conservative on the stairs. In addition,
the ground truth traversability estimate on the stairs is noisy but consistent
over the full staircase.

###  IV-B Real-World Deployment

We deployed ANYmal within various environments to verify the successful sim-
to-real transfer. The robot is equipped with a front and back-facing RoboSense
RS-BPearl LiDAR sensors. First, we illustrate the advantage of using a more
expressive occupancy map representation and compare it to an elevation
mapping-based approach. Second, we deploy the robot within a confined space to
showcase the effectiveness of action aware traversability risk assessment.
Last, we apply the network to data collected within a forest and an
underground mine to showcase the effectiveness of our method for planning in
complex real-world scenarios.

####  IV-B1 Overhanging Obstacles

Figure 4: Comparison of our network to an elevation mapping based
traversability estimation approach. Our approach correctly labels the area
underneath the desk as traversable. Elevation mapping based approach [30] is
unable to produce a correct label for this area due to the corrupted elevation
map. Figure 5: Compromised elevation map in inclined low-ceiling mine shaft.
Traversability prediction using the occupancy map. The ceiling is manually
removed for visualization. Figure 6: Path planning in complex natural
environment based on traversability estimate. Planned path is illustrated in
purple. Forest (e,h,i), Cave (a,b,c,d,f,g). In the cave images the ceiling
voxels are manually removed for illustration purposes.

This experiment is designed to highlight the need for using a more expressive
3D volumetric environment representation than the typically used 2.5D
elevation maps for traversability estimation, especially in the presence of
overhanging obstacles. Figure 4 shows the experimental setup with an elevated
desk placed at the center of the room, allowing the ANYmal to pass underneath
the desk without a special motion or gait. The elevation and occupancy map are
generated using the same LiDAR measurements and robot pose estimates provided
by CompSLAM [13]. The elevation map is generated with a spatial resolution of
$5\text{\,}\mathrm{c}\mathrm{m}$ on a Nvidia Jetson Xavier AGX [17]. The
occupancy map is generated using VoxBlox [20] with a spatial resolution of
$10\text{\,}\mathrm{c}\mathrm{m}$ on the CPU. The occupancy map adequately
captures the free space underneath the desk, while the elevation map of the
same scene is unable to represent the area under the desk as free given that
points on the tabletop are integrated. We compare our work to the
traversability estimation approach proposed in [30], which uses an elevation
map as an input. This network is trained on simulated data for ANYmal using a
different data generation procedure. We use the estimated risk output of [30]
to conduct a qualitative comparison solely. No quantitative results are
presented given that both networks are trained on different data, and no
ground truth information of the real traversability score can be obtained. The
network of [30] is queried for each cell within the elevation map to assess
the risk for 8 different headings when moving forward. In total 204.800 start-
goal pairs per elevation map are evaluated and the average risk is color-coded
per cell in Figure 4. Red indicates untraversable terrain and blue risk-free
traversable terrain. Given the corrupted input representation, the area
underneath the desk is classified as non-traversable, which prohibits planning
underneath the desk. As a comparison, the output of our total risk network
configuration is shown on the left in Figure 4 for the same scene. Here it is
essential to distinguish between the color-coding used. In our work red
implies traversable with a high risk ([30] untraversable) and blue risk-free
traversable ([30] risk-free traversable). Our network predicts a
traversability risk more adequately for the presented scene. Close to
obstacles, the risk increases, and the occupied space is labeled as
untraversable. Additionally, in Figure 5 we illustrate a scenario of ANYmal
traversing in an inclined mine shaft with a low ceiling. Due to low ceiling
clearance the elevation map gets corrupted as depth measurements on the
ceiling are integrated. Our network, on the other hand, utilizing a 3D
occupancy map can adequately predict the traversability cost. The prediction
only includes minor errors (bottom left, predicting traversable region on top
of the mine shaft), which can optionally be filtered out using a flood-fill
algorithm. Here it is important to mention that scenarios like long steep
inclined tunnels are not specifically covered within the training data.

####  IV-B2 Confined Spaces

Figure 7: Deployment of ANYmal within a confined corridor. Individual output
channels for 4-bin motion direction risk estimation.

In the following experiment, we deploy ANYmal within a confined space in an
L-shaped corridor. The experimental setup is visualized from two view-points
in Figure 7. The output of the network is illustrated for the 4-bin motion
direction risk configuration. The bins can be associated with moving up
[$315\text{\,}\mathrm{\SIUnitSymbolDegree}$,$45\text{\,}\mathrm{\SIUnitSymbolDegree}$],
right
[$45\text{\,}\mathrm{\SIUnitSymbolDegree}$,$135\text{\,}\mathrm{\SIUnitSymbolDegree}$],
down
[$135\text{\,}\mathrm{\SIUnitSymbolDegree}$,$225\text{\,}\mathrm{\SIUnitSymbolDegree}$]
and left
[$225\text{\,}\mathrm{\SIUnitSymbolDegree}$,$315\text{\,}\mathrm{\SIUnitSymbolDegree}$].
A higher risk is associated with actions that lead to a collision with the
environment. Within the corridor entry, the cost of moving down and up
strongly depends on the robot’s position. This additional information can be
helpful for safe path planning and allows to limit or rank motion commands at
a spatial location within a terrain according to the excepted risk. During
path planning this additional information can guide the search and may improve
the planning performance in terms of risk or reduce planning time.

####  IV-B3 Path Planning in Complex and Natural Environments

To illustrate the real-world application of our proposed method, we deployed
it on data gathered within a forest and the data from the DARPA subterranean
challenges final run of team CERBERUS. We use the total risk network
configuration to predict the traversability based on the generated occupancy
map. Within the forest and underground cave, no non-rigid obstacles as e.g.
high-grass or bushes, are present. Therefore the generated occupancy map
representation correctly reflects the support surfaces.

To showcase the application of our traversability prediction network we
implement a simple local path planner capable to find a safe path between the
current robot position at the center and an externally provided target
position. Given the predicted traversability cost a directional graph is
generated. Each traversable voxel is node within the graph connected with its
26-neighborhood by an edge with a cost of $c_{s-g}$:

$c_{s-g}=\left\|\mathbf{p}_{g}-\mathbf{p}_{s}\right\|_{2}+\lambda(1-\hat{\mathbf{T}}(\mathbf{p}_{g})),$
(5)

where $\mathbf{p}_{s/g}$ denotes the position of the start and goal node. The
hyper-parameter $\lambda$ trades off the length of the path and the total
risk. By increasing $\lambda$, lower total risks paths are preferred at the
cost of a geometrical longer path between start and goal node.

The Dijkstra algorithm is used to find the shortest path between the current
robot position and the provided target position. The resulting path for the
forest and cave environment with $\lambda=0.1$ is illustrated in purple for
all 5 example scenes in Figure 6. Within the forest, the generated path
traverses with a safety margin around all nearby trees. In the cave
environment, the generated path traverses along the center of the corridor,
maximizing the distance between both walls. Finally, we used the dataset
collected within the forest to assess the average inference time of the
network on a Nvidia Jetson Xavier AGX. The inference time depends on the
number of occupied voxels. On average, the network runs at
$16.5\text{\,}\mathrm{h}\mathrm{z}$$\pm$$1.2\text{\,}\mathrm{h}\mathrm{z}$
processing 12181$\pm$1611 input voxels. The inference time is measured over
100 selected key-frames in the forest dataset. Given the update rate of each
LiDAR is $10\text{\,}\mathrm{h}\mathrm{z}$ the network is capable for real-
time operation.

##  V CONCLUSION

In this work we showed that the proposed sparse 3D convolution neural network
trained on traversability data collected in simulation on randomly generated
terrains can accurately predict the traversability for unseen terrains using a
fully geometric occupancy representations as input. Therefore the network
trained exclusively on simulated data generalizes successfully to the real
world. It successfully handles overhanging obstacles, outperforming elevation
map based approaches in these scenarios, while being computational efficient
to run in real-time. To achieve this we implemented a pipeline to gather
locomotion policy-specific traversability information in simulation. For
future work we will investigate integration of visual semantic information to
assess friction coefficients or the quality of individual footholds.
Additionally, no temporal consistency is enforced between individual
predictions in the current implementation. This leads to changing predictions
for individual voxels over time. Expanding the created robot-specific
traversability estimation module to different robot categories will allow to
develop a unified planning system operating across robot types without
requiring robot-specific heuristics. Finally, we will investigate the utility
of the generated traversability estimates for a variety of downstream tasks,
including active switching between locomotion policies or gait patterns, local
and global path planning, and exploration.

## References

  * [1] A. Bouman and et. al, “Autonomous spot: Long-range autonomous exploration of extreme environments with legged locomotion,” _CoRR_ , vol. abs/2010.09259, 2020.
  * [2] R. O. Chavez-Garcia, J. Guzzi, L. M. Gambardella, and A. Giusti, “Learning ground traversability from simulations,” _CoRR_ , vol. abs/1709.05368, 2017\.
  * [3] C. B. Choy, J. Gwak, and S. Savarese, “4d spatio-temporal convnets: Minkowski convolutional neural networks,” _CoRR_ , vol. abs/1904.08755, 2019.
  * [4] D. D. Fan and et. al, “Step: Stochastic traversability evaluation and planning for risk-aware off-road navigation,” _Robotics: Science and Systems XVII_ , 2021.
  * [5] P. Fankhauser and et. al, “Robot-centric elevation mapping with uncertainty estimates,” in _International Conference on Climbing and Walking Robots (CLAWAR)_ , 2014.
  * [6] J. Frey, H. Blum, F. Milano, R. Siegwart, and C. Cadena, “Continual learning of semantic segmentation using complementary 2d-3d data representations,” _CoRR_ , vol. abs/2111.02156, 2021.
  * [7] B. Graham, M. Engelcke, and L. van der Maaten, “3d semantic segmentation with submanifold sparse convolutional networks,” _CoRR_ , vol. abs/1711.10275, 2017.
  * [8] J. Guzzi, R. O. Chavez-Garcia, M. Nava, L. M. Gambardella, and A. Giusti, “Path planning with local motion estimations,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 2586–2593, 2020.
  * [9] J. Gwak, C. B. Choy, and S. Savarese, “Generative sparse detection networks for 3d single-shot object detection,” in _European conference on computer vision_ , 2020.
  * [10] L. Han, F. Gao, B. Zhou, and S. Shen, “FIESTA: fast incremental euclidean distance fields for online motion planning of aerial robots,” _CoRR_ , vol. abs/1903.02144, 2019.
  * [11] A. Hornung, K. M. Wurm, M. Bennewitz, C. Stachniss, and W. Burgard, “OctoMap: An efficient probabilistic 3D mapping framework based on octrees,” _Autonomous Robots_ , 2013.
  * [12] L. D. Jackel and et. al, “The darpa lagr program: Goals, challenges, methodology, and phase i results,” _Journal of Field Robotics_ , vol. 23, no. 11-12, pp. 945–973, 2006.
  * [13] S. Khattak and et al., “Complementary multi–modal sensor fusion for resilient robot pose estimation in subterranean environments,” in _IEEE International Conference on Unmanned Aircraft Systems_ , 2020.
  * [14] P. Krüsi and et al., “Driving on point clouds: Motion planning, trajectory optimization, and terrain assessment in generic nonplanar environments,” _Journal of Field Robotics_ , vol. 34, pp. 940–984, 2017.
  * [15] M. Kulkarni and et. al, “Autonomous teamed exploration of subterranean environments using legged and aerial robots,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2022.
  * [16] V. Makoviychuk _et al._ , “Isaac gym: High performance GPU based physics simulation for robot learning,” in _Conference on Neural Information Processing Systems (NeurIPS)_ , 2021.
  * [17] T. Miki and et al., “Elevation mapping for locomotion and navigation using gpu,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2022), under review_ , 2022.
  * [18] T. Miki, J. Lee, J. Hwangbo, L. Wellhausen, V. Koltun, and M. Hutter, “Learning robust perceptive locomotion for quadrupedal robots in the wild,” _Science Robotics_ , vol. 7, no. 62, p. eabk2822, 2022.
  * [19] K. Museth, “Vdb: High-resolution sparse volumes with dynamic topology,” _ACM Trans. Graph._ , vol. 32, no. 3, 2013.
  * [20] H. Oleynikova and et. al, “Voxblox: Incremental 3d euclidean signed distance fields for on-board mav planning,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2017.
  * [21] K. Otsu and et. al, “Autonomous terrain classification with co- and self-training approach,” _IEEE Robotics and Automation Letters_ , vol. 1, no. 2, pp. 814–819, 2016.
  * [22] P. Papadakis, “Terrain traversability analysis methods for unmanned ground vehicles: A survey,” _Engineering Applications of Artificial Intelligence_ , vol. 26, no. 4, pp. 1373–1385, 2013.
  * [23] N. Rudin and et. al, “Learning to walk in minutes using massively parallel deep reinforcement learning,” in _Proceedings of the 5th Conference on Robot Learning_ , 2022, pp. 91–100.
  * [24] L. N. Smith and N. Topin, “Super-convergence: Very fast training of residual networks using large learning rates,” _CoRR_ , vol. abs/1708.07120, 2017\.
  * [25] M. Tranzatto and et. al, “Cerberus: Autonomous legged and aerial robotic exploration in the tunnel and urban circuits of the darpa subterranean challenge,” 2022.
  * [26] P. Wang, Y. Liu, Y. Guo, C. Sun, and X. Tong, “O-CNN: octree-based convolutional neural networks for 3d shape analysis,” _CoRR_ , vol. abs/1712.01537, 2017.
  * [27] L. Wellhausen and et. al, “Where should i walk? predicting terrain properties from images via self-supervised learning,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 2, pp. 1509–1516, 2019.
  * [28] L. Wellhausen and M. Hutter, “Rough terrain navigation for legged robots using reachability planning and template learning,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2021.
  * [29] M. Wermelinger and et. al, “Navigation planning for legged robots in challenging terrain,” in _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2016, pp. 1184–1189.
  * [30] B. Yang and et. al, “Real-time optimal navigation planning using learned motion costs,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 9283–9289.