# MegLoc: A Robust and Accurate Visual Localization Pipeline

Shuxue Peng111Equally contributed, order determined by a team lottery. ,
Zihang He111Equally contributed, order determined by a team lottery. , Haotian
Zhang111Equally contributed, order determined by a team lottery. , Ran
Yan111Equally contributed, order determined by a team lottery. , Chuting
Wang111Equally contributed, order determined by a team lottery. , Qingtian
Zhu111Equally contributed, order determined by a team lottery. , Xiao Liu
Megvii 3D

###### Abstract

In this paper, we present a visual localization pipeline, namely MegLoc, for
robust and accurate 6-DoF pose estimation under varying scenarios, including
indoor and outdoor scenes, different time across a day, different seasons
across a year, and even across years. MegLoc achieves state-of-the-art results
on a range of challenging datasets, including winning the Outdoor and Indoor
Visual Localization Challenge of ICCV 2021 Workshop on Long-term Visual
Localization under Changing Conditions, as well as the Re-localization
Challenge for Autonomous Driving of ICCV 2021 Workshop on Map-based
Localization for Autonomous Driving.

##  1 Introduction

Visual localization is a key technology for applications such as Augmented,
Mixed, and Virtual Reality, as well as for robotics. It aims to estimate 6-DoF
camera poses for query images by aid of a series of mapping images with given
ground truth poses. Mapping images and query images suffer large appearance
variations caused by seasonal and illuminational changes.

MegLoc follows a conventional two-stage localization fashion, namely mapping
and localization. The two stages will get respectively elaborated in the
following sections. MegLoc can function as a complete visual localization
pipeline, while also follows a modular design pattern, where each module can
be replaced individually. For example, it is possible to use MegLoc’s
localization module together with another mapping method. The overall pipeline
of MegLoc is illustrated in Fig. 1.

Figure 1: Illustration of the whole pipeline of MegLoc. Specifically, the
pipeline is divided into two stages, namely, mapping and localization. For
mapping stage, we first do reference matching, then triangulate the matches
using COLMAP [20], and finally we refine the 3D map by depth verification and
uncertainty rejection. For localization stage, we retrieve those reference
images similar to the current query image by global descriptors fusion. We
perform local feature matching between reference images and the query image.
After that we leverage PnP algorithm to compute initial poses and rerank the
poses. A final pose optimization step is also performed to obtain the final
poses.

##  2 Features

###  2.1 Image Pre-processing

##### Resizing

Empirically, each image is resized to make the larger side 1600 pixels while
remaining its original aspect ratio. Some pixels of the bottom-right border
might get cropped, to ensure each dimension to be a multiple of 8.

##### Masking

Frequently dynamic objects introduce unreliable and unrepeatable keypoints to
our solution. To alleviate this problem, we perform semantic segmentation on
images and mask out dynamic classes. The selected model is a DeepLab-V3+
network [4] pretrained on CityScapes dataset [5]. The definition of dynamic
classes is provided by CityScapes officially. For cases that the bottom part
of the camera view is always occluded, _e.g_. by the car shell, local features
within this area impede subsequent feature matching and triangulation. We
remove local features in this area.

##### Undistortion

Due to positional encoding, SuperGlue [19] is known to be sensitive to lens
distortion. Undistortion of the images before matching is proven by
experiments to generate more accurate matches.

###  2.2 Local Feature

SuperPoint [7] and ASLFeat [14] are selected as local feature extractors to
fully utilize their complimentary characteristics. Specifically, SuperPoint is
trained on COCO dataset [12], and it generally works fairly well on different
datasets. In contrast, ASLFeat has a more powerful backbone which involves DCN
[6, 24] for shape-awareness local feature extraction. Also, ASLFeat performs
multi-level feature fusion.

The Non-Maximum Suppression (NMS) radius and keypoints threshold are tuned
empirically based on the experiments conducted in the CVPR 2021 Image Matching
Challenge [2]. As a result, approximately 1500 keypoints are extracted for
each image on average.

###  2.3 Global Feature

##### NetVLAD

Image retrieval has long been a fundamental task in computer vision. Given a
query image, a robust image retrieval algorithm is supposed to retrieve the
most similar images in the image database, which are called reference images
or candidates. This task is challenging because the retrieved results are
sensitive to the change of illuminations, shape or even the surroundings under
certain circumstances. Efforts have been made to improve the performance.

Traditional algorithms are basically resorting to the local features such as
SIFT [13] or ORB [18], by clustering algorithms such as K-Means. And then,
they embed the current image by Bag of Words or TF-IDF using the clustering
result of local features to a specified length of vector.

VLAD is also an excellent traditional image retrieval algorithm. Given an
image and its local features of the shape $N\times D$, where $N$ and $D$
respectively denote the number of local features and the length of each local
feature vector, VLAD embeds the local feature matrix to a certain
representation matrix of shape $K\times D$. $K$ is the predefined number of
clustering centroids. Formally, the representation matrix is computed as
follows:

$V(j,k)=\sum_{i=1}^{N}a_{k}(x_{i})(x_{i}(j)-c_{k}(j)),$

where $x_{i}(j)$ and $c_{k}(j)$ are the $j$-th dimension of $i$-th descriptor
and $k$-th centroid respectively. $a_{k}(x_{i})$ is a binary signal function
whose value is 1 if the current descriptor belongs to the centroid or cluster,
or 0 otherwise.

NetVLAD [1], literally, leverages CNNs to obtain a global feature of input
image in an end-to-end fashion. The original VLAD breaks the continuities of
feature extraction because of the hard assignment of $a_{k}(x_{i})$. To ensure
the whole process of feature extraction differentiable, NetVLAD reformulates
the extracting equation in original VLAD as follows:

$\overline{a_{k}}(x_{i})=\frac{e^{-\alpha\|x_{i}-c_{k}\|^{2}}}{\sum_{k^{\prime}}{{e^{-\alpha\|{x_{i}-c_{k{{}^{\prime}}}}\|^{2}}}}}.$

The equation above assigns the weight of descriptor $x_{i}$ to cluster $c_{k}$
according to the proximity. The more the descriptor approximates the
centroids, the larger the value $\overline{a_{k}}(x_{i})$ will be. Apparently,
$\overline{a_{k}}(x_{i})$ ranges between $0$ and $1$, so the function provides
a soft way to obtain the global information and thus can be trained end-to-
end. Furthermore, one can expand the square in the function, and we can get
the final mathematical definition of NetVLAD:

$V(j,k)=\sum_{i=1}^{N}\frac{e^{w^{T}_{k}+b_{k}}}{\sum_{k^{\prime}}e^{w^{T}_{k^{\prime}}x_{i}+b_{k^{\prime}}}}(x_{i}(j)-c_{k}(j)),$

where $w_{k}$, $b_{k}$ and $c_{k}$ are learnable parameters.

##### Fusion of Multiple Global Features

According to our experiments, we find that the retrieval results are only
partial because of the nature of NetVLAD. For example, in the dataset of
RobotCar, NetVLAD tends to retrieve the images that are under the similar
illumination condition as the query image, which hinders the final
performance. To alleviate this problem, we resort to a fusion result of global
features to make global description more robust. We fuse NetVLAD [1], DELG
[3], APGeM [17, 21] and OpenIBL [8] features as the final global
representation.

##### Reranking of Global Features

Inspired by DELG [3] and Patch-NetVLAD [9]. We use local feature matching to
rerank the retrieval results obtained. To be more specific, given a query
image, we first retrieve $M$ images by global descriptors. Then, we leverage
local features of these images and SuperGlue [19] to conduct local feature
matching. We sort the retrieved images by the number of valid correspondence
between the retrieved images and the query image. Finally, we select the top
$N$ retrieval images as the final retrieval result.

###  2.4 Matching

We retrain SuperGlue [19] together with its official feature extractor
SuperPoint [7] in an end-to-end manner on MegaDepth dataset [11]. More
specifically, we split the original SuperPoint [7] into two sub-networks. The
first one is frozen with the official weights to extract keypoints from
images, while the other one is fine-tuned to provide feature descriptions.
However, we find this adjustment only advances the model performance slightly,
since SuperGlue [19] can already match the given points pretty well.

####  2.4.1 Guided Pyramid Matching

For those corner cases that the number of matches found by SuperGlue [19] is
less than 100, a pyramid extraction strategy is then applied, _e.g_. multiple
scales and/or multiple orientations and we might combine the matches in
different scales (ALL) or trust the one with the most number of matches (MAX).

##  3 Mapping

This section introduces MegLoc’s mapping strategy under different
circumstances.

###  3.1 Sparse Reconstruction

For circumstances without given SfM models, we perform sparse reconstruction
on given mapping images and utilize the reconstructed SfM model for
localization. The reconstruction consists of 4 stages, namely retrieval of
image pairs, local feature matching, triangulation of 2D keypoints and final
map refinement.

####  3.1.1 Image Pairs Retrieval

In order to triangulate 2D points into 3D space, we need to acquire image
pairs with shared observations. The retrieval strategies vary according to
different scenarios.

##### By Poses

Since mapping images are provided with ground truth poses, we can assume that
spatially adjacent poses are more likely to have shared observations.

##### By Co-visibility

The most naive strategy is to retrieve image pairs by making use of image-
level co-visibility information.

##### By Global Descriptors

