# Anomalous Coherence Length of Majorana Zero Modes at Vortices in
Superconducting Topological Insulators

Bo Fu and Shun-Qing Shen  Department of Physics, The University of Hong Kong,
Pokfulam Road, Hong Kong, China

(September 3, 2024)

###### Abstract

The coherence length of two Majorana zero energy modes in a p-wave topological
superconductor is inversely proportional to the superconducting order
parameter. We studied the finite size effect of the Majorana zero modes at
vortices in a topological insulator/superconductor heterostructure in the
presence of a vortex and found that the the coherence length of the two zero
energy modes at the terminals of a vortex line is independent of
superconducting order parameter, and determined by the intrinsic properties of
the topological insulator. This anomalous property illustrates that the
superconducting topological insulator is topologically distinct, contrary to a
$p$-wave topological superconductor.

#### Introduction

Search for Majorana zero modes in topological phases has generated extensive
interest in condensed matter physics and material science (Wilczek-09np, ;
Stern-10nature, ; Qi-11rmp, ; Alicea-12rpp, ; Shen-12book, ;
Beenakker-13arcmp, ; Sato-16jpsj, ; Lutchyn-18nrm, ). The Majorana zero mode
in topological superconductor carry zero energy and obey non-Abelian
statistics. Their occupancy can form the topological degeneracy of the ground
states of the system, which are expected to have potential application for
fault tolerant topological quantum computation (Kitaev-03ap, ; Freedman-03, ;
Nayak-08rmp, ). In their pioneering work, Fu and Kane (Fu-08prl, ) proposed
that the proximity effect between an s-wave superconductor and the surface
electrons of a strong topological insulator leads to a time-invariant
superconducting state resembling a spinless $p_{x}\pm ip_{y}$ superconductor.
They proposed that this interface supports Majorana bound modes at vortices.
Over the last decade, this proposal has attracted significant attention and
become one of the main prototypes to construct and to engineer a physical
system to host the topological excitations (Lutchyn-10prl, ; Oreg-10prl, ;
Sau2010prl, ; Sau10prb, ; Mourik-15science, ; Das-12np, ; Nadj-
Perge-14science, ; XuJP-15prl, ; Sun-16prl, ; Hu-16prb, ), and to understand
the zero energy modes observed in iron-based superconductors
(Zhang_-18Science, ; WangDF-18Science, ; Machida-19nm, ; KongLY, ; LiuQ-18prx,
; LiuWY-20nc, ; ZhuSY-20Science, ; Chiu-20sa, ; LiM-22nature, ; Konig-19prl, ;
zhang-21prl, ; Kheirkhah-21prb, ; Hu-arxiv, ). However the search for Majorana
zero modes is meeting a great of difficulty and challenge especially in
experiments.

The time-reversal-invariant superconductor with spin-orbit coupling belongs to
symmetry class DIII and in two-dimensions is characterized by a
$\mathbb{Z}_{2}$ topological invariant(Schnyder-08prb, ; Kitaev-09aip, ;
Chiu-16rmp, ). The topologically nontrivial phase hosts a pair of helical
Majorana modes on its edge (Fu2010odd, ; Qi-10prb, ; Zhang-13prl, ;
Haim-19PhysRep, ). Theoretically, it can be realized by considering the spin-
triplet superconducting pairing with odd-parity or extended-s wave pairing
which flips its sign when it evolves across the Brillouin zone. Based on the
odd-parity superconductivity criterion (Fu2010odd, ; Sato2010Topological, ),
the strong topological insulator in contact with an s-wave superconductor is
topologically trivial without edge modes, which is distinctly different from
the chiral $p$-wave topological superconductor (Green-00prb, ; Ivanov-01prl,
). The existence of the zero energy vortex mode in this system is associated
with the Atiyah-Singer index theorem which clarifies the correspondence
between the vorticity of the vortex and the number of the localized zero-
energy mode for the surface states (Atiyah-63bam, ; Jackiw-81npb, ;
Jackiw-07prl, ). It heavily relies on the validity of the topological
insulator’s surface states and the presence of the chiral symmetry.
Furthermore, the chemical potential enters into the Bogoliubov-de Gennes (BdG)
equation in a nontrivial way and breaks the chiral symmetry explicitly that
the index theorem does not apply here. There arises the question how the
tunneling between two surfaces lifts the degeneracy of the Majorana modes in a
thin film for a finite chemical potential.

In the present work, we investigate the finite size effect of the Majorana
zero energy modes in vortices in the topological insulator/superconductor
(TI/SC) heterostructure depicted as Fig. 1(a). In the presence of the
superconducting vortex the two Majorana modes are present and connected
through the bulk topological insulator along the vortex when the chemical
potential $\mu$ is lower than a critical value $\mu<\mu_{c}$. The energy
splitting of the two modes decays exponentially with the thickness and the
coherence length only depends on the intrinsic properties of the topological
insulator, and is independent of the superconducting order parameter. As a
comparison, we also present the results for the semimagnetic topological
insulator/ superconductor (SMTI/SC) heterostructure depicted as Fig. 2(a),
which is equivalent to a $p$-wave topological superconductor when the chemical
potential locates within the magnetic gap of the top surface states. The two
Majorana modes reside at the vortex core and at the boundary separately, and
their coherence length is equal to the superconducting coherence length, which
is a typical signature of a p-wave topological superconductor. Thus the
anomalous coherence length in the TI/SC heterostructure indicates that the the
pair of the zero energy modes at the vortex core is attributed to the winding
number of the superconducting order parameter, not to the p-wave topological
superconductivity.

Figure 1: The TI/SC heterostructure in the presence of a superconducting
vortex. (a) Schematic with two zero energy modes bound to a vortex line.
$L_{z}$ denotes the thickness of the sample. (b) $\mathbb{Z}_{2}$ topological
invariant as a function of $\mu$. (c) The dispersions of quasi-1D system at
different values of $\mu$. The blue-to-red color gradient indicates the radial
probability distribution for each band with the vortex line at the origin. The
blue colored bands correspond the dispersions for the vortex bound states. (d)
Energy spectrum and the wave-function for the zero energy modes for open
boundary conditions. (e) Plot of the finite size induced gap between the
surface states at $k_{x}=k_{y}=0$ as a function of the thickness
$L_{z}=N_{z}a$ with $a$ as the lattice constant. (f) For $\mu=0$, the
evolution of the energy of the vortex bound states with respect to the
thickness $N_{z}$ (g) The the energies of the vortex states as a function of
$\mu$. For $N_{z}=8$, $N_{x}=N_{y}=50$ and for $N_{z}=9,10$, $N_{x}=N_{y}=40$.
Parameters are $\Delta=0.05$, $\hbar v=0.4$, $B=0.5$, and $m=0.28$.

#### The TI/SC heterostructure

We start with a minimal bulk model for a three dimensional topological
insulator $H_{\mathbf{k}}$ which supports gapless surface states (Qi-11rmp, ).
The numerical calculation is based on a tight-binding model on a cubic
lattice, and the analytical study is based on the continuum model in the long
wavelength approximation that the lattice Hamiltonian is expanded in terms of
the wave vector $\mathbf{k}$ to the second order around $\Gamma$ point,
$H_{\mathbf{k}}=v\rho_{x}\mathbf{k}\cdot\boldsymbol{\sigma}+M(\mathbf{k})\rho_{z}-\mu$
where $M(\mathbf{k})=m-B\mathbf{k}^{2}$, the Pauli matrices $\sigma$ and
$\rho$ are acting on spin and orbit space respectively(ZhangH-09np, ;
Qi-11rmp, ; Shen-12book, ), $v,m,B$ are the material parameters and $\mu$ is
the chemical potential. In proximity to an $s$-wave superconductor, a finite
superconducting pairing $\Delta$ is induced in the three-dimensional
topological insulator which leads to the BdG Hamiltonian
$H_{\mathbf{k}}^{\mathrm{BdG}}$. Based on the odd-parity superconductivity
criterion, $H_{\mathbf{k}}^{\mathrm{BdG}}$ is topologically trivial without
gapless edge modes around the system (Fu2010odd, ). In the presence of a
vortex in $z$ direction, $\Delta\to\Delta(r)e^{i\theta}$ with $(r,\theta)$ as
the in-plane polar coordinates with respect to the vortex core, the
translational invariance along $z$ direction persists and $k_{z}$ is still a
good quantum number. The problem becomes to classifying the gapped phases in
quasi-1D whose unit cell consists all the sites in $xy$ plane. Due to the
lacking of time reversal symmetry, it belongs to symmetry class D and is
characterized by a $\mathbb{Z}_{2}$ invariant $\nu$ which is defined as the
product of signs of Pfaffians of the antisymmetric and real BdG Hamiltonians
in the Majorana representation at two time reversal invariant momenta
$k_{z}=0,\pi$ (Kitaev-09aip, ; Wimmer-12ACM, ). With increasing $\mu$, the
quasi-1D system transitions into the trivial phase via a quantum critical
point $\mu=\mu_{c}\simeq v\sqrt{\frac{m}{B}}$ at which $\nu$ changes from $-1$
to $+1$ as shown in Fig. 1(b). In Fig. 1(c), we also present the numerical
results of the evolution of the dispersion under the variation of $\mu$. The
results reveal that the vortex line is fully gapped except at $\mu=\mu_{c}$
where a vortex phase transition takes place(Hosur-11prl, ). For $\nu=-1$
($\mu<\mu_{c}$), the quasi-1D system is topological nontrivial and there
exists a single 0D zero mode at each end of the termination along $z$
direction(Hosur-11prl, ; Chiu-11prb, ) as shown in Fig. 1(d). In order to
estimate the energy splitting of two vortex line end states for finite
thickness $L_{z}$, we derive the effective Hamiltonian for the vortex line. By
projecting onto the two states centered at the vortex line which is denoted by
the darkest of blue color in Fig. 1(c), we obtain the effective dispersions
(See Appendix B)

$H_{\mathrm{eff}}=-\mathcal{F}(k_{F}^{2}\xi^{2})\left[(\widetilde{m}-Bk_{z}^{2})\nu_{z}+vk_{z}\nu_{y}\right]$
(1)

where $\widetilde{m}=m-Bk_{F}^{2}-\frac{B}{\xi^{2}}$ is the renormalized mass
with $k_{F}=\mu/v$ and $\xi=v/\Delta$, $\nu_{i}$ are Pauli matrices acting on
the projected two bands and
$\mathcal{F}(x)=[-\frac{\mathrm{E}(-x)}{x+1}+\mathrm{K}(-x)]/[\mathrm{E}(-x)-\mathrm{K}(-x)]$
is monotonically decreasing function where $\mathrm{K}$ and $\mathrm{E}$ are
the complete elliptic integral of the first and second kind, respectively. The
superconducting pairing only enters into an overall energy renormalization
function $\mathcal{F}$ which can be factored out without changing any of the
topological properties. Then, the effective Hamiltonian (1) resembles the 1D
Su–Schrieffer–Heeger model (Su-80prb, ) instead of the 1D topological Kitaev
chain in which that the p-wave pairing is linear in momentum (Kitaev06ap, ).
Consequently, when $\widetilde{m}B>0$, the vortex line Hamiltonian is
topologically nontrivial. At $\widetilde{m}=0$ which corresponds a critical
chemical $\mu=\mu_{c}$, the gap vanishes signaling a topological phase
transition which is consistent with the $\mathbb{Z}_{2}$ topological
invariant. From Eq. (1), we can obtain the energy splitting for the zero
energy modes

$\delta
E=\frac{4\widetilde{m}v\mathcal{F}(k_{F}^{2}\xi^{2})}{\sqrt{4B\widetilde{m}-v^{2}}}\left|\sin\left(\frac{\sqrt{4\widetilde{m}B-v^{2}}}{2B}L_{z}\right)\right|\exp\left(-\frac{vL_{z}}{2B}\right).$
(2)

For the chemical potential $\mu\sim 0$ or in the strong pairing limit
$k_{F}\xi\ll 1$, we have $\mathcal{F}(k_{F}^{2}\xi^{2})\simeq 1$ and
$\widetilde{m}\simeq m$, the energy splitting is independent on the
superconducting pairing and recovers the finite size effect for the surface
state of topological insulator (Zhou-08prl, ; Linder-09prb, ). As shown in
Fig. 1(f), the energy splitting of the zero energy states based on a tight-
binding numerical calculations quickly saturates when the size of the slab is
much large than the superconducting coherence length $L_{x},L_{y}\gg\xi$ and
features an oscillating exponential decay with increasing the thickness of the
sample, which agrees well with the analytic expression. For $\mu\neq 0$, since
the chemical potential enters into oscillating function
$\sin\left(\frac{\sqrt{4\widetilde{m}B-v^{2}}}{2B}L_{z}\right)$ through the
renormalized mass $\widetilde{m}$, the energy splitting is also sensitive to
$\mu$ besides thickness of the sample as shown in Fig. 1(g). For large $\mu$
or the weak pairing limit $k_{F}\xi\gg 1$,
$\mathcal{F}(k_{F}^{2}\xi^{2})\simeq\frac{\ln(4k_{F}\xi)-1}{k_{F}^{2}\xi^{2}}$
shows a power-law decrease of $k_{F}\xi$. Due to the presence of this pre-
factor, the finite size effect is strongly suppressed.

