Curvature-induced quantum spin-Hall effect on a M\"obius strip

TL;DR

This paper demonstrates through numerical methods that a monolayer graphene M"obius strip exhibits a curvature-induced spin-Hall effect, highlighting the influence of topology on quantum spin transport.

Contribution

It provides the first numerical evidence of a spin-Hall effect induced solely by the curvature and topology of a graphene M"obius strip.

Findings

Curvature induces a spin-Hall current in graphene M"obius strips.

Despite no Hall current, a spin-Hall effect naturally arises from topology.

Symmetry considerations predict the existence of this effect.

Abstract

The quantum Hall effect has been predicted and discovered in various condensed-matter systems. A promising quantum material for such topological effects is graphene. We report the numerical observation of a curvature-induced spin-Hall effect in a monolayer graphene M\"obius strip. The solution of the Dirac equation on the nontrivial and non-Euclidean manifold reveals that despite the absence of a Hall current, a spin-Hall current is a natural consequence for such a topology, as predicted from symmetry arguments.

Figures

Equations

$$H_D=-i\int {\Psi }^{†}{\sigma }^a{e}_a^{\text{leavevmode}i}({\partial }_i+{\Gamma }_i-i{A}_i)\Psi \sqrt{g}{d}^{2}x,$$

$$H_D^{\ast }=-i{v}_f\int {\Psi }^{†}{\sigma }^a(v_a^{\ast i}{\partial }_i+{\Gamma }_a^{\ast }-i{A}_a^{\ast })\Psi d^{2}x,$$

$$h^{\alpha }(\theta ,r)=$$

$$g_{ij}=\frac{\partial h^{\alpha }(\theta ,r)}{\partial x^i}\frac{\partial h^{\beta }(\theta ,r)}{\partial x^j}{\eta }_{\alpha \beta },$$

$$g_{ij}=\left(\begin{array}{cc}{G}_{11}^2& 0\\ 0& G_{22}^2\end{array}\right).$$

$$\Psi =(\mathbf{v}\cdot \tau )\otimes (\mathbf{n}\cdot \sigma ),\text{leavevmode}\mathbf{n}\perp {\mathbf{n}}_B,$$

$$(i{\gamma }^{\mu }{D}_{\mu })\Psi =0,$$

$$\Psi\rightarrow\Psi\exp\bigg{(}i\int_{\mathbf{r}_{c}}^{\mathbf{r}_{c}+\delta\mathbf{r}}\Gamma_{i}d\mathbf{r}\bigg{)}$$

$$\Psi (\mathbf{r},\mathbf{k})=\frac{1}{\sqrt{2\pi {\sigma }^2}}\left(\begin{array}{c}1\\ 0\\ 0\\ -1\end{array}\right)e^{i\int {\Gamma }_{i}d\mathbf{r}}{e}^{-\frac{|\mathbf{r}{|}^2}{4{\sigma }^2}+i\mathbf{k}\cdot \mathbf{r}}.$$

$${\mathcal{A}}_n^i(\mathbf{R})=i⟨\Psi (\mathbf{R})|{\partial }_{\mathbf{R}}|\Psi (\mathbf{R})⟩$$

$$\gamma ={\oint }_0^{2\pi }\mathcal{A}(\mathbf{R})g_{\mathbf{R}}^{ij}d\mathbf{R}.$$

$${\Psi }_r(R-r)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\left(\begin{array}{c}1\\ 0\\ 0\\ -1\end{array}\right)e^{i\int {\Gamma }_{r}dr}{e}^{-\frac{|R-r{|}^2}{4{\sigma }^2}+ik(R-r)}.$$

$$\gamma ={\oint }_0^{2\pi }\mathrm{Tr}⟨{\Psi }_r|{\partial }_r{\Psi }_r⟩g^{11}d\mathbf{r},$$

$$\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).$$

$${\psi }_P(\mathbf{k},\mathbf{r})=$$

$$\bar{J}^{r}=\bigg{(}\int\limits^{x,y}J^{r}(x,y)\sqrt{g}dxdy\bigg{)}$$

$$\frac{A_P}{A_H}=\frac{{\bar{J}}_P^r-{\bar{J}}_{P,k=0}^r}{{\bar{J}}_H^r-{\bar{J}}_{H,k=0}^r},$$