The retrieval of image pairs can naturally fit into a typical image retrieval
pipeline, where images are retrieved by ranking pairwise cosine similarity of
global descriptors. NetVLAD [1] is chosen as the global descriptor. We also
find that fusing multiple global features together in a similar way as [10] is
helpful.

##### By Temporal Sequence

For datasets delivered in a temporally-sequential order, we can apparently
make the assumption that temporally neighboring images tend to share
observations.

Figure 2: An example of 4Seasons dataset [22]. Left: original reconstruction
of map. Right: reconstructed map with outliers rejection by removing large
uncertainty map points.

####  3.1.2 3D Map Point Uncertainty & Map Refinement

A 3D map point $p_{w}$ in the world frame can be observed in multiple cameras
$C_{i}(i=1...N)$ from different views. The camera pose $C_{i}$ and its
intrinsic parameters can be represented as $[R_{wc_{i}},t_{wc_{i}}]$ and
$[f_{x},f_{y},c_{x},c_{y}]$, respectively. We assume that a 3D map point
$P_{c}=[x,y,z]$ with respect to the camera reference frame corresponds to a
observation in the image plane of camera $C_{i}$, and its 2D location can be
represented as $p_{uv}$.

The Jacobian of 2D observation $dp_{uv}$ to the 3D map point $dp_{w}$ is

$\displaystyle J_{i}$
$\displaystyle=\frac{dp_{uv}}{dp_{w}}=\frac{\frac{dp_{uv}}{dp_{c}}\cdot
dp_{c}}{dp_{w}}$ (1)
$\displaystyle=\begin{bmatrix}\frac{f_{x}}{z}&0&-\frac{x\cdot f_{x}}{z^{2}}\\\
0&\frac{f_{y}}{z}&-\frac{y\cdot f_{y}}{z^{2}}\end{bmatrix}\cdot R_{wc_{i}}.$

Further assuming the observation uncertainty of each pixel on the image plane
is an identity matrix, $\Sigma_{uv}=\begin{bmatrix}1&0\\\ 0&1\end{bmatrix}$,
then the information matrix of 3D point $p_{w}$ is

$\Sigma^{-1}_{w_{i}}=J^{T}_{i}\ast\Sigma^{-1}_{uv}\ast J_{i}.$ (2)

The total uncertainty is an addition of the uncertainty of all observations,
defined as

$\Sigma^{-1}_{w}=\Sigma^{-1}_{w_{1}}+\Sigma^{-1}_{w_{2}}+\ldots$ (3)

The uncertainty of three orthogonal directions can be obtained by eigen-
decomposition, and those map points with large uncertainty can be removed from
the map by setting a threshold. The value of the threshold should be
determined based on many factors, such as the scale of the reconstructed
scene, the number of mapping images, resolution of images. Fig. 2 shows the
effect of map refinement.

##  4 Localization

###  4.1 Image Pairs Retrieval

The retrieval of image pairs at localization stage is identical to Sec. 3.1.1,
but without ground truth poses.

###  4.2 Cluster-wise Localization

In this section, we illustrate our cluster-wise camera localization strategy.

####  4.2.1 Camera Clustering and Reranking

To resolve the issues caused by repetitive patterns, we purposed two ways of
database image clustering, by ground truth 3D poses or by their 3D points co-
visibility.

For clustering by pose, the database images are clustered simply based on
their spatial proximity. We only consider the positions while ignore the
orientation, since two database images with no overlapping viewing area might
still share different co-visible areas with the query image respectively. For
clustering by co-visibility, if the number of 3D points co-visible by two
database images above a certain threshold, those two images are clustered.
This strategy is illustrated in Fig. 3.

Figure 3: Candidate images clustering by pose (left) by 3D points co-
visibility (right).

Many clusters may occur and we estimate the query image pose for each, rather
than leveraging all 2D-3D matches into the calculation of PnP.

After that, all potential poses get reranked according to certain criteria,
such as number of inliers, and only those database images within the top-1
cluster are kept for final pose refinement. Here we present a specific
strategy of camera reranking used under indoor scenes.

In image retrieval, image-wise similarity is measured by global features.
These descriptors may suffer from weak texture and repeated pattern and get no
longer distinguishable especially under indoor scenes. In these circumstances,
the reliability of image-retrieval-based neighborhood searching strategy is
seriously challenged. To alleviate this problem, we first align views to be
reranked by warping them to the estimated camera pose. In this way, each
candidate cluster should have an aligned view to the query image. Then we
leverage perceptual similarity [23] which measures L2 distances between deep
feature maps and tends to measure structural similarity. The distance is
defined as

$d(x,x_{0})=\sum_{l}\frac{1}{H_{l}W_{l}}\sum_{h,w}||w_{l}\odot(\hat{y}^{l}_{hw}-\hat{y}^{l}_{0hw})||^{2}_{2},$
(4)

where $\hat{y}^{l}$ denotes the output of a given CNN encoder at layer $l$ and
$w$ is a vector for channel-wise scaling. We perform this metric between each
one of the aligned candidate views and the query image, and the cluster, whose
estimated candidate view is assigned with the lowest score, is selected as the
neighborhood for following refinement.

Figure 4: Illustration of our ICP-based and rendering-based pose refinement
strategy. (a) Rendered depths by the initial camera pose (framed with blue),
by the ICP-based refined pose (framed with red), and the given reference depth
(frameless). (b) Camera poses observed under the 3D scene. (c) Rendered RGB
images by the initial camera pose (framed with blue), by the rendering-based
refined pose (framed with red), and the given reference image (frameless).

###  4.3 Pose Refinement

Localization provides preliminary camera poses, which requires further
refinement to get more accurate pose. The strategies applied depend on
different datasets.

####  4.3.1 Iterative Closest Point (ICP)

Given dense depth maps of query image, we can refine the initial camera pose
using ICP (Interactive Closest Point). The fundamental goal of ICP is to
estimate an optimal transformation $R$ and $t$ from a point cloud to another.
By using ICP, we can optimize those preliminarily estimated camera poses
towards the ground truth geometry.

Fig. 4 (a) illustrates the camera pose refinement by ICP.

####  4.3.2 Refinement by rendering

We adopt differentiable rendering to optimize the camera pose, whose basic
principle is similar to ICP, but the optimization objective is a measurement
of photometric similarity instead. Fig. 4 (c) illustrates the camera pose
refinement by rendering.

##  5 Other Attempts

Apart from attempts made with a conventional visual localization pipeline, we
also make other attempts that take advantage of recent studies in the field of
dense reconstruction and neural rendering.

###  5.1 Multi-view Stereo (MVS)

Multi-view Stereo (MVS) aims to obtain a dense representation of an object or
a scene given a series of images. A denser representation is able to provide
more anchor points and make localization more robust.

For datasets not providing per-view depth maps, we attempt to estimate a dense
depth map for each mapping image to densify the map. Given the camera
intrinsics and extrinsics of mapping images, the problem is formulated as a
typical multi-view stereo dense reconstruction.

###  5.2 NeRF-W

Though traditional techniques such as Structure-from-Motion (SfM) have
achieved appealing success in 2D-3D mapping. The field of 3D reconstruction
has recently great progress in neural rendering. NeRF [16] models an object or
a scene implicitly by fitting a function from location and direction to color
and density, and then applies volume rendering to obtain the novel view.

So we also attempt to represent the map implicitly with the power of NeRF
[16], with which we are capable of rendering realistic synthetic images at
arbitrary views. NeRF-W [15] is proposed to synthesize novel views with
unconstrained images. It divided dynamic scene into two parts. One is “static”
and the other is “transient”. NeRF-W models the two components separately,
making it be robust to transient objects. For task of long-term visual
localization, we supposed that NeRF-W is able to handle more complicated
situations.

##  6 Conclusion

In this paper, we present MegLoc, a robust and accurate 6-DoF pose estimation
pipeline. We introduce the strategies applied by us for feature extraction,
feature matching, map refinement, image retrieval, localization and pose
refinement in detail. MegLoc has achieved state-of-the-art performance for
many visual localization benchmarks.