This energy splitting can be understood from aspect of top and bottom surface
states with superconducting pairing by means of the index theorem
(Weinberg-81prd, ; Fukui-10jpsj, ; Roy-14prb, ). The Dirac surface states of
strong topological insulator thin films can be described by
$h_{\mathbf{k}_{\shortparallel}}^{\mathrm{surf}}=v\varrho_{z}(\mathbf{k}_{\shortparallel}\times\boldsymbol{\sigma})_{z}+t\varrho_{x}-\mu$
where $\mathbf{k}_{\shortparallel}=(k_{x},k_{y},0)$ denotes the in-plane wave
vector, $t$ is the inter-surface tunneling and the Pauli matrices $\varrho$
and $\sigma$ denote surface and spin degrees of freedom, respectively.. In
combination with the superconducting pairing, the BdG Hamiltonian is
$H_{\mathrm{BdG}}^{\mathrm{surf}}(\mathbf{k}_{\shortparallel})=\tau_{z}h_{\mathbf{k}_{\shortparallel}}^{\mathrm{surf}}+\Delta\tau_{x}$
which belongs to the Altland-Zirnbauer symmetry class DIII and is classified
by a $\mathbb{Z}_{2}$ topological invariant. This system has the mirror
symmetry $M_{z}=i\varrho_{x}\sigma_{z}$ with $M_{z}^{2}=-1$ which reflects the
top surface to the bottom surface. After a uniform $\pi/2$ rotation around
$\varrho_{y}\sigma_{z}$, the full BdG Hamiltonian can be decoupled into the
direct sum of two mirror sectors $H_{\mathbf{k}_{\shortparallel}}^{\chi}$ with
mirror eigenvalue as $i\chi$. The two subblocks are particle-hole partners of
each other
$\tau_{y}\sigma_{x}H_{\mathbf{k}_{\shortparallel}}^{\chi*}\sigma_{x}\tau_{y}=-H_{-\mathbf{k}_{\shortparallel}}^{-\chi}$.
Each subblock breaks particle-hole symmetry explicitly and possesses the
chiral symmetry $\\{\mathcal{C},H_{\mathbf{k}_{\shortparallel}}^{\chi}\\}=0$
with $\mathcal{C}=\tau_{y}$ thus belongs to the class AIII. In two spatial
dimensions, the topological classification for class AIII is trivial
(Schnyder-08prb, ). In the presence of a vortex, the pertinent Hamiltonian for
$\mu=0$ becomes

$H^{\chi}=\chi\tau_{z}\left[-iv\left(\sigma_{y}\partial_{x}-\sigma_{x}\partial_{y}\right)+t\sigma_{z}\right]+\Delta\left(\cos\theta\tau_{x}+\sin\theta\tau_{y}\right)$
(3)

where $\mathbf{k}_{\shortparallel}$ is replaced by
$-i(\partial_{x},\partial_{y})$. Note that the five four-dimensional Hermitian
matrices anticommute with each other and the interface tunneling term enters
into the Hamiltonian as the fifth anti-commuting matrix. For $t=0$, there is
additional chiral symmetry $\tau_{z}\sigma_{z}$ in $H^{\chi}$ which can be
expressed as $\\{\tau_{z}\sigma_{z},H^{\chi}\\}=0$ and ensures the spectral
symmetry. As a consequence, the zero energy states $|\Psi_{0}^{\chi}\rangle$
of $H^{\chi}$ becomes eigenstates of $\tau_{z}\sigma_{z}$ with eigenvalue as
$+1$. The analytic index of the chiral symmetric model is defined by
$\mathrm{ind}H^{\chi}=n_{+}-n_{-}$where $n_{\pm}$ are the number of zero-
energy states with chirality $\pm 1$. The index theorem states that the
analytic index is identical to the winding number of the the order parameter
in the two-dimensional space and there are exactly $n$ number of zero modes
for the vorticity $n$ (Weinberg-81prd, ; Fukui-10jpsj, ; Roy-14prb, ). Also as
pointed in Ref. (Teo-10prl, ), the zero modes are associated with hedgehogs in
the complex vector fields of the superconducting order parameter
$\mathbf{n}(\mathbf{r})=(\Delta\cos\theta,-\Delta\sin\theta)$. In particular,
when the vorticity is one there exist single state at zero energy for
$H^{\chi}$. After including the chirality symmetry breaking term $\chi
t\tau_{z}\sigma_{z}$, we have $(H^{\chi}+\chi
t\tau_{z}\sigma_{z})|\Psi_{0}^{\chi}\rangle=\chi
t\tau_{z}\sigma_{z}|\Psi_{0}^{\chi}\rangle=\chi t|\Psi_{0}^{\chi}\rangle$ that
the energy will be shifted from zero to $\chi t$. The existence of the zero
energy solution heavily relies on the assumption that the intersurface
tunneling is negligible. When the tunneling effects are taken into account,
the zero energy bound states are actually shifted away from zero.

Figure 2: The SMTI/SC heterostructure. (a) Schematic of the zero energy modes
residing in the vortex core and the boundary. (b) In the absence of vortex,
the Chern number $\mathcal{N}$ of the quasi-2D system as a function of $\mu$.
(c) The quasi-one-dimensional band structure for $\mu$ for $N_{z}=10$ and
$N_{y}=100$. The distribution of the wave function along the $z$ direction
indicated by the scale from blue to yellow color. (d) The energy spectra for
the open boundary condition and the illustration of the spatial probability
distribution $|\Psi(x,y,z)|^{2}$ of zero energy state and the chiral state for
a $30\times 30\times 18$ lattice. (e) At $\mu=0$, the energy for the lowest
energy mode as a function of $N_{z}$ for different $N_{x}$ and $N_{y}$. (f)
For finite $\mu$, the energy for the lowest energy mode as a function of
$N_{x}$ and $N_{y}$ for $N_{z}=10$. The singularity points in the logarithmic
plot due to the oscillation are indicated by the vertical dashed lines.
Parameters are $\Delta=0.05$, $\hbar v=0.4$, $b=0.5$, and $m=0.28$. The
exchange field $V_{z}=0.1$ are added to top three layers.

#### The SMTI/SC heterostructure

For comparison, we now turn to the SMTI/SC heterostructure. The exchange
interaction between the magnetic ion and the surface electrons leads to
nonzero magnetization and makes the top surface electrons open an energy gap
$2|V_{z}|$. When the Fermi level intersects gapless bottom surface states and
locates within the magnetic gap of the top surface states, i. e.
$|\mu|<\sqrt{V_{z}^{2}-\Delta^{2}}$, the quasi-2D system is topologically
equivalent to a chiral topological superconductor with nonzero Chern number
supporting chiral Majorana modes on its boundary (Fu-08prl, ; Akhmerov-09prl,
; Fu-09prl, ; qi-10prb-b, ; Teo-10prb, ; Qi-13prb-c, ; Wang-15prb, ;
wang-16prb, ; Lian-18PNAS, ; Yan2021majorana, ). A topological phase
transition occurs at $|\mu|=\sqrt{V_{z}^{2}-\Delta^{2}}$ in Fig. 2(b). It is
verified by numerical calculations for different $\mu$ as shown in Fig. 2(c).
The presence of edge states is consistent with the bulk band topology. After
introducing a vortex, the existence of the Majorana zero energy state is
governed by the BdG Hamiltonian for the surface states. We then map the
surface states onto a 2D plane with the intersection point of the vortex line
with the bottom surface mapped to the center of the plane. In this situation,
the exchange field only exists outside a disk radius $R$, i. e.
$M(r)=V_{z}\Theta(r-R)$. For $\sqrt{V_{z}^{2}-\Delta^{2}}>\mu>0$, we find a
zero energy solution $|\psi_{\mathrm{core}}\rangle$ localized at the vortex
core and a bound state $|\psi_{\mathrm{inter}}\rangle$ localized at the
interface between the magnetic and nonmagnetic regions (Ivanov-01prl, ;
Green-00prb, ; Law-09prl, ). We can construct the approximate eigenstate wave
function as
$|\Psi_{\pm}\rangle=\frac{1}{\sqrt{2}}\left(|\psi_{\mathrm{core}}\rangle\pm|\psi_{\mathrm{inter}}\rangle\right)$
with the energies $E_{+}=-E_{-}=\delta E$ (Cheng-09prl, ; Cheng-10prb, ).
These two wave functions satisfy the particle-hole symmetry of the BdG
equations, $\Xi|\Psi_{+}\rangle=|\Psi_{-}\rangle$. As shown in Fig. 2(d),
based on the tight binding calculations, we present the energy spectra and the
wave functions for the lowest two energy states. The zero energy state is
constituted by two parts of contributions: one part is exponentially localized
at the vortex core while the other part is localized at the interface. In
addition to the zero modes, there are chiral modes peaked only at the
interface within the superconducting gap. By considering the overlapping of
two zero energy modes, the energy splitting can be obtained as (Appendix C)

$\delta E\approx 2\Delta e^{-R/\xi}|\sin(k_{F}R-\delta)|$ (4)

where
$\delta=\mathrm{arctan}\sqrt{\frac{V_{z}+\mu}{V_{z}-\mu}}+\frac{\pi}{4}$. The
energy splitting decays exponentially as a function of $R$. The coherence
length is simply the superconducting coherence length $\xi=v/\Delta$, which is
proportional inversely to the superconducting order parameter $\Delta$. The
analytic result (4) is in qualitative agreement with the numerical results as
shown in Fig. 2(e) and (f).

In order to gain better insight into the difference between the two present
cases, we derive an effective model to capture the main physics in the SMTI/SC
heterostructure. We start from the effective Hamiltonian for the strong
topological insulator thin films in contact with the magnetic insulator on its
top,
$h_{\mathbf{k}_{\shortparallel}}^{\mathrm{mag}}=-v(\mathbf{k}_{\shortparallel}\times\boldsymbol{\sigma})_{z}+M(\mathbf{k}_{\shortparallel})\sigma_{z}$
(Mogi-21np, ; Zou-22-arxiv, ). At low energy ($k_{\shortparallel}<k_{c}$ with
$k_{c}\simeq\sqrt{\frac{m}{B}}$), $M(\mathbf{k}_{\shortparallel})=0$,
$h_{\mathbf{k}_{\shortparallel}}^{\mathrm{mag}}$ turns out to be the massless
Dirac Hamiltonian which describes the bottom surface states. At high energy
regime ($k_{\shortparallel}>k_{c}$), $M(\mathbf{k}_{\shortparallel})\neq 0$
originates from the surface states evolve into the bulkin the high energy
regime and breaks time reversal symmetry explicitly. Thus,
$M(\mathbf{k}_{\shortparallel})$ behaves as a regularization term and change
the band topology. The wave function with the positive eigenvalue can be
solved as
$|\psi_{c}\rangle=\left(i\cos\frac{\varphi_{\mathbf{k}_{\shortparallel}}}{2},e^{i\theta_{\mathbf{k}_{\shortparallel}}}\sin\frac{\varphi_{\mathbf{k}_{\shortparallel}}}{2}\right)$
for $M(\infty)>0$ and
$|\psi_{c}\rangle=\left(ie^{-i\theta_{\mathbf{k}_{\shortparallel}}}\cos\frac{\varphi_{\mathbf{k}_{\shortparallel}}}{2},\sin\frac{\varphi_{\mathbf{k}_{\shortparallel}}}{2}\right)$
for $M(\infty)<0$ where
$\cos\varphi_{\mathbf{k}_{\shortparallel}}=M(\mathbf{k}_{\shortparallel})/\epsilon\text{($\mathbf{k}_{\shortparallel}$)}$
with
$\epsilon\text{($\mathbf{k}_{\shortparallel}$)}=\sqrt{v^{2}k_{\shortparallel}^{2}+M^{2}(\mathbf{k}_{\shortparallel})}$.
Note that the angular factor $e^{\pm i\theta}$ must accompany with the
component vanishing at $k_{\shortparallel}\to\infty$ to ensure the single-
valueness of the wave function. After the inclusion of the superconducting
pairing, the BdG Hamiltonian can be projected onto the two bands which are
intersected with the Fermi energy,

$H_{\mathrm{\mathrm{BdG}}}^{\mathrm{mag}}=\left(\begin{array}[]{cc}\epsilon\text{($\mathbf{k}_{\shortparallel}$)}-\mu&\frac{vk_{\shortparallel}}{\epsilon\text{($\mathbf{k}_{\shortparallel}$)}}\Delta
e^{-\mathrm{sgn}[M(\infty)]i\theta_{\mathbf{k}_{\shortparallel}}}\\\
\frac{vk_{\shortparallel}}{\epsilon\text{($\mathbf{k}_{\shortparallel}$)}}\Delta
e^{\mathrm{sgn}[M(\infty)]i\theta_{\mathbf{k}_{\shortparallel}}}&-\epsilon\text{($\mathbf{k}_{\shortparallel}$)}+\mu\end{array}\right).$
(5)