$${\rho }_{\text{LDOS}}(x,E)=\frac{1}{\pi }\sum _n|{\varphi }_n(x){|}^2\mathrm{Im}\frac{1}{E-E_n-i\delta ϵ}.$$

$$h=$$

$$(i{\gamma }^{\mu }{D}_{\mu }+m)\Psi =0,$$

$${\partial }_{t}f+v^i{\partial }_{x^i}f=\mathcal{C}[f]-F^a{\partial }_{v^a}f,$$

$${\partial }_t\Psi +{\sigma }^a{\partial }_a\Psi =\mathcal{C}\Psi +\mathcal{F}\psi .$$

$$\mathcal{C}=-(im{\gamma }^0+{\sigma }^a{e}_a^i{\Gamma }_i),$$

$$\mathcal{F}=-\sigma^{a}(e_{a}^{i}-\d{_}a^{i})\partial_{i}.$$

$$X_a^{†}{\sigma }^a{X}_a=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& -1\end{array}\right)={\gamma }^0\text{for }a=0,1,2,$$

$$X_1$$

$$X_3$$

$$\Psi (t+\frac{\Delta t}{D})$$

$$\Psi (t+\frac{2\Delta t}{D})$$

$$\Psi (t+\Delta t)$$

$$X_a^{†}\Psi (t+\frac{\Delta t}{D})$$

$${\tilde{\Psi }}_a:=X_a^{†}\Psi ,\text{leavevmode}\text{leavevmode}{\tilde{\mathcal{F}}}_a:=\frac{1}{2}{X}_a^{†}\mathcal{F}{X}_a,\text{leavevmode}\text{leavevmode}{\tilde{\mathcal{C}}}_a:=\frac{1}{2}{X}_a^{†}\mathcal{C}{X}_a,$$

Full text

I Introduction

The discovery of the quantum Hall effect has triggered a surge of research in topological states of matter with a multitude of applications in quantum engineering and condensed matter physics haldane_review ; topo_insulators_1 ; topologicalstatesofmatter . The topological invariants chernnumber of a system provide simplicity in the implementation, as well as, protection from disorder topo_disorter . Naturally, graphene has been extensively investigated in this context, inheriting its remarkable characteristics graphenerev2 ; graphenerev1 ; graphenerev3 into a topologically protected system. For example, 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 graphene_Hall ; graphene_spin_Hall ; graphene_haldane . Nevertheless, a much more promising system would require a breaking of the time invariance in a more intrinsic way without the need of an external field edge-state_withou_hflield ; qhe_without_hfield .

One such system is the simplest non-trivial fiber bundle, the famous Möbius strip calvao ; polymer . Theoretical analysis has already hypothesized that a graphene Möbius strip can exhibit topological insulator properties owning to the zig zag edge states mobius_strip ; ci_qah . In this paper we show that a graphene Möbius strip inherently exhibits a quantum spin-Hall current. The strains induced in curved graphene are modeled via a non-Euclidean space description. The result is alluring, because a pure, curved space formulation of the lower energetic states of graphene would inevitably create a spin-Hall current and not a Hall current as theoretically predicted via symmetry arguments Nontrivial . Furthemore, curvature has immense potential as a novel and straightforward control mechanism in quantum devices molfetta2013 ; Cortijo2007 .

In this paper, we carry out the a numerical study of the quantum spin Hall effect on the curved space solution of a graphene Möbius strip. We simulate Dirac particles and relate them to graphene through the low energy approximation of the band structure. The effect of curvature on transport and the energy levels of confined Dirac particles on such a manifold are investigated. The non-relativistic quantum mechanics on a the Möbius strip has been studied mobius_graphene_QM and non-trivial effects on the Möbius topology have also been proposed for molecular devices mobius_molecular_devices . It has already been shown that in graphene, the quantum Hall effect can be realized with careful strain engineering in the zero field case zerofield_qhf .

To this end numerical studies were curried out with a solver of the Dirac equation in curved space JD_thesis ; cdp . The method is based on the conceptual similarities between the Dirac and Boltzmann equations and is an extension of the quantum lattice Boltzmann method succi_qlbm to curved space. Through our simulations we observe a quantum spin-Hall current in the bulk of a Möbius graphene strip. Additionally, we compare the result to a simpler curved space topology, namely the torus and we propose a simple valid current conserving and time-reversal symmetric boundary condition for the method. Furthemore, a specific illustration of the equivalence between the Berry and Ricci curvature is presented analytically through a traveling wave-packet around a Möbius strip.