## References

  * [1] Relja Arandjelovic, Petr Gronat, Akihiko Torii, Tomas Pajdla, and Josef Sivic.  Netvlad: Cnn architecture for weakly supervised place recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 5297–5307, 2016.
  * [2] Xiaopeng Bi, Yu Chen, Xinyang Liu, Dehao Zhang, Ran Yan, Zheng Chai, Haotian Zhang, and Xiao Liu.  Method towards cvpr 2021 image matching challenge.  arXiv preprint arXiv:2108.04453, 2021.
  * [3] Bingyi Cao, Andre Araujo, and Jack Sim.  Unifying deep local and global features for image search.  In European Conference on Computer Vision, pages 726–743. Springer, 2020.
  * [4] Liang-Chieh Chen, Yukun Zhu, George Papandreou, Florian Schroff, and Hartwig Adam.  Encoder-decoder with atrous separable convolution for semantic image segmentation.  In Proceedings of the European conference on computer vision (ECCV), pages 801–818, 2018.
  * [5] Marius Cordts, Mohamed Omran, Sebastian Ramos, Timo Rehfeld, Markus Enzweiler, Rodrigo Benenson, Uwe Franke, Stefan Roth, and Bernt Schiele.  The cityscapes dataset for semantic urban scene understanding.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 3213–3223, 2016.
  * [6] Jifeng Dai, Haozhi Qi, Yuwen Xiong, Yi Li, Guodong Zhang, Han Hu, and Yichen Wei.  Deformable convolutional networks.  In Proceedings of the IEEE international conference on computer vision, pages 764–773, 2017.
  * [7] Daniel DeTone, Tomasz Malisiewicz, and Andrew Rabinovich.  Superpoint: Self-supervised interest point detection and description.  In Proceedings of the IEEE conference on computer vision and pattern recognition workshops, pages 224–236, 2018.
  * [8] Yixiao Ge, Haibo Wang, Feng Zhu, Rui Zhao, and Hongsheng Li.  Self-supervising fine-grained region similarities for large-scale image localization.  In European Conference on Computer Vision, pages 369–386. Springer, 2020.
  * [9] Stephen Hausler, Sourav Garg, Ming Xu, Michael Milford, and Tobias Fischer.  Patch-netvlad: Multi-scale fusion of locally-global descriptors for place recognition.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 14141–14152, 2021.
  * [10] Martin Humenberger, Yohann Cabon, Nicolas Guerin, Julien Morat, Jérôme Revaud, Philippe Rerole, Noé Pion, Cesar de Souza, Vincent Leroy, and Gabriela Csurka.  Robust image retrieval-based visual localization using kapture.  arXiv preprint arXiv:2007.13867, 2020.
  * [11] Zhengqi Li and Noah Snavely.  Megadepth: Learning single-view depth prediction from internet photos.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 2041–2050, 2018.
  * [12] Tsung-Yi Lin, Michael Maire, Serge Belongie, James Hays, Pietro Perona, Deva Ramanan, Piotr Dollár, and C Lawrence Zitnick.  Microsoft coco: Common objects in context.  In European conference on computer vision, pages 740–755. Springer, 2014.
  * [13] David G Lowe.  Distinctive image features from scale-invariant keypoints.  International journal of computer vision, 60(2):91–110, 2004.
  * [14] Zixin Luo, Lei Zhou, Xuyang Bai, Hongkai Chen, Jiahui Zhang, Yao Yao, Shiwei Li, Tian Fang, and Long Quan.  Aslfeat: Learning local features of accurate shape and localization.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 6589–6598, 2020.
  * [15] Ricardo Martin-Brualla, Noha Radwan, Mehdi SM Sajjadi, Jonathan T Barron, Alexey Dosovitskiy, and Daniel Duckworth.  Nerf in the wild: Neural radiance fields for unconstrained photo collections.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 7210–7219, 2021.
  * [16] Ben Mildenhall, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoorthi, and Ren Ng.  Nerf: Representing scenes as neural radiance fields for view synthesis.  In European conference on computer vision, pages 405–421. Springer, 2020.
  * [17] Jerome Revaud, Jon Almazán, Rafael S Rezende, and Cesar Roberto de Souza.  Learning with average precision: Training image retrieval with a listwise loss.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 5107–5116, 2019.
  * [18] Ethan Rublee, Vincent Rabaud, Kurt Konolige, and Gary Bradski.  Orb: An efficient alternative to sift or surf.  In 2011 International conference on computer vision, pages 2564–2571. IEEE, 2011.
  * [19] Paul-Edouard Sarlin, Daniel DeTone, Tomasz Malisiewicz, and Andrew Rabinovich.  Superglue: Learning feature matching with graph neural networks.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 4938–4947, 2020.
  * [20] Johannes Lutz Schönberger and Jan-Michael Frahm.  Structure-from-motion revisited.  In Conference on Computer Vision and Pattern Recognition (CVPR), 2016.
  * [21] Giorgos Tolias, Ronan Sicre, and Hervé Jégou.  Particular object retrieval with integral max-pooling of cnn activations.  In International Conference on Learning Representations (ICLR), pages 1–12, 2016.
  * [22] P. Wenzel, R. Wang, N. Yang, Q. Cheng, Q. Khan, L. von Stumberg, N. Zeller, and D. Cremers.  4Seasons: A cross-season dataset for multi-weather SLAM in autonomous driving.  In Proceedings of the German Conference on Pattern Recognition (GCPR), 2020.
  * [23] Richard Zhang, Phillip Isola, Alexei A Efros, Eli Shechtman, and Oliver Wang.  The unreasonable effectiveness of deep features as a perceptual metric.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 586–595, 2018.
  * [24] Xizhou Zhu, Han Hu, Stephen Lin, and Jifeng Dai.  Deformable convnets v2: More deformable, better results.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 9308–9316, 2019.


Moebius graphene ribbon]Quantum spin-Hall effect on the Möbius graphene ribbon

Kyriakos Flouris
Miller Mendoza Jimenez

Hans J. Herrmann

Zürich, Computational Physics for Engineering Materials, Institute
for Building Materials, Wolfgang-Pauli-Str. 27, HIT, CH-8093 Zürich

Topological phases of matter have revolutionized quantum engineering. Implementing a curved space Dirac equation solver based on the quantum Lattice Boltzmann method, we study the topological and geometrical transport properties of a Möbius graphene ribbon. In the absence of a magnetic field, we measure a quantum spin-Hall current on the graphene strip,  originating from topology and curvature, whereas a quantum Hall current is not observed. In the torus geometry a Hall current is measured. Additionally, a specific illustration of the equivalence between the Berry and Ricci curvature is presented through a travelling wave-packet around the Möbius band.

§ INTRODUCTION

Topological states of matter have been proposed theoretically and observed experimentally <cit.> triggering a surge of research in the context of quantum engineering and condensed matter. The enthusiasm is motivated by the robustness of the topological properties and their tolerance to disorder <cit.>. Especially useful in technological applications are topological currents as they are robust even in impure systems. Topological protection is a consequence of topological invariants, the most notable being the Chern number <cit.>.

A large fraction of the research is  concerned with topological insulators in 2-D <cit.>, which is gradually extending to higher dimensions <cit.> through combined systems. Graphene, Weyl semi-metals and fractional statistics particles are just a part of the materials associated to the field <cit.>. Simultaneously in condensed matter, another field of great interest  are particles constrained to move on curved surfaces.  Curvature has immense potential as a novel and pure control mechanism in quantum devices, examples of this have been proposed in graphene systems <cit.>. Graphene exhibits fascinating properties as a material but also as the platform for two-dimensional physics <cit.>. Furthermore, the honeycomb bipartite lattice of graphene has been shown to exhibit the quantum Hall and the quantum spin-Hall effects in the presence of a magnetic field <cit.>. Additionally, a graphene Möbius ribbon can exhibit topological insulator properties owning to the zig zag edge states <cit.>.

For the Möbius band, although simple to visualize, the geometry is topologically non-trivial. The Möbius strip is the simplest non-trivial fiber bundle, a bundle of the line segment over the circle with a twist. Here, we raise the question of the energy levels of a confined particle on the Möbius graphene strip and the consequence of the curvature to the transport of these particles. We simulate Dirac particles and more specifically we relate them to graphene through the low energy approximation of the band structure.  The non-relativistic quantum mechanics on a the Möbius ring have been studied  <cit.> and also non-trivial effects in the Möbius topology have  been proposed for molecular devices <cit.>. Through symmetry arguments it was proposed that the quantum Hall effect cannot exist on the Möbius but the quantum spin-Hall is possible <cit.>. This result is also confirmed in this work through a direct simulation.

A recurrent theme in recent condensed matter literature is the existence of the quantum Hall effect in the absence of a magnetic field <cit.>.  By avoiding the use of an external field, it is shown that the effect is purely a geometrical consequence.  The symmetry of the system is not, as usual, broken by a magnetic field but by an equivalent gauge field originating from the curvature effects of the manifold. In reality, these curvature effects are translated as forces due to the pseudo-potential in a strained graphene sheet.

To this end the electronic and topological properties of a Möbius graphene ribbon are explored numerically using a solver of the Dirac equation in curved space <cit.>. The method is based on the conceptual similarities between the Dirac and the Boltzmann equations and is an extension of the quantum lattice Boltzmann  method  <cit.> to curved space. Through our simulations we observe a quantum spin-Hall current in the bulk of a Möbius graphene strip.

Firstly we introduce the Dirac equation and its extension to curved space and specifically deformed graphene. In the next section, the simulations and results are presented for the Möbius and the torus geometries. The paper finishes with a summary and outlook section. An appendix can be found at the end including a more detailed description of the numerical model (Appendix <ref>) and a calculation of the Berry phase relevant to a Möbius band (Appendix <ref>).

§ THE DIRAC EQUATION IN CURVED SPACE AND GRAPHENE

The original Dirac equation <cit.> can be written in compact notation:

\begin{equation}
\label{eq:Dirac}
(i \g^{\mu}\del_{\mu}-m)\Psi=0,
\end{equation}

in natural units such that $\hbar=c=1$ for $\hbar$, Planck's constant and $c$, the speed of light, where $m$ is the particle mass, $\mu={0,1,2}$ for 2D space-time. where $\Psi = (\Psi^+, \Psi^-) = (\psi_1^+,\psi_2^-,\psi_1^-,\psi_2^+) \in \C^4$ denotes the spinor, and $\g^\mu = \g^\a e_\a^{\ \mu}$ are the generalized $\g$-matrices, where $\g^\a \in \C^{4\times 4}$ are the standard $\g$-matrices (in Dirac representation). This can be naturally extended to curved space with the addition of a covariant derivative as