It is an effective chiral $p$-wave BdG Hamiltonian and the chirality crucially
depends on the sign of the $M(\infty)$. The projected Hamiltonian cannot be
determined without ambiguity in the absence of the regulator. The Chern number
for the occupied band is $\mathcal{N}=\mathrm{sgn}[M(\infty)]$. The presence
of a vortex will leads to the anti-periodical condition for the wave function,
$\Psi(r,\theta+2\pi)=e^{i\pi}\Psi(r,\theta)$(Green-00prb, ; Ivanov-01prl, ).
By using the Ansatz for the wave function
$\Psi_{l}(\mathbf{r})=\frac{e^{il\theta}}{\sqrt{2\pi
r}}[e^{-i\theta/2}u_{l}(r),e^{i\theta/2}v_{l}(r)]$, we can obtain the the
radial BdG Hamiltonian for $[u_{0}(r),v_{0}(r)]$

$H_{\mathrm{BdG}}^{\mathrm{radial}}=[\epsilon(-\partial_{r}^{2})-\mu]\tau_{z}-i\frac{\Delta}{2\mu}v\partial_{r}\tau_{x}$

which is equivalent to 1D Kitaev chain. The superconducting pairing only
enters into the off-diagonal terms in sharp contrast with Eq. (1).

#### Summary

The pairing patterns of the Majorana modes in the TI/SC and SMTI/SC
heterostructure in the presence of a vortex have different features. Both
numerical simulation and analytical analysis show that the energy splitting of
the modes decays exponentially in the thickness $L_{z}$ of the TI layer in the
TI/SC structure and the size $R$ of the SMTI/SC structure. The coherence
length in the TI/SC is independent of the superconducting order parameter,
which is contrary to that in the SMTI/SC structure or a p-wave topological
superconductor. The distinct behaviors of the coherence lengths in two cases
reveal that the microscopic origins and the topological nature of the vortex
Majorana zero modes are different.

###### Acknowledgements.

This work was supported by the National Key R&D Program of China under Grant
No. 2019YFA0308603 and the Research Grants Council, University Grants
Committee, Hong Kong under Grant Nos. C7012-21GF and 17301220.

## Appendix

### A. The tight-binding model for numerical simulations

In this section, we give the explicit model for the tight-binding
calculations. The microscopic model for the bulk model of the hybrid system
with a single vortex is defined by the Hamiltonian,

$\hat{H}_{\mathrm{tot}}=\hat{H}_{\mathrm{TI}}+\hat{H}_{\mathrm{SC}}+\hat{H}_{\mathrm{Z}}.$

The $\hat{H}_{\mathrm{TI}}$ term describes the topological electronic
structure of the bulk system. In the basis of
$|\mathrm{P}1_{z}^{+},\uparrow\rangle$,$|\mathrm{P}1_{z}^{+},\downarrow\rangle$
$|\mathrm{P}2_{z}^{-},\uparrow\rangle$,
$|\mathrm{P}2_{z}^{-},\downarrow\rangle$, it can be written as (ZhangH-09np, ;
Shen-12book, ; Qi-11rmp, )

$\displaystyle\hat{H}_{\mathrm{TI}}$
$\displaystyle=\sum_{\mathbf{r}}\psi_{\mathbf{r}}^{\dagger}\left[\left(M-\sum_{\boldsymbol{\delta}}B_{\boldsymbol{\delta}}\right)\rho_{z}\sigma_{0}-\mu\rho_{0}\sigma_{0}\right]\psi_{\mathbf{r}}$
$\displaystyle+\sum_{\mathbf{r},\boldsymbol{\delta}}\left[\psi_{\mathbf{r}}^{\dagger}\frac{1}{2}\left(B_{\boldsymbol{\delta}}\rho_{z}\sigma_{0}-\frac{iv_{\boldsymbol{\delta}}}{a}\rho_{x}\sigma_{\delta}\right)\psi_{\mathbf{r}+\boldsymbol{\delta}}+h.c.\right]$

where $\rho_{i}$ and $\sigma_{i}$ ($i=0,x,y,z$) are Pauli matrices act on the
orbit and spin spaces, respectively.
$\psi_{\mathbf{r}}=[c_{\mathbf{r}1\uparrow},c_{\mathbf{r}1\downarrow},c_{\mathbf{r}2\uparrow},c_{\mathbf{r}2\downarrow}]^{T}$
are annihilation operators of the four-component spinor at position
$\mathbf{r}$. $2M$ is the band gap at $\Gamma$ point.
$\frac{\hbar}{a}v_{\boldsymbol{\delta}}$ and $B_{\boldsymbol{\delta}}$
describes the spin-dependent and spin-independent hoppings on the cubic
lattice along $\boldsymbol{\delta}$ direction with
$\boldsymbol{\delta}=x,y,z$. $a$ is the lattice constant. $\mu$ is the
chemical potential. For simplicity, we take $B_{x}=B_{y}=B_{z}=B$ then the
above model describes a strong topological insulator phase with gapless
surface states for $2>M/B>0$.

The proximity-induced superconductivity $\hat{H}_{\mathrm{SC}}$ can be
described by,

$\hat{H}_{\mathrm{SC}}=\sum_{\mathbf{r}}\left[\psi_{\mathbf{r}}^{\dagger}\Delta(\mathbf{r})e^{-i\phi(\mathbf{r}_{\parallel})}i\sigma_{y}\psi_{\mathbf{r}}^{\dagger}+h.c.\right]$

where
$\phi(\mathbf{r}_{\parallel})=\mathrm{Arg}(\mathbf{r}_{\parallel}-\mathbf{r}_{\parallel}^{0})$
with $\mathbf{r}_{\parallel}=(x,y)$ as the planar position vector,
$\mathbf{r}_{\parallel}^{0}$ as the coordinate for the vortex line and
$\mathrm{Arg}$ representing the the argument of the vector. The pairing
function is written as
$\Delta(\mathbf{r})=\Delta_{0}f(z)\tanh(r_{\parallel}/\xi)$ where $\Delta_{0}$
is the amplitude of the pairing function, $f(z)$ and
$\tanh(r_{\parallel}/\xi)$ are the distribution function along $z$ and the
planar direction respectively, and $\xi=v/\Delta_{0}$ is the coherent length
of the superconductor. Here we think the thickness of topological insulator
film is much less than the coherent length of the superconductor such that
$f(z)\simeq 1$.

$\hat{H}_{Z}$ describes the Zeeman term which is modeled as

$\hat{H}_{\mathrm{Z}}=\sum_{\mathbf{r}}V_{z}(z)\psi_{\mathbf{r}}^{\dagger}\rho_{0}\sigma_{z}\psi_{\mathbf{r}}$

with $V_{z}$ as the amplitude of the exchange field. The exchange field is
only restricted to several layers near the top surface of the system,
$V_{z}(z)=V_{z}\Theta(z_{m}-z)$ where $\Theta$ is Heaviside step function and
$z_{m}$ is the thickness of magnetic layers.

### B. Derivation of the effective vortex Hamiltonian

The inclusion of the s-wave superconductivity in three dimensional topological
insulator leads to the Bogoliubov–de Gennes (BdG) Hamiltonian as

$H_{\mathbf{k}}^{\mathrm{BdG}}=\left(\begin{array}[]{cc}H_{\mathbf{k}}&\Delta\\\
\Delta&-\sigma_{y}H_{-\mathbf{k}}^{*}\sigma_{y}\end{array}\right),$

where $\Delta$ is the superconducting pairing. In the presence of a vortex in
$z$ direction, the superconducting pairing $\Delta\to\Delta(r)e^{i\theta}$
with respect to the vortex core. We have used cylindrical coordinates
$(x,y,z)=(r\cos\theta,r\sin\theta,z)$. In view of the rotational symmetry
about $z$ axis, we can assign the quantum numbers $(k_{z},l,n)$ for the bulk
of the system, which are the momentum in $z$-direction, the angular momentum
in $xy$ plane and the radial quantum number respectively. These emergent
symmetries of the effective Hamiltonian allow us to obtain the radial BdG
equation
$H_{k_{z},l}^{\mathrm{BdG}}\Psi_{k_{z},l,n}=E_{k_{z},l,n}\Psi_{k_{z},l,n}$ at
a given $k_{z}$ and $l$ with

$H_{k_{z},l}^{\mathrm{BdG}}=\left(\begin{array}[]{cc}H_{k_{z},l}&\Delta(r)\\\
\Delta(r)&-\sigma_{y}H_{-k_{z},-l}^{*}\sigma_{y}\end{array}\right)$

where

$H_{k_{z},l}=\left(\begin{array}[]{cccc}M_{k_{z}}^{1-l}-\mu&0&vk_{z}&-viD_{r}^{l}\\\
0&M_{k_{z}}^{l}-\mu&-viD_{r}^{1-l}&-vk_{z}\\\
vk_{z}&-viD_{r}^{l}&-M_{k_{z}}^{1-l}-\mu&0\\\
-viD_{r}^{1-l}&-vk_{z}&0&-M_{k_{z}}^{l}-\mu\end{array}\right).$

Here we have introduced $D_{r}^{l}=\partial_{r}+l/r$ and
$M_{k_{z}}^{l}=m-Bk_{z}^{2}+BD_{r}^{1-l}D_{r}^{l}$. To find the effective
Hamiltonian for the vortex line, we first rewrite the Hamiltonian as
$H(k_{z},l)=H_{1}(\partial_{r},l)+H_{2}(k_{z},l)$, then solve solutions for
$k_{z}$-independent part $H_{1}$and finally project $k_{z}$-dependent part
$H_{2}$ onto the relevant bands. $H_{1}$ can be expressed as a direct sum of
$+$ and $-$ sectors
$H_{1}^{\pm}(\partial_{r},l)=-iv(\partial_{r}+\frac{1}{2r})\tau_{z}\nu_{x}-\mu\tau_{z}+\Delta(r)\tau_{x}\mp\frac{v}{2r}\tau_{0}\nu_{y}\pm\frac{vl}{r}\tau_{z}\nu_{y}$.
For $l=0$, there exists a chiral symmetry
$\tau_{y}\nu_{x}H_{1}^{\pm}(\partial_{r},l=0)\tau_{y}\nu_{x}=-H_{1}^{\pm}(\partial_{r},l=0)$
which plays an essential role in the determination of the zero energy
solutions. The zero energy solutions for $H_{1}^{\pm}$ can be solved as

$\displaystyle|\phi_{+}\rangle$
$\displaystyle=Ne^{-\int_{0}^{r}dr^{\prime}\frac{\Delta(r^{\prime})}{v}}\left(\begin{array}[]{c}-J_{1}(k_{F}r)\\\
0\\\ 0\\\ iJ_{0}(k_{F}r)\\\ J_{0}(k_{F}r)\\\ 0\\\ 0\\\
iJ_{1}(k_{F}r)\end{array}\right),$ $\displaystyle|\phi_{-}\rangle$
$\displaystyle=Ne^{-\int_{0}^{r}dr^{\prime}\frac{\Delta(r^{\prime})}{v}}\left(\begin{array}[]{c}0\\\
J_{0}(k_{F}r)\\\ iJ_{1}(k_{F}r)\\\ 0\\\ 0\\\ J_{1}(k_{F}r)\\\
-iJ_{0}(k_{F}r)\\\ 0\end{array}\right)$

with the normalization factor defined as

$4\pi
N^{2}\int_{0}^{\infty}drre^{-2\int_{0}^{r}dr^{\prime}\frac{\Delta(r^{\prime})}{v}}[J_{0}^{2}(k_{F}r)+J_{1}^{2}(k_{F}r)]=1.$

The expectation values of the remaining terms in $H_{2}$ with respect to the
zero energy solutions for $H_{1}^{\pm}$ can be calculated as

$\displaystyle\langle\phi_{s}|\mathcal{M}_{k_{z}}^{0}|\phi_{s^{\prime}}\rangle$
$\displaystyle=s\delta_{ss^{\prime}}\mathcal{F}(k_{F}^{2}\xi^{2})\left[m-B(k_{z}^{2}+k_{F}^{2}+\frac{1}{\xi^{2}})\right],$
$\displaystyle\langle\phi_{s}|vk_{z}\tau_{z}\rho_{x}\sigma_{z}|\phi_{s^{\prime}}\rangle$
$\displaystyle=-s\delta_{s,-s^{\prime}}ivk_{z}\mathcal{F}(k_{F}^{2}\xi^{2})$

where $s,s^{\prime}=\pm$ and
$\mathcal{M}_{k_{z}}^{0}=\rho_{z}\otimes\mathrm{diag}(M_{k_{z}}^{1},M_{k_{z}}^{0},-M_{k_{z}}^{0},-M_{k_{z}}^{1})$,
which leads to Eq. (1) in main text.