II The Dirac equation in curved space

By minimally coupling the Dirac equation to curved space and a vector potential, $A_{i}(x)$, the Dirac Hamiltonian takes the form:

[TABLE]

in natural units such that $\hbar=c=1$ where $\hbar$, is Planck’s constant and $c$, the speed of light, here $m$ is the particle mass, $\mu={0,1,2}$ for two-dimensional space-time. Here $\Psi=(\Psi^{+},\Psi^{-})=(\psi_{1}^{+},\psi_{2}^{-},\psi_{1}^{-},\psi_{2}^{+})\in\mathbbm{C}^{4}$ denotes the spinor, and $\gamma^{\mu}=\gamma^{\alpha}e_{\alpha}^{\ \mu}$ are the generalized $\gamma$-matrices, where $\gamma^{\alpha}\in\mathbbm{C}^{4\times 4}$ are the standard $\gamma$-matrices (in Dirac representation). $e_{\alpha}^{\ \mu}$ is the tetrad (first index: flat Minkowski, second index: curved space-time). Here, the tetrad is defined by $e_{\alpha}^{\mu}g_{\mu\nu}e_{\beta}^{\nu}=\eta_{\alpha\beta}$, where $g_{\mu\nu}$ denotes the metric tensor and $\eta_{\alpha\beta}$ is the Minkowski metric. The tetrad basis is chosen such that the standard Dirac matrices can be utilized with no need to transform to a new coordinate basis. $\Gamma_{\mu}$ denotes the spin connection matrices given by $\Gamma_{\mu}=-\frac{i}{4}\omega_{\mu}^{\alpha\b{}}\sigma_{\alpha\b{}},$ where $\sigma_{\alpha\b{}}=\frac{i}{2}[\gamma_{\alpha},\gamma_{\b{]}}$ and $\omega_{\mu}^{\alpha\b{}}=\quad e_{\nu}^{\alpha}{\nabla}_{\mu}e^{\nu\b{}}$. The Dirac equation in curved space describes quantum relativistic Dirac particles (e.g., electrons ) moving on arbitrary manifold trajectories.

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 graphene_tb1 . In the low energy limit it has been shown that the tight binding Hamiltonian converges to the Dirac Hamiltonian in the continuum limit. Therefore, for graphene the effective Hamiltonian is OLIVALEYVA

[TABLE]

where v_{a}^{*i}=\d{_}{ai}+u_{ai}-\beta\varepsilon_{ai} is the space dependent Fermi velocity, $\Gamma^{*}_{a}=\frac{1}{2v_{f}}\partial_{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}(\varepsilon_{xx}-\varepsilon_{yy},-2\varepsilon_{xy}$), $\beta$ is the material dependent electron Grueneisen parameter, $a$ the lattice spacing and $\varepsilon_{i\jmath}=u_{i\jmath}+\frac{1}{2}\partial_{i}h\partial_{j}h$ the general strain tensor with in-plane, $u_{i\jmath}$ and out of plane, $h$ deformations. The term $u_{ai}$ in $v_{a}^{*i}$ can be interpreted as the deformation potential 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 (1) we can match both Hamiltonians $H_{D}$ and $H^{*}_{D}$. The numerical solutions are obtained with the quantum lattice Boltzmann method as described in Appendix and Ref. debus_qlb .

III Method and Results

III.1 Quantum spin-Hall effect on a Möbius strip

III.1.1 Geometry and boundary conditions

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

[TABLE]

with $\theta\in\{-\pi,\pi\}$, $r\in\{-L_{r}/2,L_{r}/2\}$, half-width $w$, mid-circle of radius $R$, and number of turns $\eta$ 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 boundary conditions in one direction on a grid of size $l_{\theta}\times l_{r}=512\times 512$ or 100 $\times$ 100 nm2, $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:

[TABLE]

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

[TABLE]

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

[TABLE]

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, we observe a small $(<1\%)$, resolution convergent error related to the finite size effect of the wavefunction. The error convergence plot is shown in Appendix C.

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, expressed by the Ricci scalar, see the Appendix. Additionally, $\theta$ is intentionally chosen such that the spin connection is continuous across the reverse periodic boundary, see Fig. 1. 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 achieved only for $w\ll 1$, for $w\sim 1$ there is still some discrepancy as $\omega_{i}^{jk}$ is not symmetric in $r$.