\begin{equation}
\label{eq:Dirac}
(i \g^{\mu}D_{\mu}-m)\Psi=0,
\end{equation}

where $D_{\mu}$ is the covariant spinor derivative, $e_\a^{\ \mu} $ is the tetrad (first index: flat Minkowski, second index: curved space-time). Here the tetrad is defined by

\begin{equation*}
e_{\alpha}^{\mu} g_{\mu \nu}e_{\beta}^{\nu}=\eta_{\alpha \beta},
\end{equation*}

where $g_{\mu \nu}$ denotes the metric tensor and $\eta_{\alpha \beta}$ is the
Minkowski metric. In 2D the tetrad can be computed directly from the metric. Now defining the covariant derivative as $D_{\mu} \Psi=\del_{\mu} \Psi + \G_{\mu} \Psi$, where $\G_{\mu}$ denotes the spin connection matrices given by

\begin{align*}
\Gamma_\mu = - \frac{i}{4} \w_\mu^{\a\b} \s_{\a\b},
\end{align*}

where $\s_{\a\b} = \frac{i}{2} [\g_\a,\g_\b]$ and

\begin{align*}
\w_\mu^{\a\b}=\quad&e_\nu^\a \cov_\mu e^{\nu \b}
\end{align*}

The Dirac equation in curved space describes the relativistic Dirac particles (e.g. electrons ) moving on arbitrary manifold trajectories, see Appendix <ref>.

§.§ Theory of strained graphene

To model the single layer carbon atom honeycomb lattice structure we start from the tight binding Hamiltonian which is constructed assuming superposition of local wavefunctions for isolated atoms on a honeycomb lattice <cit.>. In the low energy limit it has been shown that the tight binding Hamiltonian converges to the Dirac Hamiltonian in the continuum limit,

\begin{equation}
\label{eq:dirac_hamiltonian}
H_D=-i v_f \int \Psi^\dagger \gamma^0 \gamma^i \del_i \Psi d^2x,
\end{equation}