### C. Zero energy solutions for SMTI/SC heterostructure

In this section, we calculate the zero energy solutions for SMTI/SC
heterostructure and the energy splitting due to their overlapping. In polar
coordinates, by using the Ansatz for the wavefunction
$\Psi_{l,n}(r,\theta)=e^{i[l-\frac{1}{2}(\tau_{z}+\sigma_{z})]\theta}\Psi_{l,n}(r)$
with $l$ being an integer to make it single-valuedness, we can separate
angular and radial variables. The radial Hamiltonian is rewritten as
$H_{l}(r)\Psi_{l,n}(r)=E_{l,n}\Psi_{l,n}(r)$ with

$H_{l}(r)=\left(\begin{array}[]{cc}h_{l}(r)&\Delta(r)\\\
\Delta^{*}(r)&-\sigma_{y}h_{-l}^{*}(r)\sigma_{y}\end{array}\right)$ (6)

where $h_{l}(r)$ is given by

$h_{l}(r)=\left(\begin{array}[]{cc}-\mu+M(r)&-ivD_{r}^{l}\\\
-ivD_{r}^{1-l}&-\mu-M(r)\end{array}\right).$ (7)

Due to the presence of the particle-hole symmetry $H_{l}=-\Xi
H_{-l}\Xi^{-1}$with $\Xi=\tau_{y}\sigma_{y}\mathcal{K}$, if $E$ is an
eigenvalue with the eigenfunction
$\Psi_{l,E}(r)=[u_{\uparrow},u_{\downarrow},v_{\downarrow},v_{\uparrow}]^{T}$,
then $-E$ is also an eigenvalue and the corresponding eigenfunction is
$\Psi_{-l,-E}(r)=[-v_{\uparrow}^{*},v_{\downarrow}^{*},u_{\downarrow}^{*},-u_{\uparrow}^{*}]^{T}$.
The zero energy state only exists for $l=0$ and needs to be an eigenstate of
$\Xi$, i.e. $\Xi\Psi_{0,0}(r)=\zeta\Psi_{0,0}(r)$, which gives constraints on
components of the eigenfunction $u_{\uparrow}=-\zeta v_{\uparrow}^{*}$ and
$u_{\downarrow}=\zeta v_{\downarrow}^{*}$.By redefining the spinor
$u_{\uparrow}=\widetilde{u}_{\uparrow}e^{-i\pi/4-\zeta\int_{0}^{r}dr^{\prime}\Delta(r^{\prime})/\hbar
v}$ and
$u_{\downarrow}=\widetilde{u}_{\downarrow}e^{i\pi/4-\zeta\int_{0}^{r}dr^{\prime}\Delta(r^{\prime})/\hbar
v}$ , the four coupled differential equations in radial BdG equations is
reduced to two real equations

$\displaystyle-[\mu-M(r)]\widetilde{u}_{\uparrow}+v\partial_{r}\widetilde{u}_{\downarrow}$
$\displaystyle=0,$ $\displaystyle
v\left(\partial_{r}+\frac{1}{r}\right)\widetilde{u}_{\uparrow}+[\mu+M(r)]\widetilde{u}_{\downarrow}$
$\displaystyle=0.$ (8)

For $0<\mu<V_{z}$, by solution Eq. (8), we find a solution localized at the
vortex core,

$\left(\begin{array}[]{c}u_{\uparrow}\\\
u_{\downarrow}\end{array}\right)_{\mathrm{core}}\sim
e^{-i\pi/4\sigma_{z}-\int_{0}^{r}dr^{\prime}\frac{\Delta(r^{\prime})}{\hbar
v}}\left(\begin{array}[]{c}J_{1}(k_{F}r)\\\ -J_{0}(k_{F}r)\end{array}\right),$

with the positive eigenvalue $\zeta=+1$ of $\Xi$ and a bound state localized
at the interface

$\displaystyle\left(\begin{array}[]{c}u_{\uparrow}\\\
u_{\downarrow}\end{array}\right)_{\mathrm{inter}}$
$\displaystyle\sim\frac{e^{-i\pi/4\sigma_{z}}}{\sqrt{r}}e^{\int_{R}^{r}dr^{\prime}\frac{\Delta(r^{\prime})}{\hbar
v}}$ (11)
$\displaystyle\times\begin{cases}\left(\begin{array}[]{c}\sqrt{\frac{\mu+V_{z}}{2V_{z}}}\\\
\sqrt{\frac{V_{z}-\mu}{2V_{z}}}\end{array}\right)e^{-\frac{\sqrt{m^{2}-\mu^{2}}(r-R)}{\hbar
v}},&r>R,\\\
\left(\begin{array}[]{c}-\sin\left[k_{F}(r-R)-\delta+\frac{\pi}{4}\right]\\\
\cos\left[k_{F}(r-R)-\delta+\frac{\pi}{4}\right]\end{array}\right),&r<R,\end{cases}$
(12)

with the negative eigenvalue $\zeta=-1$ of $\Xi$ where
$\delta=\arctan\sqrt{\frac{V_{z}+\mu}{V_{z}-\mu}}+\frac{\pi}{4}$ is a phase
determined by matching the wavefunction at the interface. To estimate the
energy splitting for Majorana modes, we multiply
$\langle\psi_{\mathrm{core}}|$ to the BdG equation
$H_{0}(r)|\Psi_{+}\rangle=E_{+}|\Psi_{+}\rangle$ which yields
$E_{+}=\frac{\langle\psi_{\mathrm{core}}|H_{0}(r)|\Psi^{+}\rangle}{\langle\psi_{\mathrm{core}}|\Psi^{+}\rangle}$.
Then using the relation
$H_{0}(r)|\psi_{\mathrm{inter}}\rangle=-iv\tau_{z}\sigma_{x}|\psi_{\mathrm{inter}}\rangle\delta(r-R)$
and $H_{0}(r)|\psi_{\mathrm{core}}\rangle=0$, we arrive the expression for the
energy splitting for the zero mode in the main text [Eq. (4)].