III.1.2 Berry connection and symmetry arguments

Topological currents are traditionally described as being a consequence of the integral of the phase space connection, the Berry phase berry_paper . Analogous 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 here, we show that a Gaussian wave-packet completing a circle around a Möbius strip will attain a phase $\pi$ equivalent to a Berry phase. This is a consequence of the topology and represents a specific illustration of the relation between Berry and real space curvature QAC .

To solve the Dirac equation, minimally coupled to curvature,

[TABLE]

where $D_{\mu}$ is the covariant derivative as $D_{\mu}\Psi=\partial_{\mu}\Psi+\Gamma_{\mu}\Psi$ we assume 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

[TABLE]

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}$

[TABLE]

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

[TABLE]

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

[TABLE]

In a manner similar to the treatment of the Aharonov-Bohm effect from the Berry connection berry_paper , we define the fast and slow coordinates as $R$ and $r$ respectively such that $\Psi(R)\rightarrow\Psi(r-R)$. For the Möbius strip, choosing the real space to be 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

[TABLE]

From Eq. 9 the explicit form of the wavefunction implies that $\mathcal{A}^{i}=\operatorname{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 strip can be calculated from the Berry phase. Integrating naively around the band

[TABLE]

where $\operatorname{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

[TABLE]

$\Gamma_{r}$ is discontinuous when $\theta=2\pi\rightarrow 0$. To make the function single valued we can perform a gauge transformation such that $\Psi\rightarrow\Psi^{\prime}=\Psi e^{i\frac{r}{2}}$. This implies that $\gamma^{\prime}=\gamma-\oint dr/2=-\pi$. Therefore the wavefunction picks up a phase of $\pi$ as it moves around the band. This solutions clarifies the connection between the phase and real space curvatures and emphasizes that a topological current is expected from geometrically non-trivial manifolds.

In the case of low energy graphene, the calculated phase will be equal and opposite between the $k<0$ and $k>0$ states for a wavepacket on the central line. To observe a non-trivial current for a graphene strip we need to depart from the center-line approximation and the edges should be taken into consideration. Furthermore since the phase will oscillate the total Hall currents are expected to be zero, which is not the case for the spin-Hall currents, as explained also in Beugeling et al. Nontrivial . The argument follows from the non-orientability of the Möbius strip implying that any pseudovectorial field can not be defined globally and smoothly on such a surface. The obvious consequence of this is the definition of a Chern number within such topology. The introduction of spin $1/2$ degrees of freedom creates a unique orientable manifold that has an injective map to the original (base) space, i.e. an orientable double cover (ODC). Subjecting this ODC to a Haldane-like flux, by virtue of the geometry alone, two Haldane fluxes of opposite chirality are expected to develop for the particle-hole pair. This would be an intrinsic SO coupling inducing spin-Hall currents; that is, the two edge modes will counter-propagate.

III.1.3 Numerical results

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 a Möbius graphene strip. Additionally, a magnetic field is difficult to keep tangential to such a manifold and unnecessary, as shown later. Consequently, a wave function with non-zero $k_{\theta}$ is initialized in the domain for a particle and a hole respectively,

[TABLE]

where $\hbar=v_{f}=1$, $E=|k|$, $\mathbf{r}=(\theta,r)$ and $\mathbf{k}=(k_{\theta},k_{r}),\leavevmode\nobreak\ k_{r}=0$.

In this section the subscript $\theta$ is dropped for brevity, $k_{\theta}\rightarrow k$. The time evolution of the radial space-integrated Dirac current

[TABLE]

for $\psi_{P}$ is shown in Fig. 2(a); for the zero-momentum and $k=2\pi/l_{\theta}$ cases, for both the particle and hole wave-functions. The dashed and the dotted 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 time evolution of $\bar{J}^{r}_{P}$ and $\bar{J}^{r}_{P,k=0}$ for longer times are plotted in the inset of Fig. 2(a). The oscillations are a result of the geometry and the closed boundary condition because 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.

The asymmetry ratio, $A_{P}/A_{H}$, is plotted in the inset of Fig. 2(b) for two different half-widths $w$, defined by

[TABLE]

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$. $\bar{J}^{r}_{P,k=0}$ and $\bar{J}^{r}_{H,k=0}$ are equal to each other. 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 curvature of the manifold and implies again a net anomalous velocity effect.

In Fig. 2(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 solution. Since a Möbius strip has only one boundary, it is easier to visualize the result when it is projected to a rectangle (simulational domain), where the current is flowing towards one direction.

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 OLIVALEYVA . The correspondence between a magnetic field and curvature is also evident in Sec. III.1.2, as the Berry phase is directly related to the Aharanov-Bohm effect aharanov-bohm .

The space integrated spin-Hall current is then a result of the topology, i.e. the twist. This can be understood theoretically as a space-integrated current would not be affected by the ’sharpness’ of the twist and any ’flat’ regions are irrelevant. Moreover, assuming geometrical disorder is not stronger than the cumulative effect of the twist, the bulk current will not be affected more than a perturbative correction.

Furthermore, it has been established that the quantum levels on the Möbius strip are quantized in $k_{\theta}$ for the Schrödinger equation, see Ref. mobius_graphene_QM . Equivalently this quantisation is expected also for the Dirac equation indeed this is evident from the local density of states (LDOS) as calculated from the energy spectrum $E_{n}$ of the system and its normalized eigenfunctions $\phi_{n}(x)$ according to the following relation:

[TABLE]

$\delta\epsilon\approx 0.02eV$ is the approximate broadening of the energy spectrum peaks. The result is plotted in Fig. 3, where the discrete energy levels of the system are clearly visible. 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 the Hall conductivity. Through symmetry arguments it was proposed that the quantum Hall effect cannot exist on the Möbius strip but the quantum spin-Hall effect should be possible Nontrivial . This result is also confirmed here through this direct simulation.

Experimental realisations of Möbius strips of single crystals of NbSe3 have been reported Tanda2002 and recently a supra-molecular strategy for fabrication was proposed fabrication . Graphene Möbius strips could be realized, by similar methods or in combination with graphene shape engineering tools such as optical forging opticalforging . It can be shown that a pulsed laser beam can forge a graphene sheet into controlled three-dimensional shapes in the nanoscale, therefore two halves of a Möbius ribbon can be forged, and further connected together selfassembly . Additionally, due to the topological nature of the current the result will be robust to any possible temperature fluctuation. $\bar{J}^{r}_{P,H}$ can be measured across $r$ at a fixed $k$.

III.2 Curvature-induced Hall current on the torus

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

[TABLE]

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. 4 the time-evolution of the azimuthal current $\bar{J}^{\phi}$ is plotted for a particle, (dotted line) and hole (dashed line) wavefunctions as before. The solid line represents a particle wavefunction with $k_{\theta}=0$. Here, it is evident that an equivalent Hall current develops for both particles and holes in the presence of a $\theta$ velocity. The different geometry results in a simpler behaviour, where a transverse current develops as a consequence of curvature. The net current, and spin-hall currents would be zero on a graphene torus. This result clarifies that the spin-hall current observed in the previous section is of topological nature, developed from the introduction of the non-trivial topology, the Möbius strip.

IV Conclusions and Outlook

We have presented a study of the topological and geometrical transport properties of a Möbius graphene strip. 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 strip.

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. 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 wave-packet traveling around a Möbius strip.

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.

Acknowledgements.

The authors are grateful for the financial support Swiss National Science Foundation, under Grant No. 200021 165497 and Prof. S. Succi and I. Petrides for helpfull conversations.

Appendix A 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. succi_qlbm . 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.

A.1 The Dirac Equation

The Dirac equation can be naturally extended to curved space described by a metric tensor $g_{\mu\nu}$ with a covariant derivative $D_{\mu}$ as

[TABLE]

where $\gamma^{\mu}$ denote the Dirac matrices.

The classical Boltzmann equation for a particle density distribution function $f(x_{a},v_{a},t)$ is given by:

[TABLE]

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 Eq. 19 can be cast into a kinetic theory form,

[TABLE]

Therefore, similarly to the Boltzmann equation, the left-hand side represents the ’free streaming’ step along matrix valued ’velocities’ $\sigma^{i}$ while the right-hand side contains a ’collision’ and a ’forcing’ term.

The collision term of Equation (21) is represented by:

[TABLE]

where $m$ is the fermion mass. The ’forcing term’ is given by:

[TABLE]

where 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 classical streaming operator of the form $\partial_{t}+v^{a}\partial_{a}$, where $v^{a}\in\mathbb{Z}$. The forcing term is a consequence of the generalized Dirac matrices $\gamma^{i}=e_{a}^{\leavevmode\nobreak\ i}\gamma^{a}$ and captures the bulk of the curvature effects. The partial derivative in Equation (23) is approximated by a local finite difference scheme .

A.2 Diagonal Streaming Operator

In order to obtain a diagonal streaming operator, the complex $\sigma$-matrices have to be diagonalized first, which yields a diagonal velocity matrix with eigenvalues $v^{a}=\pm 1$. The diagonalization is achieved by suitable "rotation matrices":

[TABLE]

where the unitary transformation matrices $X_{1},X_{2},X_{3}$ are given by:

[TABLE]

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

[TABLE]

where $D=3$ denotes the spatial dimensions. Each streaming step can be diagonalized by left multiplying with $X_{a}^{\dagger}$.

[TABLE]

with the definitions:

[TABLE]

for $a=1,2,3$ (no Einstein summation is used here). The exponential is approximated as:

[TABLE]

The expansion of the collision operator $e^{{\Delta t}\tilde{\mathcal{C}}_{a}}$ is unitary and thus conserves exactly the probability of the wavefunction. The streaming $e^{-{\Delta t}\gamma^{0}\partial_{a}}$ and forcing $e^{{\Delta t}\tilde{\mathcal{F}}_{a}}$ operators are not expanded, 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 implies an error of order $\mathcal{O}({\Delta t}^{2})$, as $e^{{\Delta t}X}\cdot e^{{\Delta t}Y}=e^{{\Delta t}(X+Y)+1/2{\Delta t}^{2}[X,Y]}=e^{{\Delta t}(X+Y)}+\mathcal{O}({\Delta t}^{2})$.

The manifold is described by a chart $h$ defined in linear space, discretized on a regular rectangular lattice. The curved space Quantum Lattice Boltzmann Method evolves the four-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 in sequence 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}$,

[TABLE] 2. 2.

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

[TABLE]

where $\mathbb{I}$ is the identity matrix and $\tilde{\Psi}_{a}^{*}(x,t)$ denotes an auxiliary field,

[TABLE]

[TABLE]

where $n_{i}$ is the lattice direction and $\mathcal{C}$ is the collision term, Equation (22). The upper sign applies to the spin-up components $(\Psi_{1}^{+},\Psi_{2}^{+})$ and the lower sign to the spin-down components $(\Psi_{1}^{-},\Psi_{2}^{-})$. 3. 3.

Streaming: The spinor components are streamed to the closest grid points along the lattice direction $\pm n_{a}$,

[TABLE] 4. 4.

Inverse Rotation: The spinor is rotated back via $X_{a}$,

[TABLE] 5. 5.

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 Equation (26), such that:

[TABLE]

The simulation for strained graphene is carried out with modified Equations (26) and (27), according to the following scheme:

[TABLE]

The additional factor $\sqrt{g}$ originates from the volume element of the Hamiltonian Equation (2).

Appendix B 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{\partial\mathbf{h}}{\partial x^{1}},...,\frac{\partial\mathbf{h}}{\partial 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

[TABLE]

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^{i}g_{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 $\partial$ 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{\partial\mathbf{h}}{\partial x^{i}}$, the Christoffel symbols are related to the metric by

[TABLE]

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

[TABLE]

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}$,

[TABLE]

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

[TABLE]

The geodesic equation can be interpreted as the generalization of Newtons law of inertia to curved space. The solutions of Eq. (35) 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.

[TABLE]

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

[TABLE]

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 $\nabla^{i}f=g^{ij}\partial_{j}f$, the divergence as $\nabla_{i}v^{i}=\frac{1}{\sqrt{g}}\partial_{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^{D}x$ 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.

Appendix C Convergence of the closed Dirac boundary conditions

The convergence of the proposed quantum lattice Boltzmann boundary condition is shown on Figure 5.

Appendix D Ricci scalar for the Möbius strip

The Ricci scalar for the Möbius strip is non-zero across the whole domain for the specific choice of parameters as shown on Fig. 6.