in natural units, where $\Psi$ is in the chiral representation, $v_f$ is the Fermi velocity. In the context of graphene, the general Dirac spinor is defined as $\Psi=(\Psi_a^K,\Psi_a^{K'})=(\psi_A^K,\psi_B^{K'},\psi_A^{K'},\psi_B^{K})$, for sub-lattices $A,B$ and valleys $K,K'$.

The convergence from the tight binding Hamiltonian to  Eq. (<ref>) can be seen as the Dirac cones in graphene with linear dispersion relation at the conduction and valence band connecting point $E=p$ for $E$, energy and $p$, momentum.

The equation of motion from this Hamiltonian is simply the  Dirac equation.
In this work, we consider a static space-time metric with trivial time components

\begin{align*}
g_{\mu\nu} =
\begin{pmatrix}
1 & 0 \\
0 & -g_{ij}
\end{pmatrix},
\end{align*}

where the latin indices run over the spatial dimensions. This simplifies the Dirac equation Eq. (<ref>)

\begin{equation}
\del_t \Psi + \s^a e_a^{~i}(\del_i + \G_i)\Psi = 0-i\g^0 m \Psi,
\end{equation}

with $\s^a = \g^0 \g^a$. After addition of external vector and scalar potentials $A_i(x)$ and $V(x)$ respectively as explained in Ref. <cit.>, the Dirac equation takes the following form:

\begin{equation}
\label{eq:Diracfull}
\del_t \Psi + \s^a e_a^{~i}(\del_i + \G_i- i A_i)\Psi =-i\g^0 (m-V) \Psi.
\end{equation}

Defining the Dirac current with $J^\mu = \overline \Psi \g^\mu \Psi$, the conservation law for can be written as $\del_t \rho + \nabla_i J^i = 0$,
where $\r = \Psi^\dagger \Psi \in \R$ and the $J^i = \overline{\Psi}\gamma^i \Psi \in \R$.

The standard Dirac Hamiltonian for Eq. (<ref> equation is

\begin{equation}
\label{eq:hamiltoniandirac}
H_D=-i \int \Psi^\dagger \sigma^a e_a^{~i}( \del_i + \G_i-i A_i)\Psi \sqrt{g} d^2x,
\end{equation}

For graphene the effective Hamiltonian looks like<cit.>:

\begin{equation}
\label{eq:hamiltoniangraphene}
H^*_D=-i v_f \int \Psi^\dagger \sigma^a (v_a^{*i} \del_i + \G_a^*-i A^*_a)\Psi d^2x,
\end{equation}

where $v_a^{* i}=\d_{a i} + u_{a i} -\beta \e_{a i}$ is the space depended Fermi velocity, $\G^*_a=\frac{1}{2 v_f}\del_j v_a^{* j}$ is a complex gauge vector field which guarantees the hermiticity of the Hamiltonian and $A^*_a$ is a strain-induced pseudo-vector potential given by $A^*_a=(A_x^*,A_y^*)=\frac{\beta}{2a}(\e_{xx}-\e_{yy},-2\e_{xy}$), $\beta$ is the material dependent electron Grueneisen parameter, $a$ the lattice spacing and $\e_{i\jmath}= u_{i\jmath} +\frac{1}{2}\del_i h \del_j h$ the general strain tensor with in-plane, $u_{i\jmath}$ and out of plane, $h$ deformations.The term $u_{a i}$ in $v_a^{* i}$ can be interpreted as the deformation potential term and is purely a geometric consequence due to lattice distortion, it does not depend on the material as long as it has the same topology. Comparing this to the standard Dirac Hamiltonian in curved space Eq. (<ref>) we can match both Hamiltonians $H_D$ and $H^*_D$ by fulfilling the following relations:

\begin{equation}
\label{eq:effectivefields}
v_a^{*i}= v_f \sqrt{g}e_a^{~i}, \ \ \G_a^*= v_f \sqrt{g}e_a^{~i}\G_i, \ \ A^*_a= v_f \sqrt{g}e_a^{~i}A_i.
\end{equation}

All three can be simultaneously fulfilled by an effective metric tensor derived from the explicit expression of the tetrad <cit.>. The effective Dirac model for non-uniformly strained graphene as explained in Ref. <cit.> does not distinguish between the graphene valleys K and K' as it relies on the basic principle that the theory for graphene under nonuniform strain should describe the particular case of a uniform strain, where both Dirac points in the  Brillouin zone are affected the same way and thus considered equivalent. The numerical solutions are obtained with the Quantum Lattice Boltzmann Method as described in Appendix <ref> and Ref. <cit.>.

§ MÖBIUS GRAPHENE STRIP

The system geometry is initialized to the Möbius ribbon by the discreet mapping (or chart),

\begin{align}
h^\alpha(\theta,r)	= & \begin{pmatrix}
(R+ w r \cos(\eta\theta/2))cos(\theta) \\
(R+ w r \cos(\eta\theta/2))sin(\theta) \\
w r \sin(\eta\theta/2)
\end{pmatrix}
\end{align}

with $\theta \in \{-\pi, \pi \}$, $ r \in \{-L_r/2, L_r/2 \}$,  half-width $w$, mid-circle of radius $R$, number of turns $\eta$ and at height $z=0$. $L_r$ is the domain size across the radial direction, for simplicity we set $L_r/2=1$. In the simulations, we consider a square sheet with reverse periodic and closed boundary conditions on a grid of size $l_\theta \times l_r= 512 \times 512$ or $100nm \times 100nm$, $A_a$, the external potential is set to zero.

The metric tensor can be computed from the $h^\alpha(x,y) $ relating the positions of the atoms from the three dimensional flat space (laboratory frame with Minkowski-metric) to the curved space by:

\begin{equation}
g_{ij}= \frac{\del h^\alpha(\theta,r)}{\del x^i}\frac{\del h^\beta(\theta,r)}{\del x^j}\eta_{\alpha \beta},
\end{equation}

With this parametrisation the metric is cast to a diagonal form,

\begin{align*}
g_{i j} =
\begin{pmatrix}
A^2 & 0 \\
0 & B^2
\end{pmatrix}.
\end{align*}

where $A=\sqrt{R+wr\cos(\eta\theta / 2)+w^2r^2/4}$ and $B=w$. The reversed periodic boundary at $\theta=\pi$ is implemented as $\Psi(\theta=-\pi,r)=\Psi(\theta=\pi,-r)$. From Ref. <cit.> a valid current conserving and time-reversal symmetric boundary condition for the Dirac equation can be written as:

\begin{equation}
\Psi=(\mathbf{v} \cdot \mathbf{\tau}) \otimes (\mathbf{n} \cdot \mathbf{\sigma}), ~ \mathbf{n} \bot \mathbf{n}_B,
\end{equation}

where $\mathbf{n}_B$ is the unit vector in the $r-\theta$ plane normal to the boundary, $\mathbf{v}$ and $\mathbf{n}$ are three dimensional unit vectors  and $\mathbf{\tau}$, $\mathbf{\sigma}$ are Pauli matrices. In the case of a mass-less  particle inside a one-dimensional box, assuming perpendicular reflection, it is sufficient to set $\Psi_{1,2}(r=-1)=\Psi_{1,2}(r=1)=0$. This closed boundary condition ensures probability conservation for the decoupled Weyl equations such that $\rho(r=-1)=\rho(r=1)$ and $J^\mu(r=-1)=J^\mu(r=1)=0$. Numerically, for the present simulational timescales, there is some small $(< 1\%)$, resolution convergent error related to finite size effect of the wavefunction.

The metric tensor, although diagonal, for the typical choice of $w\sim 1$ imposes some large gradients on the spin connection $\Gamma^i$ which introduce numerical instabilities. From the form of the metric it can be shown that for $w \ll 1$, the metric variation is minimized resulting in more stability, physically this results in a long and thin strip. It should be noted that this choice of parametrization results in non-zero curvature across the whole domain, seen from the Ricci scalar, Appendix <ref>. Additionally, the choice of $\theta$ is intentional such that the spin connection is continues across the reverse periodic boundary, see Fig. <ref>.  This is not the case for the most common convention $\theta \in \{0, 2\pi \}$. The spin connection (or the Christoffel symbols) is not a gauge independent quantity but the metric and Ricci tensors are. In fact, the complete continuity of the spin connection is only achieved for $w \ll 1$, for $w \sim 1$ there is still some discrepancy as $\omega_i^{jk}$ is not symmetric in $r$.

Correspondingly to the real space curvature, the Berry connection is not gauge independent but the Berry phase and the Berry curvature are. The similarities of the Ricci and Berry curvatures are further investigated in Appendix <ref>, where we show that a Gaussian wave-packet completing a circle around a Möbius band will attain a phase $\pi$ equivalent to a Berry phase. This is a consequence of the topology and it represents a specific illustration of the relation between Berry and real space curvature. Furthermore, in Appendix <ref> we present an example of a periodic, quantized alternating current device where a spacial equivalent of Bloch oscillations are observed. The observed oscillations are validated analytically by a set of semi-classical equations of motion dependent on Berry curvature.

The non-trivial spin connection component on the Möbius strip for half-width $w=0.1$ and mid-circle of radius $R=1$.

To realize the energy levels of a Dirac particle on a Möbius graphene ribbon we initialize the wavefunction to a plane wave solution,

\begin{align}
\psi(\mathbf{k,\mathbf{r}})=  & \begin{pmatrix}
i \sin(k_r r) \\
0 \\
0 \\
\frac{\hbar c \mathbf{k} }{E} \cos(k_r r)
\end{pmatrix}
\end{align}

where $\hbar=c=1$, $E=|k|$, $\mathbf{r}=(\theta, r)$ and  $\mathbf{k}=(k_\theta,k_r),~k_\theta=0$. The wave-functions at some later time-step are shown in Fig. <ref>.

The ground (a), 2$^{nd}$ (b) and 3$^{rd}$ (c) state wavefunctions at $t=1000$ on the Möbius graphene strip for half-width $w=0.01$ and mid-circle of radius $R=1$.

To study the topological properties and observe a Hall-type effect it is necessary to realize a forward moving wavefunction that explores the complete domain. An electric field is both experimentally challenging and theoretically inconsistent across a periodic system such as the Möbius graphene strip. Additionally, a magnetic field is difficult to keep tangential to such a manifold and unnecessary, as shown later. Consequently a wavefunction with non-zero $k_\theta$ is initialized in the domain for a particle and a hole respectively,

\begin{align}
\psi_P(\mathbf{k},\mathbf{r})=  & \begin{pmatrix}
i \sin(k_\theta \theta) \\
0 \\
0 \\
\cos(k_\theta \theta)
\end{pmatrix}
\psi_H(\mathbf{k},\mathbf{r})=  & \begin{pmatrix}
0 \\
\cos(k_\theta \theta) \\
i \sin(k_\theta \theta)\\
\end{pmatrix}.
\end{align}

In the course of this section the subscript $\theta$ is dropped for brevity, $k_\theta \rightarrow k$. The time evolution of the radial space-integrated Dirac current $\bar{J}^r$  for $\psi_P$ is shown in  Fig. <ref>(a); for the zero-momentum and $k=2\pi/l_\theta$ cases, both for the particle and hole wave-functions. The two dashed lines lie on opposite sides relative to the $k_\theta=0$ scenario, indicating a net anomalous velocity effect (in the $r$-direction) as a consequence of a $\theta$ velocity, i.e. a Hall current. The asymmetry ratio, $A_P/A_H$, is plotted in the inset of Fig. <ref>(a) for two different half-widths $w$, defined by

\begin{equation}
\frac{A_P}{A_H}=\frac{\bar{J}^r_{P}-\bar{J}^{r}_{P,k=0}}{\bar{J}^{r}_{H}-\bar{J}^{r}_{H,k=0}},
\end{equation}

where $\bar{J}^r_{P,H}$ denote the particle, hole currents for $k=2\pi/l_\theta$ and $\bar{J}^r_{P,H,k=0}$ denote the particle, hole currents for  $k=0$.  In practice  $\bar{J}^r_{P,k=0}$ and $\bar{J}^r_{H,k=0}$ are equal. The dependence of the asymmetry ratio on $w$ suggests that the difference between $\bar{J}^r_{P}$ and $\bar{J}^r_{H}$ is a consequence of the manifold.

The time evolution of $\bar{J}^r_{P}$ and  $\bar{J}^r_{P,k=0}$ for  longer time  are plotted in the inset of Fig. <ref>(b). The oscillations are a result of the geometry and the closed boundary condition as they are also present for $\bar{J}^r_{P,k=0}$. There is no obvious non-zero average $\bar{J}^r_{P}$, or quantum Hall current. Nevertheless,  the non-zero asymmetry ratio and its dependence on $w$ imply again a net anomalous velocity effect. In Fig. <ref>(b) the time evolution of the average particle-hole space-integrated current $\langle \bar{J}^r_{PH} \rangle= |\bar{J}^r_P-\bar{J}^r_H|$, is plotted for two different non-zero $k$. The result clearly shows a quantum spin-Hall current (or particle- antiparticle current), which is $k_\theta$ dependent. The effect is a transient result as there is no constant electric field in order to measure a stationary state solution. Since a Möbius band has only one boundary, it is easier to visualize the result when projected to a rectangle (simulational domain), where the current is flowing towards one direction.

(a) The time evolution of the radial space-integrated Dirac current  $\bar{J}^r$ the zero-momentum particle wavefunction and non-zero momentum particle and antiparticle wave-functions.  The inset shows the asymmetry ratio for two different half-widths. (b)  The time evolution of the average particle-hole space-integrated current for zero and two different non-zero $k_\theta$. The inset shows the time evolution of $\bar{J}^r$  for a particle wavefunction for zero and non-zero momentum.

The spin-Hall current is a consequence of the curvature and topology of the system. The gauge field breaking the symmetry of the system is the spin connection matrix $\Gamma_i$, similar to a twisting magnetic field. For graphene, these curvature effects are realized as forces due to the pseudo-potential in a strained sheet <cit.>. The correspondence between a magnetic field and curvature is also evident from Appendix <ref>, as the Berry phase is directly related to the Aharanov-Bohm effect.

The space integrated spin-Hall current is then a result of the topology, i.e. the twist. This can be understoond theoretically as a space-averaged current would not be affected by the 'sharpness' of the twist and any 'flat' regions are irrelevant. Moreover, assuming geometrical disorder is not greater than the cumulative effect of the twist,  the bulk current will not be affected more than a perturbative correction. Furthermore,  it is established that the quantum levels on the Möbius strip are quantized in $k_\theta$, see Ref. <cit.>.  The current measured here is then quantized with respect to $k_\theta$. Therefore, the response is identified as a quantum spin-Hall current, which can be used to define a Hall conductivity.

§ THE TORUS

To investigate the relevance of topology and curvature, a similar but simpler shape is also simulated. The system is initialized to the torus geometry by the discreet mapping,

\begin{align}
h	=  & \begin{pmatrix}
(R+ w  \cos(\phi)) cos(\theta) \\
(R+ w  \cos(\phi))sin(\theta) \\
w  \sin(\theta)
\end{pmatrix},
\end{align}

with $\theta,\phi \in \{-\pi, \pi \}$,  width $w$, mid-circle of radius $R$,
$R,w \in \mathbb{R}_{>0}$, both boundaries are periodic. In this parametrization the metric is diagonal and $\phi$ dependent $g_{ij}=\mathrm{diag}((R+w\cos(\phi))^2,w^2)$. In  Fig. <ref> the time-evolution of the azimuthal current $\bar{J}^\phi$ is plotted for a particle, (dotted line) and hole (dashed line) wavefunctions as in Sec. <ref>. The solid line represents a particle wavefunction with $k_\theta=0$. Here, it is evident that a Hall current develops in the presence of a $\theta$ velocity. The different topology, i.e. the absence of the twist and closed boundaries, results in a simpler behaviour, where a transverse current develops as a consequence of curvature.

The time-evolution of the azimuthal current $\bar{J}^\phi$ for a non-zero momentum particle, (dotted line) and hole (dashed line) wavefunctions. The solid line represents a particle wavefunction with $k_\theta=0$

§ CONCLUSIONS AND OUTLOOK

We have presented a study of the topological and geometrical transport properties of a Möbius graphene ribbon. The challenges and resolutions in simulating such a system were outlined. The continuity of the spin connection was achieved by a gauge transformation and limiting the half-width of the ribbon. We provide some possible numerical wave-function solutions of the Möbius riboon.

In the absence of a magnetic field, we measure a quantum spin-Hall current on the graphene ribbon originating from topology and curvature, whereas a quantum Hall current is not observed. This result also represents an example of the correspondence between a magnetic field and a curved manifold. The torus geometry is simulated for comparison, where a Hall current is measured. Additionally, a concrete illustration of the equivalence between the Berry and Ricci curvature is presented through a travelling wave-packet around the Möbius band.  Finally we propose a periodic quantized alternating current device with a curved graphene sheet.

Building on these results, higher order and different types of topologies can be further investigated in the context of curvature. Without the need of a magnetic or electric field and by further understanding the properties of topological graphene sheets, simple and topologically robust quantum devices could be developed just by exploiting their geometry.

The authors are grateful for the financial support Swiss National Science Foundation, under Grant No. 200021 165497

§ CURVED-SPACE QUANTUM LATTICE BOLTZMANN

The quantum lattice Bolzmann (QLB) method used for solving the Dirac equation as minimally coupled to curved space is an extension of the original method developed by Succi et al.  <cit.>. The method exploits the conceptual similarities between the Dirac equation and the Boltzmann equation on the lattice.  We present here the QLB method for a three-dimensional manifold, with straight forward usage to lower dimensional systems.

§.§ The Dirac equation

The classical Boltzmann equation for a particle density distribution function $f(x_a, v_a, t)$ is given by

\begin{equation}
\del_t f + v^i \del_{x^i}f=\mathcal{C}[f]-F^a\del_{v^a}f,
\end{equation}

the left-hand side describes the advection of the distribution function, velocity $v^a$, whereas the right-hand side describes the collisions
between particles and the effect of external forces $F^a$. Furthermore, The Dirac equation in curved space in Eq. <ref> can be cast into a kinetic theory form,

\begin{equation}
\label{eq:dirac_qlbm}
\del_t \Psi + \s^a\del_a \Psi = \mathcal{C} \Psi + \mathcal{F} \psi.
\end{equation}

Therefore, similarly to the Boltzmann equation,  the left hand side represents the 'free streaming' step along matrix valued 'velocities' $\s^i$ and the right hand site contains a 'collision term' and a 'Forcing term'.

The collision term of Eq. <ref> is represented by

\begin{equation}
\label{eq:collision}
\mathcal{C}=-(i m \g^0 + \s^a e_a^i\G_i),
\end{equation}

where m is the fermion mass, and the 'forcing term' by

\begin{equation}
\label{eq:forcing_qlbm}
\mathcal{F}=-\s^a(e_a^i-\d_a^i) \del_i.
\end{equation}

The symbols have their usual meaning.
The partial derivative of the Dirac equation is distributed between the streaming part and the forcing term resulting in a lattice compatible streaming operator of the form $\del_t + v^a\del_a$, where $v^a \in \mathbb{Z}$. The forcing term is a consequence of the generalized Dirac matrices $\g^i=e_a^{~i}\g^a$ and captures the bulk of the curvature effects. The partial derivative in Eq. <ref> is approximated by a local finite difference scheme on the lattice.

§.§ Diagonal streaming operator

In order to obtain a diagonal streaming operator the complex $\s$-matrices have to be diagonalized first, which also yields a diagonal velocity matrix with eigenvalues $v^a=\pm 1$. The digitalization is achieved by:

\begin{align*}
X_a^\dagger \,\s^a\, X_a
= \begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & -1 & 0 \\
0 & 0 & 0 & -1
\end{pmatrix}
= \g^0 \qquad \text{for } a=0,1,2,
\end{align*}

with unitary transformation matrices $X_1, X_2, X_3$ given by

\begin{align*}
\setlength{\arraycolsep}{2pt}
\renewcommand{\arraystretch}{0.8}
X_1 &= \T\frac{1}{\sqrt 2}
\begin{pmatrix}
1 & 0 & -1 & 0 \\
0 & 1 & 0 & -1 \\
0 & 1 & 0 & 1 \\
1 & 0 & 1 & 0
\end{pmatrix},
X_2 = \T\frac{1}{\sqrt 2}
\begin{pmatrix}
0 & i & 0 & 1 \\
-i & 0 & i & 0 \\
-1 & 0 & -1 & 0 \\
0 & -1 & 0 & -i
\end{pmatrix},
\\
X_3 &= \T\frac{1}{\sqrt 2}
\begin{pmatrix}
1 & 0 & 0 & -1 \\
0 & 1 & 1 & 0 \\
1 & 0 & 0 & -1 \\
0 & 1 & 1 & 0
\end{pmatrix}.
\end{align*}

The streaming and collision operations are performed in successive steps using operator splitting as simultaneous diagonalization of three $\s$ matrices is not possible:

\begin{align*}
\Psi(t+\frac{\dt}{D}) &=
\exp\big(-\dt\s^1\del_1+\frac{\dt}{D}(\mathcal{C}+\mathcal{F})\big) \Psi(t), \\
\Psi(t+\frac{2\dt}{D}) &=
\exp\big(-\dt\s^2\del_2+\frac{\dt}{D}(\mathcal{C}+\mathcal{F})\big) \Psi(t+\frac{\dt}{D}), \\
\Psi(t+\dt) &=
\exp\big(-\dt\s^3\del_3+\frac{\dt}{D}(\mathcal{C}+\mathcal{F})\big) \Psi(t+\frac{2\dt}{D}),
\end{align*}

where $D=3$ denotes the dimension. Each streaming step can be diagonalized by left multiplying with $X_a^\dagger$.

\begin{align}
\label{eq:diagonalization}
X_a^\dagger	\Psi(t+\frac{\dt}{D}) &=
\exp\big(-\dt\s^a\del_a+\dt(\tilde{\mathcal{C}}_a+\tilde{\mathcal{F}}_a)\big) \tilde{\Psi}_a(t),
\end{align}

with the definitions:

\begin{align*}
\Tilde{\Psi}_a \coloneqq X_a^\dagger \Psi, ~~ \tilde{\mathcal{F}}_a \coloneqq \frac{1}{2}X_a^\dagger \mathcal{F} X_a, ~~ \tilde{\mathcal{C}}_a \coloneqq \frac{1}{2} X_a^\dagger \mathcal{C} X_a,
\end{align*}

for $a=1,2,3$ (no Einstein summation is used here), and the exponential approximated as

\begin{align*}
& \exp\big(-\dt\s^a\del_a+ \dt(\tilde{\mathcal{C}}+\tilde{\mathcal{F}})\big) \\
& \approx \big(\mathbb{1} -\dt\s^a \del_a+\dt(\tilde{\mathcal{C}_a}+\dt\tilde{\mathcal{F}}_a)+(\mathbb{1} - \frac{\dt}{2}\tilde{\mathcal{C}}_a)^{-1}(\mathbb{1} + \frac{\dt}{2}\tilde{\mathcal{C}}_a)\big)
\end{align*}

The expansion of the collision operator $e^{\dt\tilde{\mathcal{C}}_a}$ is unitary and thus conserves exactly the probability of the wavefunction. The streaming $e^{-\dt\g^0\del_a}$ and forcing $e^{\dt\tilde{\mathcal{F}}_a}$ operators are not expanded in analogous as this is prohibited by the derivative. A simple $2^{nd}$-order expansion is performed limiting the probability norm to $\Delta t^2$ accuracy. The operator splitting contains an error of $\mathcal{O}(\dt^2)$, as $e^{\dt X}\cdot e^{\dt Y}=e^{\dt(X+Y)+1/2\dt^2[X,Y]}=e^{\dt(X+Y)}+\mathcal{O}(\dt^2)$.

The manifold is described by a chart h defined in linear space (see Appendix. <ref>) discretized on a regular rectangular lattice. The curved space quantum Lattice Boltzmann method evolves the spinor $\Psi = (\Psi^+, \Psi^-) = (\Psi_1^+,\Psi_2^-,\Psi_1^-,\Psi_2^+)$ from $t$ to $t+\delta t$. Once the operators are split, the following algorithm is performed consecutively for each lattice direction $n_a$, where $n_1=(1,0)$, $n_2=(0,1)$ and $a=1,2$.

* Rotation: The spinor is rotated by $X_a$

\begin{equation}
\tilde{\Psi}_a(x,t)=X^\dagger_a \Psi(x,t).
\end{equation}

* Collisions and curvature: The collision and force operators are applied on the rotated spinor,

\begin{equation*}
\tilde{\Psi}_a^*(x,t)=\big( \Delta t \tilde{\mathcal{F}}_a + ( \mathbb{1} - \frac{\Delta t}{2} \tilde{\mathcal{C}}_a)^{-1}  ( \mathbb{1} + \frac{\Delta t}{2} \tilde{\mathcal{C}}_a)   \big) \tilde{\Psi}_a(x,t),
\end{equation*}

where $\tilde{\Psi}_a^*(x,t)$ denotes an auxiliary field,

\begin{equation}
\label{eq:collisiontilde}
\tilde{\mathcal{C}}_a=\frac{1}{2}X^\dagger_a \mathcal{C}X_a=-\frac{i}{D}m(X^\dagger_a \gamma^0 X_a)- \gamma^0e_a^i\Gamma_i,
\end{equation}

\begin{equation}
\label{eq:forcingtilde}
\tilde{\mathcal{F}}_a\tilde{\Psi}_a(x,t)=\big(e_a^i-\delta_a^i\big) \Big(\tilde{\Psi}_a(x \mp n_i \Delta t,t)-\tilde{\Psi}_a(x,t) \Big),
\end{equation}

$n_i$ the lattice direction and $\mathcal{C}$ the collision term, Eq. (<ref>).
The upper sign applies for the spin-up components $(\Psi_1^+,\Psi_2^+)$ and the lower sign for the spin-down components $(\Psi_1^-,\Psi_2^-)$.

* Streaming: The spinor components are streamed to the closest grid points along the lattice direction $\pm n_a$,

\begin{equation}
\tilde{\Psi}_a(x,t+\frac{\Delta t}{2})=\tilde{\Psi}_a^*(x \mp n_a\Delta t,t).
\end{equation}

* Inverse Rotation: The spinor is rotated back by $X_a$,

\begin{equation}
\Psi_a(x,t+\frac{\Delta t}{2})=X_a \tilde{\Psi}_a(x,t+\frac{\Delta t}{2}).
\end{equation}

* Repeat steps 2-4 for the next spatial direction

The external potentials $V(x)$, scalar, and $A(x)$, vector are added to the collision operator Eq. <ref> such that

\begin{equation}
\tilde{\mathcal{C}}_a=\frac{1}{2}X^\dagger_a \mathcal{C}X_a=-\frac{i}{D}(m-V)(X^\dagger_a \gamma^0 X_a)- \gamma^0e_a^i(\Gamma_i-iA_i).
\end{equation}

The simulation for strained graphene is carried out with modified Eqs. (<ref>,<ref>) such that

\begin{equation*}
\tilde{\mathcal{C}}_a \rightarrow \sqrt{g} \tilde{\mathcal{C}}_a, ~ e_a^i\rightarrow \sqrt{g}e_a^i.
\end{equation*}

The additional factor of $\sqrt{g}$ originates from the volume element of the Hamiltonian Eq. (<ref>).

§ RIEMANNIAN GEOMETRY

The Latin indices run over the spatial dimensions and Einstein summation convection is used for repeated indices.

A $D$ dimensional curved space is represented by a Riemannian manifold M, which is locally described by a smooth diffeomorphism $\mathbf{h}$, called the chart. The set of tangential vectors attached to each point $\mathbf{y}$ on the manifold is called the  tangent space $T_\mathbf{y} M$. In the fluid model, all the vector quantities are represented as elements of $T_\mathbf{y} M$. The derivatives of the chart $\mathbf{h}$ are used to define the standard basis $(\textbf{e}_1,...,\textbf{e}_D)=\frac{\del\mathbf{h}}{\del x^1},...,\frac{\del \mathbf{h}}{\del x^D}$.

The metric tensor $g$ can be used to measure the length of a vector or the angle between two vectors. In local coordinates, the components of the metric tensor are given by

\begin{equation}
g_{ij}(x)= \textbf{e}_i(x)\cdot \textbf{e}_j(x)= \frac{\del \mathbf{h}}{\del x^i} \cdot \frac{\del \mathbf{h}}{\del x^j},
\end{equation}

where $\cdot$ is the standard Euclidean scalar product.

For a given metric tensor, the vector $v=v^i\textbf{e}_i \in T_\mathbf{y} M$ has a norm $||v||_g=\sqrt{v^ig_{ij}v^j}$ and a corresponding dual vector $v^*=v^i\textbf{e}_i \in T^*_\mathbf{y} M$ in the cotangent space, which is spanned by the differential 1-forms $dx^i=g(\textbf{e}_i,\cdot)$. The coefficients $v_i$ of the dual vector are typically denoted by a lower index and are related to the upper-index coefficients $v^i$ by contraction with the metric tensor $v_i = g_{ij}v^j$ or equivalently, $ v^i=g^{ij}v_j$, where $g^{ij}$ denotes the inverse of the metric tensor. The upper-index coefficients $v^i$ of a vector $v$ are typically called  contravariant components, whereas the lower-index coefficients $v_i$ of the dual vectors $v^*$ are known as the covariant components.

A necessary feature for the description of objects moving on the manifold is parallel transport of vectors along the manifold.  The tangent space is equipped with a covariant derivative $\nabla$ (Levi-Civita connection), which connects the tangent spaces at different points on the manifold and thus allows to transport a tangent vector  from one tangent space to the other along a given curve  $\gamma(t)$. The covariant derivative can be viewed  as the orthogonal projection of the Euclidean derivative $\del$ onto the tangent space, such that the tangency of the vectors is preserved during the transport. In local coordinates, the covariant derivative is fully characterized by its connection coefficients $\Gamma^i_{jk}$  (Christoffel symbols), which are defined by the action of the covariant derivative on the basis vector, $\nabla_j \textbf{e}_k= \Gamma^i_{jk}$. In the standard basis, $\textbf{e}_i = \frac{\del \mathbf{h}}{\del x^i}$, the Christoffel symbols are related to the metric by

\begin{equation}
\Gamma^i_{jk}=\frac{1}{2}g^{ij}(\del_j g_{kl} + \del_k g_{jl} - \del_l g_{jk}).
\end{equation}

Acting on a general vector $v=v^i \mathbf{e}_i,$ the covariant derivative becomes:

\begin{equation}
\nabla_k v =(\del_k v^i + \Gamma^i_{kj}v^j)\mathbf{e}_i,
\end{equation}

where the product rule has been applied, using that the covariant derivative acts as a normal derivative on the scalar functions
$v^i$. Extending to tensors of higher rank, for example the second order tensors $T= T^{ij} $,

\begin{equation}
\nabla_k T=( \del_kT^{ij}+ \G ^i _{kl} T^{lj} + \G^j_{kl} T^{il})\mathbf{e}_i \otimes \mathbf{e}_j
\end{equation}

in this work the basis vectors $\mathbf{e}_i $ are generally dropped. Compatibility of the covariant derivative with the metric tensor implies that $\nabla_k g^{ij}=\nabla_k g_{ij} =0$. This property allows us to commute the covariant derivative with the metric tensor for the raising or lowering of tensor indices in derivative expressions.

The motion of the particle can be described by the curve $ \gamma(t)$, which parametrizes the position of the particle at time $t$. The geodesic equation,  $ \nabla_{\dot{\gamma}} \dot{\gamma} =0 $, in local coordinates $\gamma(t)=\gamma^i(t)\mathbf{e}_i$ is defined by

\begin{equation}
\label{eq:geodesic}
\ddot{\gamma}^i + \Gamma_{jk}^i \dot{\gamma^j} \dot{\gamma^k} = 0.
\end{equation}

The geodesic equation can be interpreted as the generalization of Newtons law of inertia to curved space. The solutions of Eq. (<ref>) represent lines of constant kinetic energy on the manifold, i.e. the geodesics.
The Riemann curvature tensor $R$ can be used to measure curvature, or more precisely, it measures curvature-induced change of a tangent vector $v$ when transported along a closed loop.

\begin{equation}
R(\textbf{e}_i,\textbf{e}_j)v=\nabla_i \nabla_j v-\nabla_j \nabla_i v.
\end{equation}

In a local coordinate basis $ { \textbf{e}_i } $, the coefficients of the Riemann curvature tensor are given by

\begin{align}
R^l_{ijk}&= g(R(\textbf{e}_i,\textbf{e}_j)\textbf{e}_k,\textbf{e}_l)
\\
&=\del_j \Gamma^l_{ik} - \del_k \Gamma^l_{ij} + \Gamma^l_{jm} \Gamma^m_{ik} -\Gamma^l_{km} \Gamma^m_{ij}.
\end{align}

Contraction of $R^i_{jkl}$ to a rank 2 and 1 tensor yields  the Ricci-tensor $R_{ij}=R^k_{ikj}$ and the Ricci-scalar $R=g^{ij}R_{ij}$ respectively, which can also be used to quantify curvature.

The gradient is defined as the $\nabla^i f= g^{ij} \del_j f$, the divergence as $\nabla_i v^i= \frac{1}{\sqrt{g}} \del_i (\sqrt{g} v^i) $, and the integration over curved volume as $V=\int_V QdV$, where $dV=\sqrt{g}dx^1...dx^D=:\sqrt{g}d^Dx$ denotes the volume element. $\sqrt{g}$ denotes the square root of the determinant of the metric tensor.

It should be clarified that in the simulations there is no time curvature and $g_{ij}$ denotes the curved space metric.

§ BERRY PHASE ON THE MÖBIUS RIBBON

To solve the Dirac equation, minimally coupled to curvature, Eq. <ref>, with $A_i=\mathbf{0}$ and assuming that the wave-packet has a negligible profile, $\delta \mathbf{r} \rightarrow 0$, the connection component of the covariant derivative can be absorbed into the wavefunction such that

\begin{equation}
\Psi \rightarrow \Psi \exp\bigg(i \int_{\mathbf{r}_c}^{\mathbf{r}_c+\delta \mathbf{r}} \Gamma_i d\mathbf{r} \bigg)
\end{equation}

where $\mathbf{r}_c$ is the center of mass position and $\Gamma_i$ is the spin connection matrix. For a Gaussian wavepacket with spread $\sigma$ and momentum $\mathbf{k}$

\begin{equation}
\Psi(\mathbf{r},\mathbf{k})   = \frac{1}{\sqrt{2 \pi \mathcal{\s}^2}}
\begin{pmatrix}
1 \\
0 \\
0 \\
\end{pmatrix}
e^{i \int \Gamma_i d\mathbf{r}} e^{-\frac{|\mathbf{r}|^2}{4\mathcal{\s}^2}+i\mathbf{k}\cdot\mathbf{r}}.
\end{equation}

This wave-function effectively minimally couples the standard Dirac equation to curved space through the spin connection.
Defining the Berry connection as

\begin{equation}
\label{eq:berry_conection}
\mathcal{A}^i_n(\mathbf{R})=i\langle \Psi(\mathbf{R})| \del_\mathbf{R} | \Psi(\mathbf{R}) \rangle
\end{equation}

for some parameter space $\mathbf{R}$ and eigen-function $n$. The Berry phase can be calculated from the complete loop integral of the connection

\begin{equation}
\gamma=\oint_0^{2 \pi} \mathcal{A}(\mathbf{R})g^{i j}_{\mathbf{R}}d\mathbf{R}.
\end{equation}

In a similar manner to the treatment of the Aharonov-Bohm effect from Berry <cit.>, we define the fast and slow coordinates as $R$ and $r$ respectively such that $\Psi(R) \rightarrow \Psi(r-R)$.
For the Möbius band, choosing the real space as the parameter of integration and  restricting the motion to one-dimension $r \in{0,2\pi}$, the center of the band, the wave-function takes the form

\begin{equation}
\Psi_r(R-r)   = \frac{1}{\sqrt{2 \pi \mathcal{\s}^2}}
\begin{pmatrix}
1 \\
0 \\
0 \\
\end{pmatrix}
e^{i \int \Gamma_r d r} e^{-\frac{|R-r|^2}{4\mathcal{\s}^2}+ik(R-r)}.
\end{equation}

From Eq. <ref> the explicit form of the wavefunction implies that $\mathcal{A}^i=\tr \Gamma_i$. The implication of these result is that the Berry connection and curvatures can be directly related to the real space affine connection and Ricci curvature tensor under some conditions.

As a consequence, the phase change of a wavepacket moving around a Möbius band can be calculated from the Berry phase. Integrating naively around the band

\begin{equation}
\gamma=\oint_0^{2 \pi} \tr \langle \Psi_r| \del_r  \Psi_r \rangle g^{1 1} d\mathbf{r},
\end{equation}

where the $\tr$ denotes the trace of the resulting matrix and it takes into account the spinor nature of the Dirac wavefunction, yields a trivial result: $\gamma=0$. The caveat is that, in this coordinate basis, for both half-width and half-radius equal to unity  $\Gamma_r$  simplifies to a diagonal matrix such that

\begin{equation}
\left(
\begin{array}{cccc}
\frac{1}{2} i \cos \left(\frac{r}{2}\right) & 0 & 0 & 0 \\
0 & -\frac{1}{2} i \cos \left(\frac{r}{2}\right) & 0 & 0 \\
0 & 0 & \frac{1}{2} i \cos \left(\frac{r}{2}\right) & 0 \\
0 & 0 & 0 & -\frac{1}{2} i \cos \left(\frac{r}{2}\right).
\end{array}
\right).
\end{equation}

$\Gamma_r$ is discontinues when $\theta=2\pi \rightarrow 0$. To make the function single valued we can perform a gauge transformation such that $\Psi \rightarrow \Psi'=\Psi e^{i\frac{r}{2}}$. This implies that $\gamma' = \gamma - \oint dr/2=-\pi$. Therefore the wavefucntion picks up a phase of $\pi$ as it moves around the band, a consequence of the specific topology.

§ RICCI SCALAR FOR MÖBIUS BAND

The Ricci scalar for Möbius band is non-zero across the whole domain for the specific choice of parameters

The Ricci scalar of the Möbius band for half-width $w=0.1$ and mid-circle of radius $R=1$.

§ QUANTIZED ALTERNATING CURRENT GRAPHENE STRIP

To further investigate the potential of curvature on curved graphene sheets we propose a periodic system with alternating current (AC) behaviour which is quantized according to its shape.
The system geometry is initialized by the discreet mapping (or chart),

\begin{align}
h^\alpha(x,y)	= & \begin{pmatrix}
x \\
y \\
y \sin(\eta x/2)
\end{pmatrix}
\end{align}

with $x \in \{0, 2\pi \}$, $ r \in \{-L_y/2, L_y/2 \}$,  with $L_y$ being the domain size in the $y$ dimension, see Fig. <ref>. The boundaries are periodic along the $x$-direction and closed at $-L_y/2, L_y/2$.

A Gaussian wave-packet is initialized as

\begin{align}
\Psi(\mathbf{r},\mathbf{k})= \frac{1}{\sqrt{2 \pi \mathcal{\s}^2}} & \begin{pmatrix}
1 \\
\lambda e^{i\theta}
\end{pmatrix}
\end{align}

where $\lambda=\pm 1$ is the band index, $\theta=\arctan(k_y/k_x)$, $\mathcal{\s}$ is a measure for the width, $\mathbf{r}=(x,y)$, $x,~y$ are the two space dimensions, $\mathbf{k}=(k_x,k_y)$, $k_x,~ k_y$ represent the $x$ and $y$ momenta respectively. $k_x$ is initialized to one, $k_y$ to zero and $\lambda$ to plus one. In the simulations, we consider a rectangular sheet with periodic boundary conditions on a grid of size $L_x \times L_y= 256 \times 128$ or $20nm \times 5nm$, $A_a$, the external potential is set to zero. The norm of the wave-function, $\norm{\Psi}$, i.e. the probability density, $\rho$ is plotted in Fig. <ref> for the initial and a later time-step. The domain is initialized to a close boundary in the $y$ direction and periodic across $x$.

Density of the wave-packet at the initial time-step superimposed on the real space shape of the ribbon.

The center of charge density

\begin{equation}
\bar{\rho}=\bigg[ \int^{area} r \cdot \rho dA -\int^{area} r \cdot \rho(t=0) dA \bigg] /\int^{area} \rho dA,
\end{equation}

where $dA=dxdy$, is plotted as a time evolution in Fig. <ref>.  The oscillations can be understood as the geometrical equivalent to the Bloch oscillations and they are a consequence of the sinusoidal, periodic domain. The frequency is quantized relative to $\eta$. For a slowly perturbed Hamiltonian and expanding around the wave-packet center $c$, $H=H_c+\Delta H$, assuming a periodic system that can be described with a Bloch wavefunction, the semi-classical equations of motion are given by <cit.>

\begin{align}
\dot{r}_c &=\frac{\del \e}{\del q} - ( \Omega_{qr}\cdot \dot{r}_c + \Omega_{qq}\cdot \dot{q}_c)-\Omega_{qt} \label{eq:semiclassical1} \\
\dot{q}_c &=\frac{\del \e}{\del r} - ( \Omega_{rq}\cdot \dot{q}_c + \Omega_{rr}\cdot \dot{r}_c)-\Omega_{rt}\label{eq:semiclassical2}
\end{align}

where $r_c, q_c$ are the center of mass position and momentum of the wavepacket, $\e$ is the band energy and  $\Omega_{qr}=(\Omega_{qr})_{\a\b}= \del q_\a \mathcal{A}_{r_\b}- \del r_\b \mathcal{A}_{q_\a} $ is the Berry curvature and $\mathcal{A}_{r}$ the Berry connection. As seen in Appendix  <ref>, the Berry phase, and thus Berry curvature, can be related directly to the spin connection $\G_\mu$ through
$\mathcal{A}^i_{n}(\mathbf{R})=i\langle \Psi(\mathbf{R})| \del_\mathbf{R}|\Psi(\mathbf{R}) \rangle \implies \mathcal{A}_{R_i}=\tr \langle \Psi|\Gamma_i|\Psi \rangle$ for some parameter space $\mathbf{R}$ and eigen-function $n$. The non zero terms of Eqs. (<ref>,<ref>) for the specific geometry   are $  \dot{r}_c=\del \e / \del q\approx v_f, ~ \dot{q}_c=\Omega_{rr} \cdot \dot{r}_c$, which imply:

\begin{equation}
\frac{\del q_\a}{\del r_\beta}= (\Omega_{rr})_{\a\b}.
\end{equation}

For small local wave-packet $\delta q_\a = \int (\Omega_{rr})_{\a\b} dr_\b=\int (\del r_\a \mathcal{A}_{r_\b}- \del r_\b \mathcal{A}_{r_\a} )dr_\b$ and thus $δq_y∝sin(ηx)$  and $δq_x ∝cos(ηx)$.
Therefore, the oscillations can be explained by a real space Berry curvature description of the system as well as the classical geodesic equation of a manifold.

The frequency of these Bloch-like oscillations are quantized according to $η$. Finally, as the center of mass density would be equivalent to a driven oscillating current the system can be implemented as a periodic, quantized AC device.

\begin{figure}[!ht]
\includegraphics[width=\columnwidth, height=\columnwidth]{periodic_main.pdf}
\caption{\label{fig:periodic_main} Time evolution of relative position of the center of charge in the radial direction for geometries with $\eta=1,2,3$.}
\end{figure}

\end{document}