## References

  * (1) F. Wilczek, Majorana returns, Nat. Phys. 5, 614 (2009).
  * (2) A. Stern, Non-Abelian states of matter, Nature 464, 187 (2010).
  * (3) X. L. Qi and S. C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * (4) J. Alicea, New directions in the pursuit of Majorana fermions in solid state systems, Rep. Prog. Phys. 75, 076501 (2012).
  * (5) S. Q. Shen, Topological insulators, Springer Series of Solid State Science, Vol. 174 (Springer, Heidelberg, 2012).
  * (6) C. W. J. Beenakker, Search for Majorana fermions in superconductors, Annu. Rev. Condens. Matter Phys. 4, 113 (2013).
  * (7) M. Sato and S. Fujimoto, Majorana fermions and topology in superconductors, J. Phys. Soc. Jpn 85, 072001 (2016).
  * (8) R. M. Lutchyn, E. P. A. M. Bakkers, L. P. Kouwenhoven, P. Krogstrup, C. M. Marcus, and Y. Oreg, Majorana zero modes in superconductor–semiconductor heterostructures, Nat. Rev. Mater. 3, 52 (2018).
  * (9) A. Y. Kitaev, Fault-tolerant quantum computation by anyons, Ann. Phys. 303, 2 (2003).
  * (10) M. H. Freedman, A. Kitaev, M. J. Larsen and Z. Wang, Topological quantum computation, Bull. Amer. Math. Soc. 40, 31 (2003)
  * (11) C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. D. Sarma, Non-Abelian anyons and topological quantum computation, Rev. Mod. Phys. 80, 1083–1159 (2008).
  * (12) L. Fu and C. L. Kane, Superconducting Proximity Effect and Majorana Fermions at the Surface of a Topological Insulator, Phys. Rev. Lett. 100, 096407 (2008).
  * (13) R. M. Lutchyn, J. D. Sau, and S. D. Sarma, Majorana fermions and a topological phase transition in semiconductor-superconductor heterostructures, Phys. Rev. Lett. 105, 077001 (2010).
  * (14) Y. Oreg, G. Refael, and F. von Oppen, Helical Liquids and Majorana Bound States in Quantum Wires, Phys. Rev. Lett. 105, 177002 (2010).
  * (15) J. D. Sau, R. M. Lutchyn, S. Tewari, and S. Das Sarma, Generic New Platform for Topological Quantum Computation Using Semiconductor Heterostructures, Phys. Rev. Lett. 104, 040502 (2010).
  * (16) J. D. Sau, R. M. Lutchyn, S. Tewari, and S. Das Sarma, Robustness of Majorana fermions in proximity-induced superconductors, Phys. Rev. B 82, 094522 (2010).
  * (17) V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Signatures of Majorana fermions in hybrid superconductor-semiconductor nanowire devices, Science 336, 1003 (2012).
  * (18) A. Das, Y. Ronen, Y. Most, Y. Oreg, M. Heiblum, and H. Shtrikman, Zero-bias peaks and splitting in an Al–InAs nanowire topological superconductor as a signature of Majorana fermions, Nat. Phys. 8, 887 (2012).
  * (19) S. Nadj-Perge, I. K. Drozdov, J. Li, H. Chen, S. Jeon, J. Seo, A. H. MacDonald, B. A. Bernevig, and A. Yazdani, Observation of Majorana fermions in ferromagnetic atomic chains on a superconductor, Science 346, 602 (2014).
  * (20) J. P. Xu, M. X. Wang, Z. L. Liu, J. F. Ge, X. J. Yang, C. H. Liu, Z. A. Xu, D. D. Guan, C. L. Gao, D. Qian, Y. Liu, Q. H. Wang, F. C. Zhang, Q. K. Xue, and J. F. Jia, Experimental Detection of a Majorana Mode in the core of a Magnetic Vortex inside a Topological Insulator-Superconductor $\mathrm{Bi_{2}Te_{3}/NbSe_{2}}$ Heterostructure, Phys. Rev. Lett. 114, 017001 (2015).
  * (21) H. -H. Sun, K. -W. Zhang, L. -H. Hu, C. Li, G. -Y. Wang, H. -Y. Ma, Z. -A. Xu, C. -L. Gao, D. -D. Guan, Y. -Y. Li, C. Liu, D. Qian, Y. Zhou, L. Fu, S. -C. Li, F. -C. Zhang, and J. -F. Jia, Majorana Zero Mode Detected with Spin Selective Andreev Reflection in the Vortex of a Topological Superconductor, Phys. Rev. Lett. 116, 257003 (2016).
  * (22) L.-H. Hu, C. Li, D. -H. Xu, Y. Zhou, and F. -C. Zhang, Theory of spin-selective Andreev reflection in the vortex core of a topological superconductor, Phys. Rev. B 94, 224501 (2016).
  * (23) P. Zhang, K. Yaji, T. Hashimoto, Y. Ota, T. Kondo, K. Okazaki, Z. Wang, J. Wen, G. D. Gu, H. Ding, and S. Shin, Observation of topological superconductivity on the surface of an iron-based superconductor, Science 360, 182–186 (2018).
  * (24) D. Wang, L. Kong, P. Fan, H. Chen, S. Zhu, W. Liu, L. Cao, Y. Sun, S. Du, J. Schneeloch, R. Zhong, G. Gu, L. Fu, H. Ding, and H. J. Gao, Evidence for Majorana bound states in an iron-based superconductor, Science 362, 333–335 (2018).
  * (25) T. Machida, Y. Sun, S. Pyon, S. Takeda, Y. Kohsaka, T. Hanaguri, T. Sasagawa, and T. Tamegai, Zero-energy vortex bound state in the superconducting topological surface state of Fe(Se,Te), Nat. Mater. 18, 811–815 (2019).
  * (26) L. Y. Kong, L. Cao, S. Y. Zhu, M. Papaj, G. Y. Dai, G. Li, P. Fan, W. Y. Liu, F. Z. Yang, X. C. Wang, S. X. Du, C. Q. Jin, L. Fu, H.-J. Gao, and H. Ding, Tunable Vortex Majorana Zero Modes in LiFeAs Superconductor, Nat. Commun. 12, 4146 (2021).
  * (27) Q. Liu, C. Chen, T. Zhang, R. Peng, Y.-J. Yan, C.-H.-P. Wen, X. Lou, Y.-L. Huang, J.-P. Tian, X.-L. Dong, G.-W.Wang, W.-C. Bao, Q.-H. Wang, Z.-P. Yin, Z.-X. Zhao, and D.-L. Feng, Robust and Clean Majorana Zero Mode in the Vortex Core of High-Temperature Superconductor $\mathrm{(Li_{0.84}Fe_{0.16})OHFeSe}$, Phys. Rev. X 8, 041056 (2018).
  * (28) W. Liu, L. Cao, S. Zhu, L. Kong, G. Wang, M. Papaj, P. Zhang, Y.-B. Liu, H. Chen, G. Li, F. Yang, T. Kondo, S. Du, G.-H. Cao, S. Shin, L. Fu, Z. Yin, H.-J. Gao, and H. Ding, A new Majorana platform in an Fe–As bilayer superconductor, Nat. Commun. 11, 5688 (2020).
  * (29) S. Zhu, L. Kong, L. Cao, H. Chen, M. Papaj, S. Du, Y. Xing, W. Liu, D. Wang, C. Shen, F. Yang, J. Schneeloch, R. Zhong, G. D. Gu, L. Fu, Y. Zhang, H. Ding, and H.-J. Gao, Nearly quantized conductance plateau of vortex zero mode in an iron-based superconductor, Science 367, 189–192 (2020).
  * (30) M. Li, G. Li, L. Cao, X. Zhou, X. Wang, C. Jin, C. K. Chiu, S. J. Pennycook, Z. Wang, and Hong-Jun Gao, Ordered and tunable Majorana-zero-mode lattice in naturally strained LiFeAs, Nature 606, 890–895 (2022).
  * (31) C. -K. Chiu, T. Machida, Y. Huang, T Hanaguri, and F. -C. Zhang, Scalable majorana vortex modes in iron-based superconductors, Sci. Adv. 6, eaay0443 (2020).
  * (32) E. J. K$\ddot{\mathrm{o}}$nig and P. Coleman, Crystalline symmetry protected helical Majorana modes in the iron pnictides, Phys. Rev. Lett. 122, 207001 (2019).
  * (33) R. -X. Zhang and S. Das Sarma, Intrinsic Time-Reversal-Invariant Topological Superconductivity in Thin Films of Iron-Based Superconductors, Phys. Rev. Lett. 126, 137001 (2021).
  * (34) M. Kheirkhah, Z. Yan, and F. Marsiglio, Vortex-line topology in iron-based superconductors with and without second-order topology, Phys. Rev. B 103, L140502 (2021).
  * (35) L. -H. Hu, X. Wu, C. -X. Liu, and R. -X. Zhang, Competing Vortex Topologies in Iron-based Superconductors, Phys. Rev. Lett. 129, 277001 (2022)
  * (36) D. A. Ivanov, Non-Abelian Statistics of Half-Quantum Vortices in p -Wave Superconductors, Phys. Rev. Lett. 86, 268 (2001).
  * (37) N. Read and D. Green, Paired states of fermions in two dimensions with breaking of parity and time-reversal symmetries and the fractional quantum Hall effect, Phys. Rev. B 61, 10267 (2000).
  * (38) B. Zhou, H. Z. Lu, R. L. Chu, S. Q. Shen, and Q. Niu, Finite size effects on helical edge states in a quantum spin-Hall system, Phys. Rev. Lett. 101, 246807 (2008).
  * (39) A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Classification of Topological Insulators and Superconductors, Phys. Rev. B 78, 195125 (2008).
  * (40) A. Kitaev, Periodic table for topological insulators and superconductors, AIP Conf. Proc. 1134, 22 (2009).
  * (41) C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Classification of topological quantum matter with symmetries, Rev. Mod. Phys. 88, 035005 (2016).
  * (42) A. Haim and Y. Oreg, time-reversal-invariant topological superconductivity in one and two dimensions, Physics Reports 825, 1 (2019).
  * (43) L. Fu and E. Berg, Odd-parity topological superconductors: theory and application to $\mathrm{Cu_{x}Bi_{2}Se_{3}}$, Phys. Rev. Lett. 105, 097001 (2010).
  * (44) F. Zhang, C. L. Kane, and E. J. Mele, Time Reversal Invariant Topological Superconductivity and Majorana Kramers Pairs, Phys. Rev. Lett. 111, 056402 (2013).
  * (45) X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Topological invariants for the Fermi surface of a time-reversal- invariant superconductor, Phys. Rev. B 81, 134508 (2010).
  * (46) M. Sato, Topological odd-parity superconductors, Phys. Rev. B 81, 220504(R) (2010).
  * (47) M. Atiyah, and I. Singer, The Index of Elliptic Operators on Compact Manifolds, Bull. Amer. Math. Soc. 69, 322 (1963).
  * (48) R. Jackiw, and P. Rossi, Zero modes of the vortex-fermion system, Nucl. Phys. B 190, 681 (1981).
  * (49) R. Jackiw, and S.-Y. Pi, Chiral Gauge Theory for Graphene, Phys. Rev. Lett. 98, 266402 (2007).
  * (50) P. Hosur, P. Ghaemi, R. S. K. Mong, and A. Vishwanath, Majorana Modes at the Ends of Superconductor Vortices in Doped Topological Insulators, Phys. Rev. Lett. 107, 097001 (2011).
  * (51) M. Wimmer, Algorithm 923: Efficient numerical computation of the pfaffian for dense and banded skew-symmetric matrices, ACM Trans. Math. Softw. 38, 30 (2012).
  * (52) C. -K. Chiu, M. J. Gilbert, and T. L. Hughes, Vortex lines in topological insulator-superconductor heterostructures, Phys. Rev. B 84, 144507 (2011).
  * (53) H. Zhang, C. X. Liu, X. L. Qi, X. Dai, Z. Fang, and S. C. Zhang. Topological insulators in $\mathrm{Bi}_{2}\mathrm{Se}_{3}$, $\mathrm{Bi}_{2}\mathrm{Te}_{3}$ and $\mathrm{Sb}_{2}\mathrm{Te}_{3}$ with a single Dirac cone on the surface, Nat. Phys. 5, 438 (2009).
  * (54) W. P. Su, J. R. Schrieffer and A. J. Heeger, Soliton excitations in polyacetylene, Phys. Rev. B 22, 2099 (1980).
  * (55) A. Kitaev, Anyons in an exactly solved model and beyond, Ann. Phys. 321, 2 (2006).
  * (56) J. Linder, T. Yokoyama, and A. Sudbø, Anomalous finite size effects on surface states in the topological insulator $\mathrm{Bi}_{2}\mathrm{Se}_{3}$, Phys. Rev. B 80, 205401( 2009).
  * (57) E. J. Weinberg, Index calculations for the fermion-vortex system, Phys. Rev. D 24, 2669 (1981).
  * (58) T. Fukui and T. Fujiwara, Topological Stability of Majorana Zero Modes in Superconductor–Topological Insulator Systems, J. Phys. Soc. Jpn. 79, 033701 (2010).
  * (59) B. Roy, and P. Goswami, $Z_{2}$ index for gapless fermionic modes in the vortex core of three dimensional paired Dirac fermions, Phys. Rev. B 89 144507 (2014).
  * (60) J. C. Y. Teo and C. L. Kane, Majorana Fermions and Non-Abelian Statistics in Three Dimensions, Phys. Rev. Lett. 104, 046401 (2010).
  * (61) X. -L. Qi, T. L. Hughes, and S. -C. Zhang, Chiral topological superconductor from the quantum Hall state, Phys. Rev. B 82, 184516 (2010).
  * (62) A. R. Akhmerov, J. Nilsson, C. W. J. Beenakker, Electrically detected interferometry of Majorana fermions in a topological insulator, Phys. Rev. Lett. 102, 216404(2009).
  * (63) L. Fu and C. L. Kane, Probing Neutral Majorana Fermion Edge Modes with Charge Transport, Phys. Rev. Lett. 102, 216403 (2009)
  * (64) J. C. Y. Teo, and C. L.Kane, Topological Defects and Gapless Modes in Insulators and Superconductors, Phys. Rev. B 82, 115120 (2010).
  * (65) X. -L. Qi, E. Witten, and S. -C. Zhang, Axion topological field theory of topological superconductors, Phys. Rev. B 87, 134519 (2013).
  * (66) J. Wang, Q. Zhou, B. Lian, and S.-C. Zhang, Chiral topological superconductor and half-integer conductance plateau from quantum anomalous Hall plateau transition, Phys. Rev. B 92, 064520 (2015).
  * (67) J. Wang, Electrically tunable topological superconductivity and Majorana fermions in two dimensions, Phys. Rev. B 94, 214502(2016).
  * (68) B. Lian, X. -Q. Sun, A. Vaezi, X. -L. Qi, S. -C. Zhang, Topological Quantum Computation Based on Chiral Majorana Fermions, Proc. Nat. Acad. Sci. 115, 10938 (2018).
  * (69) Q. Yan, H. Li, J. Zeng, Q. -F. Sun, and X. C. Xie, A Majorana perspective on understanding and identifying axion insulators, Commun. Phys.4, 239(2021).
  * (70) K. T. Law, P. A. Lee, and T. K. Ng, Majorana Fermion Induced Resonant Andreev Reflection, Phys. Rev. Lett. 103. 237001 (2009).
  * (71) M. Cheng, R. M. Lutchyn, V. Galitski, and S. Das Sarma, Splitting of Majorana modes due to intervortex tunneling in a p + ip superconductor, Phys. Rev. Lett. 103, 107001 (2009).
  * (72) M. Cheng, R. M. Lutchyn, V. Galitski, and S. Das Sarma, Tunneling of anyonic Majorana excitations in topological superconductors, Phys. Rev. B 82, 094504 (2010).
  * (73) M. Mogi, Y. Okamura, M. Kawamura, R. Yoshimi, K. Yasuda, A. Tsukazaki, K. S. Takahashi, T. Morimoto, N. Nagaosa, M. Kawasaki, Y. Takahashi, and Y. Tokura, Experimental signature of parity anomaly in semi-magnetic topological insulator, Nat. Phys. 18, 390 (2022).
  * (74) J. Y. Zou, R. Chen, B. Fu, H. W. Wang, Z. A. Hu, and S. Q. Shen, Half-Quantized Hall Effect at the Parity-Invariant Fermi Surface, Phys. Rev. B 107, 125153 (2023).


# Safe Control using Occupancy Grid Map-based Control Barrier Function (OGM-
CBF)

Golnaz Raja, Teemu Mökkönen, Reza Ghabcheloo  This work was jointly supported
by Academy of Finland (SA) 345517, under SA-NSF joint call on artificial
intelligence and wireless communication, and Together for RISc-V Technology
and ApplicatioNs (TRISTAN) under Grant 101095947. _(Corresponding author:
Golnaz Raja.)_ All authors are with the Faculty of Engineering and Natural
Sciences, Tampere University, 33720 Tampere, Finland. (e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>reza.ghabcheloo@tuni.fi).

###### Abstract

Safe navigation in unknown environments stands as a significant challenge in
the field of robotics. Control Barrier Function (CBF) is a strong mathematical
tool to guarantee safety requirements. However, a common assumption in many
works is that the CBF is already known and obstacles have predefined shapes.
In this letter, we present a novel method called Occupancy Grid Map-based
Control Barrier Function (OGM-CBF), which defines Control Barrier Function
based on Occupancy Grid Maps. This enables generalization to unknown
environments while generating online local or global maps of the environment
using onboard perception sensors such as LiDAR or camera. With this method,
the system guarantees safety via a single, continuously differentiable CBF per
time step, which can be represented as one constraint in the CBF-QP
optimization formulation while having an arbitrary number of obstacles with
unknown shapes in the environment. This enables practical real-time
implementation of CBF in both unknown and known environments. The efficacy of
OGM-CBF is demonstrated in the safe control of an autonomous car in the CARLA
simulator and a real-world industrial mobile robot.

###### Index Terms:

Autonomous Vehicle Navigation, Collision Avoidance, Vision-Based Navigation,
Sensor-based Control.

##  I INTRODUCTION

Applications of mobile robots in domains such as autonomous vehicles,
industrial robots, and multi-agent systems are often tasked to respect safety
while navigating in complex and unknown environments, which stands as one of
the most formidable challenges in the field of robotics [1]. One promising
approach to integrate safety into design is via Control Barrier Function
(CBF)[2] due to the strong mathematical guarantee of satisfying safety
requirements. CBFs have been applied in various fields such as legged robots
[3], autonomous vehicles [4], multi-robot systems [5] [6], decentralized
multi-agent control [7], and drones [8].

In recent years, CBF has been successfully integrated with Control Lyapunov
Function (CLF) via Quadratic Programs (QP) for safety critical systems [9],
enabling the simultaneous satisfaction of safety and stability. Moreover,
given a potentially unsafe nominal control policy, CBF can be embedded as a
QP-based safety filter to modify the unsafe control command in a minimally
invasive manner to prevent violations of the safety constraints[2].

In designing CBF, it is generally assumed that unsafe sets are already known
(usually as circular shapes), which is not applicable to unsafe sets with
arbitrary shapes in uncharted environments. Practical implementation of CBF in
unknown environments relies on a sufficient awareness of the environment,
unsafe regions, and their changes to reach a quick, yet smooth generation of a
safety certificate as a response. This requires online estimation of unsafe
sets and their integration into the CBF formulation.

One practical approach to operate safely in unknown and dynamic environments
without sufficient prior information about the environment is the integration
of sensory data into the CBF formulation. This enables the online synthesis of
CBF based on the aggregated sensor data from perceptive sensors such as LiDAR
or camera.

To construct CBF based on sensor data, one approach is to learn both the
classification of safe and unsafe regions and CBF simultaneously based on
sensor measurements. In [10], the authors used neural networks to learn the
CBF directly in the observation space using local observations from LiDAR,
without the need for the robot’s state estimation and localization. In [11],
the authors designed a neural network, called BarrierNet, to learn CBF with
respect to the user-defined safe sets, and in [12], a neural network was
trained to extract features from front-view images and learn the BarrierNet as
an end-to-end vision-based solution for lane keeping and obstacle avoidance.
In [13], CBF is synthesized based on the online approximation of the signed
distance function from LiDAR measurements using an incremental training
approach with replay memory, and in [14], CBF is constructed based on the
classification of safe and unsafe sets using a support vector machine on LiDAR
measurements. In these approaches, sensor data at each time step is
immediately employed to construct CBF, resulting in a lack of historical
memory for obstacles that have been previously seen. Moreover, these methods
are bound to specific sensor modalities and cannot combine information from
different sensors for CBF generation.

One other approach is to distinguish the tasks of classifying safe and unsafe
areas from the synthesis of CBF, treating them as separate modules.
Consequently, established classification techniques, such as semantic
segmentation and mapping, can be utilized to identify safe and unsafe areas,
upon which CBF can be defined. This modular approach can enhance the efficacy
of the final controller, particularly in the case of significant changes in
the type of environment where the safe controller is deployed. In V-CBF [15],
Control Barrier Functions are defined based on segmented unsafe sets using
RGBD cameras. Although this framework successfully navigates safely and avoids
obstacles with unknown shapes in uncharted environments, it relies solely on
instantaneous sensor inputs and lacks memory of previously encountered
obstacles. Moreover, as the robot is considered as a point, to cover the size
and shape of the robot, a handful of the most conservative points of the
constructed CBF should be selected and added as multiple constraints in QP
formulation for a successful safe navigation. It is not trivial how these
points are selected.

One practical way to approximate the shapes of obstacles is Occupancy Grid Map
(OGM). OGM, as a multi-dimensional discrete map, efficiently captures the
complexity of the surrounding environment, including the state and shape of
all detected obstacles. With accurate state estimation for the robot, the OGM
illustrates the relative positions of obstacles with respect to the robot,
providing comprehensive situational awareness. The map employs a probabilistic
formulation to store all previously seen obstacles, retaining information
about these past encounters. Additionally, obstacles within the OGM can be
inflated according to the robot’s size, enhancing safe navigation. OGM not
only embeds essential information for the robot’s safe navigation, including
the position and inflated shape of obstacles, but it also can be constructed
online using real-time sensor data. Online mapping facilitates the robot’s
navigation in unknown environments since mapping and synthesizing CBFs can be
done in real time. Moreover, various sensor configurations can be utilized to
generate OGMs, making it a versatile tool for representing unsafe zones.
However, OGM alone cannot be directly used for generating CBF since OGM
behaves in a binary manner and does not hold correct variational behavior near
the boundary of unsafe sets. Thus, a domain translation function is necessary
to transform the occupancy grid map into a representation suitable for CBF
synthesis. In [16], CBF was designed using elliptical outer approximation on
occupancy grid maps. However, the obstacles were approximated with minimum
bounding circles, which can be overly conservative. Additionally, this
approach requires a separate constraint for each obstacle.

In this paper, we introduce a novel method to define CBF based on occupancy
grid maps called Occupancy Grid Map-based Control Barrier Function (OGM-CBF).
Our contribution is as follows:

  * •

CBF is represented as a single constraint at each time step in a QP
formulation while assuming the presence of an arbitrary number of obstacles
with unknown shapes, enabling real-time and computationally light optimization
to guarantee safety.

  * •

The designed safe controller is sensor agnostic. In addition, CBF can be
constructed from hybrid sensor modalities using occupancy grid maps.

  * •

The constructed CBF has a memory of all previously seen obstacles.

  * •

The constructed CBF accounts for the size and shape of the robot using an
inflated occupancy grid map.

  * •

We evaluated the efficacy of the proposed method in the CARLA simulator in
complex scenarios as well as with an industrial mobile robot in real-world
experiments.

The outline of this paper is as follows. Section II reviews the principles of
Control Barrier Function, Occupancy Grid Map, and Signed Distance Function.
Section III describes our designed control methodology. Section IV explains
the designed smoothing function used in our CBF. Section V is dedicated to the
demonstration of the introduced control strategy in the CARLA simulator and on
a real robot. Finally, in section VI, the conclusions and future work are
provided.

##  II PRELIMINARIES AND BACKGROUND

Consider a non-linear control-affine system:

$\displaystyle\dot{\boldsymbol{x}}$
$\displaystyle=f(\boldsymbol{x})+g(\boldsymbol{x})\boldsymbol{u}$ (1)
$\displaystyle z$ $\displaystyle=O(\boldsymbol{x},\boldsymbol{q})$

where $\boldsymbol{x}\in\mathbb{X}\subseteq\mathbb{R}^{n}$,
$f:\mathbb{R}^{n}\to\mathbb{R}^{n}\quad\text{and}\quad
g:\mathbb{R}^{n}\to\mathbb{R}^{2}$ are locally Lipschitz continuous functions,
$O:\mathbb{R}^{n}\times\mathbb{R}^{q}\rightarrow\mathbb{R}^{k}$ is a
perception sensor module such as camera or LiDAR capable of generating
occupancy grid maps, $\boldsymbol{x}$ and $z$ are the state and observation
vectors, respectively, and $\boldsymbol{q}$ is environment variables. In this
paper, our system is a kinematic model of a robot with
$\boldsymbol{u}=[\mathrm{v},\boldsymbol{\omega}]^{T}$, the linear and angular
speeds as input.

######  Definition 1 (Local Lipschitz continuity)

Let $f:D\subseteq\mathbb{R}^{n}\to\mathbb{R}^{m}$ be a function. $f$ is
locally Lipschitz continuous at a point $x\in D$ if there exist constants
$M>0$ and $\delta>0$ such that for all $x^{\prime}\in D$ with $\lVert
x-x^{\prime}\rVert\leq\delta$, the inequality $\lVert
f(x)-f(x^{\prime})\rVert\leq M\lVert x-x^{\prime}\rVert$ is satisfied.

###  II-A Control Barrier Function

In this paper, we define the safe system as follows:

######  Definition 2 (Safe system)

The system (1) is categorized as safe if
$\boldsymbol{x}(t)\in\mathbb{X}_{s}\subseteq\mathbb{X}$,
$\boldsymbol{u}(t)\in\mathbb{U}_{s}\subseteq\mathbb{U}$, $\forall t\geq 0$.
Where $\mathbb{X}_{s}$ is the set of safe states, and $\mathbb{U}_{s}$ is the
set of admissible control inputs.

Consider $S\subseteq\mathbb{X}_{s}$ to be zero-superlevel set of a
continuously differentiable function
$h(\boldsymbol{x}):\mathbb{X}\rightarrow\mathbb{R}$, yielding:

$\displaystyle S=\left\\{\boldsymbol{x}\in\mathbb{X}\;|\;h(\boldsymbol{x})\geq
0\right\\}$ (2)

To keep the system (1) remain in the safe set for all time, the set $S$ should
be a forward invariant set.

######  Definition 3 (Forward invariant set)

A set $S$ is forward invariant set for the control system (1) if state
trajectories starting in $S$ will remain in $S$ for the future time. That is,
if $\boldsymbol{x}(0)\in S$, there exists
$\boldsymbol{u}(t)\in\mathbb{U}_{s}\subseteq\mathbb{U}$, such that
$\boldsymbol{x}(t)\in S$, $\forall t\geq 0$ [17].

As proven in [2], given the dynamic system (1), the function
$h(\boldsymbol{x})$ is a control barrier function if there exists an extended
class $\kappa$ function $\alpha(.)$ such that for all
$\boldsymbol{x}\in\mathbb{X}$:

$\sup_{\boldsymbol{u}\in\mathbb{U}_{s}}[\frac{\partial
h(\boldsymbol{x})}{\partial\boldsymbol{x}}(f(\boldsymbol{x})+g(\boldsymbol{x})\boldsymbol{u})]\geq-\alpha(h(\boldsymbol{x})).$
(3)

We will later use $\dot{h}(\boldsymbol{x},\boldsymbol{u})=\frac{\partial
h(\boldsymbol{x})}{\partial\boldsymbol{x}}(f(\boldsymbol{x})+g(\boldsymbol{x})\boldsymbol{u})$
notation for brevity.

######  Definition 4 (Extended class $\kappa$ function)

An extended class $\kappa$ function is a function
$\alpha:\mathbb{R}\rightarrow\mathbb{R}$ that is strictly increasing and
$\alpha(0)=0$.

###  II-B Occupancy Grid Map

Occupancy grid maps (OGMs) are fundamental tools in mobile robotics for
mapping and navigation. They represent the environment as a tessellated,
multi-dimensional grid, with each cell quantifying whether a part of the
environment is occupied or not [18]. This setup, given the state of the
system, incorporates a probabilistic approach to account for uncertainties in
measurements and the robot’s pose estimation [19], [20]. OGMs are typically
updated using sensor data from LiDAR, sonar, or vision systems and the robot’s
movements over time. This process can be employed using various frameworks,
including probabilistic [18], Bayesian [21], neural network-based [22], and
forward modeling approaches [23].

The process of dynamically generating an occupancy grid map based on the
current robot state and sensor inputs can be succinctly described as follows:

$m^{\prime}=M(\boldsymbol{x},z,m),$ (4)

where $\boldsymbol{x}$ represents the robot’s current state, $z$ denotes the
measurements from onboard perception sensors, and $m$ is the existing map. The
function
$M:\mathbb{R}^{n}\times\mathbb{R}^{k}\times\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$
updates the current map belief $m$ to $m^{\prime}$, by integrating new sensor
data at the latest state of the robot into the map.

This formulation reflects that the occupancy grid map continuously evolves,
capturing real-time changes in the environment as detected by the sensors and
the movement of the robot.

In this paper, the occupancy grid map $m$ is converted into a binary
representation. In the binary OGM, each occupied cell signifies an obstacle in
the environment. Consequently, the overall unsafe set, denoted by
$\mathbb{X}_{u}$, is defined as the union of all states corresponding to the
occupied cells of the map. Trivially, the final safe set $\mathbb{X}_{s}$ is
given by $\mathbb{X}\setminus\mathbb{X}_{u}$, which includes all states where
the robot is located in non-occupied cells of the map. The process of deriving
the safe set $\mathbb{X}_{s}$ from map $m$ can be described briefly as
follows:

$\mathbb{X}_{s}=M_{b}(m)$ (5)

###  II-C Signed Distance Function

A signed distance function (SDF) is a function that, for a given point $p$,
outputs the signed distance from $p$ to the closest boundary of a set $\Omega$
in a metric space. The sign of the distance indicates whether $p$ is inside or
outside $\Omega$ [24]. In this paper, this function is defined to take
positive values for points inside $\Omega$, decrease in value as $p$
approaches the boundary where the distance is zero, and take negative values
outside $\Omega$. If $\Omega$ is a subset of a metric space with metric $d$,
the described signed distance function is defined as follows:

$\displaystyle\phi(p,\Omega)\equiv$ $\displaystyle d(p,\partial\Omega)$ if
$p\in\Omega$ (6a) $\displaystyle\phi(p,\Omega)\equiv$
$\displaystyle-d(p,\partial\Omega)$ if $p\notin\Omega$ (6b)

where $\partial\Omega$ denotes the boundary of the set $\Omega$ and
$d(p,\partial\Omega)$ represents the minimum distance from $p$ to points on
the boundary. Additionally, $d(.,\partial\Omega)$ is the desired distance
metric to the set $\partial\Omega$.

In the case where $\Omega$ is a subset of Euclidean space $\mathbb{R}^{n}$
with a piecewise smooth boundary, and $d(.,\partial\Omega)$ is the Euclidean
distance metric, the signed distance function is differentiable almost
everywhere. Moreover, the gradient of the function satisfies the eikonal
equation $|\nabla\phi|=1$ at every point where it is differentiable [25].

##  III CONTROL METHODOLOGY

In this section, we introduce a novel method for defining CBF using occupancy
grid maps. Like traditional CBFs, the designed CBF can be applied using
various formulations, such as a safety filter or a CBF-CLF QP. For our
experimental evaluations, we have chosen the CBF-CLF QP formulation. Essential
details of this formulation are succinctly provided in this section.

###  III-A Occupancy Grid Map-based Control Barrier Function (OGM-CBF)

Assuming the online occupancy grid map $m$ updating as defined in (4), we will
show that the safety is achieved using a control barrier function defined as
follows:

$\displaystyle h(\boldsymbol{x})$
$\displaystyle=\Phi_{s}(\boldsymbol{x},m)+l_{s}+l_{a}cos(\eta)$ (7)

where:

  * •

$\Phi_{s}(\boldsymbol{x},m)=T(\phi(\boldsymbol{x},M_{b}(m)))$, and
$T:\mathbb{R}\rightarrow\mathbb{R}$ is a scalar odd function for shaping and
smoothing. In other words, $\Phi_{s}$ is a smoothened signed distance function
of the set $\mathbb{X}_{s}$, evaluated at the robot’s body frame position $p$
in the map, noting that $p$ is part of the state vector $\boldsymbol{x}$.

  * •

$\eta$ is the angle between vectors $\hat{\boldsymbol{x}}$, and
$\nabla_{p}\Phi_{s}(\boldsymbol{x},m)$, where $\hat{\boldsymbol{x}}$ is the
unit vector of the first coordinate of the body frame,
$\nabla_{p}\Phi_{s}(\boldsymbol{x},m)$ is gradient of
$\Phi_{s}(\boldsymbol{x},m)$.

  * •

$l_{s}$, and $l_{a}$ are scalar values such that $0<l_{a}\leq-l_{s}$.

Note that in above equations, all the quantities such as body frame and map
frame, and accordingly $\Phi_{s}$ are defined in the same coordinate frame,
for instance, the robot ego-centric frame.

Figure 1: Top-view snapshot of the environment in CARLA. The second row shows
the constructed occupancy grid map, and the last row shows
$\Phi_{s}(\boldsymbol{x},m)$ and its gradient field in the map frame for all
the states in the map.

Next, we provide some intuition of the proposed CBF in (7). The proposed CBF
function comprises two main components: i) the effect of the robot’s distance
to obstacles, and ii) the angle between the robot’s heading and the gradient
of the SDF function $\Phi_{s}$.

In the first component, $\Phi_{s}(\boldsymbol{x},m)$, map $m$ captures the
latest location, shape, and size of the obstacles within the map at each time
step, as outlined in (4). In the map, each occupied cell indicates an
obstacle. Although effective, this map lacks smooth transitional values. In
Fig. 1, the top row displays a top view of the scenario and obstacle status,
while the second row depicts the resulting binary occupancy grid map. Function
$\phi$ is smooth almost everywhere except at the curves where the closest
obstacle region switches. Applying $T(.)$ to $\phi$ ensures smooth $\Phi_{s}$
everywhere. The last row in Fig. 1 illustrates the smooth $\Phi_{s}$ in the
map coordinate.

The second component, $l_{s}+l_{a}cos(\eta)$, involves calculating the cosine
similarity between two vectors: $\hat{x}$, representing the robot heading, and
$\nabla_{p}\Phi_{s}(\boldsymbol{x},m)$. Notice that $\nabla_{p}\Phi_{s}$ is
the direction that $\Phi_{s}$ increases the most, thus moving in that
direction will move the robot further away from obstacles. The gradient field
depicted in the last row of Fig. 1 visualizes
$\nabla_{p}\Phi_{s}(\boldsymbol{x},m)$ within the map frame. Notice that
$l_{s}$ is a negative number and the value $-l_{s}-l_{a}\geq 0$ is the safety
margin, a design choice.

Now, we can establish how the constructed function has the necessary
properties to be a control barrier function.

###### Theorem 1

Set $S=\left\\{\boldsymbol{x}\in\mathbb{X}\;|\;h(\boldsymbol{x})\geq
0\right\\}$ with $h(\boldsymbol{x})$ defined in (7) is a forward-invariant set
and $h(\boldsymbol{x})$ is a control barrier function for system (1), with
$\boldsymbol{u}=[\mathrm{v},\boldsymbol{\omega}]^{T}$ as the control input.

###### Proof:

We will first show that $S\subseteq\mathbb{X}_{s}$. The condition
$0<l_{a}\leq-l_{s}$ ensures that $\forall\boldsymbol{x}$ we have
$l_{s}+l_{a}cos(\eta)\leq 0$. For any $\boldsymbol{x}\in
S=\\{\boldsymbol{x}:h(\boldsymbol{x})\geq 0\\}$, we have
$h(\boldsymbol{x})\geq 0$, therefore $\Phi_{s}(\boldsymbol{x})\geq 0$, and
consequently $\phi(\boldsymbol{x})\geq 0$, since $T(.)$ is an odd function.
Thus
$\boldsymbol{x}\in\mathbb{X}_{s}=\\{\boldsymbol{x}:\phi(\boldsymbol{x})\geq
0\\}$, and we have shown the first part.

We will now show that there are control signals $\boldsymbol{u}$ such that (3)
is satisfied. It is easy to see that $\dot{h}=\frac{\partial\Phi_{s}}{\partial
x}\dot{p}+l_{a}\sin(\eta)\dot{\eta}$, which can be written as
$\dot{h}=a(x)v+b(x)\boldsymbol{\omega}$, noting that $\dot{p}$ and
$\dot{\eta}$ provide the $v$ term, while $\dot{\eta}$ provides the
$\boldsymbol{\omega}$ term, assuming a kinematic model. Since
$(v,\boldsymbol{\omega})$ can be made zero, condition (3) can be satisfied and
thus $h(\boldsymbol{x})$ is a control barrier function for (1), and $S$ is an
invariant set. ∎

Considering the control system in (1), with control inputs
$\boldsymbol{u}=[\mathrm{v},\boldsymbol{\omega}]^{T}$, the control policy
outlined in (8) ensures safety across all trajectories originating within $S$.

$\pi_{ogm-
cbf}=\\{\boldsymbol{u}\in\mathbb{U}\;|\;\dot{h}(\boldsymbol{x},\boldsymbol{u})\geq-\alpha(h(\boldsymbol{x}))\\}$
(8)

It is important to recognize that the efficacy of the safety controller is
influenced by the configuration of the function $T(.)$, $\alpha(.)$, and the
values of $l_{s}$, and $l_{a}$. which can be adjusted according to the model
of the system, the dimensions of the map, and the environment.

Note that since the union of all occupied cells is considered in the function
$h(\boldsymbol{x})$, and the distance to the closest occupied cell is being
calculated in this function, the OGM-CBF automatically will find the arbitrary
safe set $S$, that considers all the detected obstacles. This means, in
contrast to other methods, one CBF construction is enough at each timestep to
guarantee staying in the safe set $S$, which makes this method computationally
efficient.

Figure 2: Top-view snapshot of simulation in the CARLA environment. The blue
car is a self-driving vehicle using OGM-CBF to navigate safely. The orange
dashed line is the path of the autonomous car. The second row shows the upper
half of the ego-centric local occupancy grid map denoted as $m$, and the last
line illustrates $\Phi_{s}(\boldsymbol{x},m)$ and its gradient field in the
map frame for all the states in the map.

###  III-B Including Robot shape in OGM-CBF

OGM-CBF is capable of considering the robot’s size and shape using established
inflation methods on occupancy grid maps. For example, the robot can be
modeled as a union of multiple spheres, similar to that described in [26]. In
addition, the occupied cells in the map can be inflated by the radius of the
spheres. Then the designed CBF can be analyzed at the centers of these
spheres.

One can also choose number of points at the boundary of the robot to be
avoided, with one constraint per point in OGM-CBF. Having enough points, the
OGM-CBF ensures that the entire boundary of the robot remains within the safe
set $S$.

###  III-C Control Lyapunov Function (CLF)

In addition to ensuring safety, it is imperative to establish a control
objective. Consider a candidate Control Lyapunov Function (CLF) for this
purpose, denoted as $V:\mathbb{X}\to\mathbb{R}$, characterized by continuous
differentiability and positive definiteness.

###### Definition 5

Consider a function $V:\mathbb{X}\to\mathbb{R}$ to be positive-definite over
the domain $\mathbb{X}$ provided that: ($i$) $V(0)=0$, and ($ii$) for any
$\boldsymbol{x}\in\mathbb{X}$ with $\boldsymbol{x}\neq 0$, it holds that
$V(\boldsymbol{x})>0$.

As mentioned in [2], the stability of (1) in equilibrium point can be achieved
with a necessary and sufficient condition as follows:

$\inf_{\boldsymbol{u}\in\mathbb{U}}[\dot{V}(\boldsymbol{x},\boldsymbol{u})=\frac{\partial
V}{\partial\boldsymbol{x}}(f(x)+g(x)\boldsymbol{u})]\leq-\gamma(V(\boldsymbol{x}))$
(9)

where $\gamma(.)$ is defined as an extended class $\kappa$ function.
Consequently, the subsequent control strategy will ensure the stability
requirements:

$\pi_{clf}=\\{\boldsymbol{u}\in\mathbb{U}_{s}\;|\;\dot{V}(\boldsymbol{x},\boldsymbol{u})\leq-\gamma(V(\boldsymbol{x}))\\}$
(10)

###  III-D Integrating CLF and OGM-CBF

In order to fulfill the criteria for both safety and stability, the proposed
control policy should represent the intersection of two preceding policies,
denoted as $\pi^{*}=\pi_{ogm-cbf}\cap\pi_{clf}$. In instances where this
intersection is empty, prioritization of safety becomes imperative.
Consequently, the following constrained optimization problem can be defined to
find an optimal control input that simultaneously satisfies the constraints
related to safety and stability, respectively as hard and soft constraints.

$\displaystyle\boldsymbol{u}^{*}$
$\displaystyle=\underset{\boldsymbol{u}}{\mathrm{argmin}}\;J(\boldsymbol{u})$
(11a) $\displaystyle\text{s.t.:}\quad\dot{V}(\boldsymbol{x},\boldsymbol{u})$
$\displaystyle\leq-\gamma(V(\boldsymbol{x}))+\delta$ (11b)
$\displaystyle\dot{h}(\boldsymbol{x},\boldsymbol{u})$
$\displaystyle\geq-\alpha(h(\boldsymbol{x}))$ (11c)
$\displaystyle\boldsymbol{u}_{lb}$
$\displaystyle\leq\boldsymbol{u}\leq\boldsymbol{u}_{ub}$ (11d)

where $J$ represents the cost function, while $\boldsymbol{u}_{lb}$ and
$\boldsymbol{u}_{ub}$ are the lower and upper bounds of the control input,
respectively. Additionally, $\delta$ is identified as the relaxation
parameter, which is engaged when the constraints are in conflict, leading to
an absence of elements in the intersection $\pi_{ogm-cbf}\cap\pi_{clf}$. Under
such circumstances, the hard constraint (related to safety) is maintained,
whereas the soft constraint (related to stability) is moderated. It is
customary to set a quadratic cost $J(\boldsymbol{u})$, which renders the
optimization a Quadratic Programs (QP), since the constraints are linear in
control $u$.

In the Algorithm (1) steps of deployment of the proposed map-based safe
control is expressed.

Input: sensor data: $z$, state of the system: $\boldsymbol{x}$

Output: Control Input: $\boldsymbol{u}$

1 while  _Robot is running_ do

2

3 receive sensor data $z$ and state of the system $\boldsymbol{x}$

4 generate occupancy grid map $m$

5

6 generate $\Phi_{s}(\boldsymbol{x},m)$ and its gradient

7 calculate $h(\boldsymbol{x})$, and $\dot{h}(\boldsymbol{x},\boldsymbol{u})$
for the current state of the robot

8 run optimization (11) to obtain control inputs: $\boldsymbol{u}^{*}$

Algorithm 1 Map-based Safety Control Algorithm

##  IV function shaping and smoothing

In this section, we aim to explore the properties of the function $T(.)$ used
in generating $h(\boldsymbol{x})$ as a control barrier function. Given that
the signed distance function is continuously differentiable almost everywhere,
as discussed in section II-C, the design of $T(\phi(\boldsymbol{x},M_{b}(m)))$
ensures that $h(\boldsymbol{x})$ maintains continuous differentiability across
the entire state space. Notably, $\phi(.)$ may lack differentiability at non-
smooth points or centers of the set $\Omega$. To address this,
$T(\phi(\boldsymbol{x},M_{b}(m)))$ translates the grid map values of
$\phi(\boldsymbol{x},M_{b}(m))$ to the state space and employs polynomial
interpolation $p(.)$, thus guaranteeing continuous differentiability at all
points in the state space. This approach ensures that $h(\boldsymbol{x})$ is a
continuously differentiable function, which is required as specified in (2).
As stated before, care needs to be taken as $T(.)$ needs to be an odd
function.

Moreover, $T(.)$ enhances the control over the shape of $h(\boldsymbol{x})$ by
adjusting the gradient behaviors and scaling the values of
$\phi(\boldsymbol{x},M_{b}(m))$ through non-linear modifications prior to
polynomial interpolation. This adjustment is tailored to fit the dynamics of
the system described in (1).

Figure 3: $h(\boldsymbol{x})$ (first),
$\dot{h}(\boldsymbol{x},\boldsymbol{u})+\alpha(h(\boldsymbol{x}))$ (second),
Linear (third) and angular velocity (forth) of the autonomous vehicle vs. time
in the CARLA simulator. The devised commands by the controller are visualized
with a solid line, while real measurements are shown using a dashed line.

Figure 4: 2D snapshot of the robot during the experiment (left), the
constructed occupancy grid map using onboard sensors (middle), and
$\Phi_{s}(\boldsymbol{x},m)$ and its gradient field in the map frame for all
the states in the map (right). The orange line shows the path of the robot and
the red gradients are the taken gradient values to the algorithm during the
experiment.

##  V Experimental Evaluation

In this section, we will demonstrate the effectiveness of the introduced
control strategy on a self-driving car in the CARLA autonomous driving
simulator[27] and on an industrial mobile robot.

###  V-A CARLA

The experiment aims to validate the safe navigation of an autonomous car using
OGM-CBF in a scenario including diverse obstacles. The control inputs are
$\boldsymbol{u}=[\mathrm{v},\boldsymbol{\omega}]^{T}$, and steering and
throttle commands are calculated using the standard kinematic bicycle model in
a 2D workspace [28]. A LiDAR sensor is mounted on the ego vehicle to generate
the occupancy grid map $m$ online, where each grid cell represents a $0.1$ by
$0.1$ meter area of the environment. The simulation involves navigating an
unknown environment cluttered with various stationary obstacles, such as cars
and pedestrians.

In this experiment, safe navigation is deployed by a controller based on
Algorithm (1), where a quadratic Lyapunov function is used to align the
vehicle’s heading with a target heading, while the optimization’s cost
function maintains a desired speed $\mathrm{v}_{d}=2m/s$. The controller’s
effectiveness can be inspected from Fig. 2. In this Figure, in the first row,
the orange line shows the path of the ego vehicle. In addition to cars,
smaller obstacles such as pedestrians, bicyclists, and motorcyclists are
spawned in the environment (marked with red circles). The second row
illustrates the upper half of the ego-centric local occupancy grid map $m$
during the experiment, and the last row displays $\Phi_{s}(\boldsymbol{x},m)$
and the gradient field of it in the map frame for all the states in the map.
It should be noted that in Fig. 2, the ego vehicle initiates its journey from
the right side of the page and proceeds towards the left. Additionally, as the
robot has a positive lower bound velocity, $T(\phi(\boldsymbol{x},M_{b}(m)))$
is only constructed for the upper half of the ego-centric map as the function
$p(a\tanh(b\phi(\boldsymbol{x},M_{b}(m))))$, where $a$ and $b$ are tuned based
on the desired behaviour, and $p(.)$ is a polynomial interpolation function as
described in section IV.

Fig. 3 quantitatively illustrates that the designed controller in fact
navigates safely during the experiment. The first row displays the value of
$h(\boldsymbol{x})$, and the second row shows
$\dot{h}(\boldsymbol{x},\boldsymbol{u})+\alpha(h(\boldsymbol{x}))$, which
remains non-negative, indicating constant satisfaction of the safety
constraint. The third and fourth rows depict the control inputs
$\boldsymbol{u}=[\mathrm{v},\boldsymbol{\omega}]^{T}$ generated by Algorithm
(1) (solid lines), alongside the actual linear and angular velocities of the
vehicle (dashed lines) with respect to time. Note the value of
$h(\boldsymbol{x})$ can become negative near obstacles due to the
approximation of the system (1) as a kinematic model.

###  V-B Real-world

Besides simulated validation, additional experiments were conducted on an
industrial mobile robot (MiR100), which is a differential wheeled robot
equipped with an Nvidia Jetson AGX Orin and two laser scanners. Mapping and
safe control are deployed using the ROS2 framework on this platform, and
localization is done internally using onboard sensors.

In this experiment, the robot utilizes a pre-established global map of the
environment, which is dynamically updated with an average frequency of
$16.372$ Hz using real-time sensor measurements collected from laser scanners.
Each grid cell in the map covers a $0.05$ by $0.05$ meter section of the
environment. Safe navigation is managed by a controller based on Algorithm (1)
using OGM-CBF, with a quadratic Lyapunov function for maintaining the desired
heading. Moreover, the optimization’s cost function is designed to keep the
speed of the robot at $\mathrm{v}_{d}=0.25$ m/s. The employed controller runs
on the platform with an average frequency of $12.498$ Hz, and the applied
function $T(\phi(\boldsymbol{x},M_{b}(m)))$ is the same as in the CARLA
experiments.

Figure 5: $h(\boldsymbol{x})$ (first),
$\dot{h}(\boldsymbol{x},\boldsymbol{u})+\alpha(h(\boldsymbol{x}))$ (second),
Linear (third) and angular velocity (forth) of the mobile robot vs. time. The
controller’s commands are displayed using a solid line, whereas the real
measurements are indicated with a dashed line.

Fig. 4 shows the robot’s path (orange lines) during the experiment. The left
image demonstrates OGM-CBF’s capabilities to avoid obstacles with arbitrary
shapes (the blue obstacle and a human). The middle image shows the global
occupancy grid map updated in real-time by laser scanners, and the robot’s
path in the map frame. The right image presents $\Phi_{s}(\boldsymbol{x},m)$
and its gradient field in the map frame for all the states in the map (blue
arrows). Moreover, the gradients of $\Phi_{s}(\boldsymbol{x},m)$ along the
robot’s path are marked by red arrows, indicating the direction of avoidance
near obstacles utilized in the controller. Note that the robot progresses from
the top to the bottom of the images throughout the experiment.

Fig. 5 quantitatively validates the robot’s safe navigation using OGM-CBF
during the experiment. The first row displays the value of
$h(\boldsymbol{x})$, the second row presents
$\dot{h}(\boldsymbol{x},\boldsymbol{u})+\alpha(h(\boldsymbol{x}))$, and the
final two rows illustrate the commanded (solid blue line) and actual (dashed
red line) linear and angular velocities of the robot with respect to time.

As values in Fig. 5 illustrate, the robot is able to satisfy the safety
constraint of OGM-CBF and navigate safely in the environment. Note that as the
system (1) is considered kinematically, $h(\boldsymbol{x})$ goes to negative
values in proximity to obstacles. However, the controller increases the value
of $h(\boldsymbol{x})$ successfully by generating proper control commands.

##  VI CONCLUSION

This paper introduces a novel method, the Occupancy Grid Map-based Control
Barrier Function (OGM-CBF), which enhances safe navigation in robotics across
both known and unknown terrains. By integrating the Control Barrier Function
(CBF) with occupancy grid maps, our approach effectively accommodates
arbitrary obstacle shapes and configurations. The designed CBF is sensor-
agnostic and can be implemented using a range of hybrid sensor modalities.
Moreover, the OGM-CBF is formulated as a single constraint in a Quadratic
Programming (QP) formulation, taking into account all previously seen
obstacles. Additionally, the constructed CBF considers the robot’s size while
applying established inflation methods to the occupancy grid map.

We validated the effectiveness of OGM-CBF through experiments in the CARLA
simulator and practical implementations on an industrial mobile robot,
confirming its capability for reliable autonomous navigation. Future work will
demonstrate the integration of different sensor modalities where safe/unsafe
segmentation is integrated to one OGM, and incorporate the effect of the
actuator dynamics into the CBF design.

## References

  * [1] G.-Z. Yang, J. Bellingham, P. E. Dupont, P. Fischer, L. Floridi, R. Full, N. Jacobstein, V. Kumar, M. McNutt, R. Merrifield, et al., “The grand challenges of science robotics,” Science robotics, vol. 3, no. 14, p. eaar7650, 2018.
  * [2] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada, “Control barrier functions: Theory and applications,” in 2019 18th European control conference (ECC), pp. 3420–3431, IEEE, 2019.
  * [3] Q. Nguyen, A. Hereid, J. W. Grizzle, A. D. Ames, and K. Sreenath, “3d dynamic walking on stepping stones with control barrier functions,” in 2016 IEEE 55th Conference on Decision and Control (CDC), pp. 827–834, IEEE, 2016.
  * [4] A. D. Ames, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs with application to adaptive cruise control,” in 53rd IEEE Conference on Decision and Control, pp. 6271–6278, IEEE, 2014.
  * [5] L. Wang, A. D. Ames, and M. Egerstedt, “Safety barrier certificates for collisions-free multirobot systems,” IEEE Transactions on Robotics, vol. 33, no. 3, pp. 661–674, 2017.
  * [6] Y. Chen, A. Singletary, and A. D. Ames, “Guaranteed obstacle avoidance for multi-robot operations with limited actuation: A control barrier function approach,” IEEE Control Systems Letters, vol. 5, no. 1, pp. 127–132, 2020.
  * [7] Z. Qin, K. Zhang, Y. Chen, J. Chen, and C. Fan, “Learning safe multi-agent control with decentralized neural barrier certificates,” arXiv preprint arXiv:2101.05436, 2021.
  * [8] A. Singletary, A. Swann, Y. Chen, and A. D. Ames, “Onboard safety guarantees for racing drones: High-speed geofencing with control barrier functions,” IEEE Robotics and Automation Letters, vol. 7, no. 2, pp. 2897–2904, 2022.
  * [9] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs for safety critical systems,” IEEE Transactions on Automatic Control, vol. 62, no. 8, pp. 3861–3876, 2016.
  * [10] C. Dawson, B. Lowenkamp, D. Goff, and C. Fan, “Learning safe, generalizable perception-based hybrid control with certificates,” IEEE Robotics and Automation Letters, vol. 7, no. 2, pp. 1904–1911, 2022.
  * [11] W. Xiao, T.-H. Wang, R. Hasani, M. Chahine, A. Amini, X. Li, and D. Rus, “Barriernet: Differentiable control barrier functions for learning of safe robot control,” IEEE Transactions on Robotics, 2023.
  * [12] W. Xiao, T.-H. Wang, M. Chahine, A. Amini, R. Hasani, and D. Rus, “Differentiable control barrier functions for vision-based end-to-end autonomous driving,” arXiv preprint arXiv:2203.02401, 2022.
  * [13] K. Long, C. Qian, J. Cortés, and N. Atanasov, “Learning barrier functions with memory for robust safe navigation,” IEEE Robotics and Automation Letters, vol. 6, no. 3, pp. 4931–4938, 2021.
  * [14] M. Srinivasan, A. Dabholkar, S. Coogan, and P. A. Vela, “Synthesis of control barrier functions using a supervised machine learning approach,” in 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 7139–7145, IEEE, 2020.
  * [15] H. Abdi, G. Raja, and R. Ghabcheloo, “Safe control using vision-based control barrier function (v-cbf),” in 2023 IEEE International Conference on Robotics and Automation (ICRA), pp. 782–788, IEEE, 2023.
  * [16] J. Liu, M. Li, J. W. Grizzle, and J.-K. Huang, “Clf-cbf constraints for real-time avoidance of multiple obstacles in bipedal locomotion and navigation,” in 2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 10497–10504, IEEE, 2023.
  * [17] F. Blanchini, “Set invariance in control,” Automatica, vol. 35, no. 11, pp. 1747–1767, 1999.
  * [18] H. Moravec and A. Elfes, “High resolution maps from wide angle sonar,” in Proceedings. 1985 IEEE international conference on robotics and automation, vol. 2, pp. 116–121, IEEE, 1985.
  * [19] G. N. DeSouza and A. C. Kak, “Vision for mobile robot navigation: A survey,” IEEE transactions on pattern analysis and machine intelligence, vol. 24, no. 2, pp. 237–267, 2002.
  * [20] T. Collins and J. Collins, “Occupancy grid mapping: An empirical evaluation,” in 2007 mediterranean conference on control & automation, pp. 1–6, IEEE, 2007.
  * [21] L. Matthies and A. Elfes, “Integration of sonar and stereo range data using a grid-based representation,” in Proceedings. 1988 IEEE International Conference on Robotics and Automation, pp. 727–733, IEEE, 1988.
  * [22] S. B. Thrun, “Exploration and model building in mobile robot domains,” in IEEE international conference on neural networks, pp. 175–180, IEEE, 1993.
  * [23] S. Thrun, “Learning occupancy grid maps with forward sensor models,” Autonomous robots, vol. 15, pp. 111–127, 2003.
  * [24] T. Chan and W. Zhu, “Level set based shape prior segmentation,” in 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), vol. 2, pp. 1164–1170, IEEE, 2005.
  * [25] C. Dapogny and P. Frey, “Computation of the signed distance function to a discrete contour on adapted triangulation,” Calcolo, vol. 49, pp. 193–219, 2012.
  * [26] A. Muller, “Energy optimal control of serial manipulators avoiding collisions,” in Proceedings of the IEEE International Conference on Mechatronics, 2004. ICM’04., pp. 299–304, IEEE, 2004.
  * [27] A. Dosovitskiy, G. Ros, F. Codevilla, A. Lopez, and V. Koltun, “Carla: An open urban driving simulator,” in Conference on robot learning, pp. 1–16, PMLR, 2017.
  * [28] R. Rajamani, Vehicle dynamics and control.  Springer Science & Business Media, 2011.