Quantum transport simulation of non-local response in Weyl semimetals

TL;DR

This paper uses numerical simulations to demonstrate non-local surface transport in Weyl semimetals, revealing current flow on opposite surfaces under magnetic fields, and introduces an exact recursive Green's function method applicable to surface transport studies.

Contribution

It introduces an exact recursive Green's function approach to simulate non-local transport in Weyl semimetals, highlighting surface interconnectedness and potential experimental observability.

Findings

Current flows on the opposite surface opposite to the driven surface under magnetic field.

Recursive Green's function method provides exact surface Green's functions.

Surface-to-surface transport effects are significant and observable.

Abstract

We numerically study the non-local transport signature in Weyl semimetal as a test for interconnectedness of the surface states, using a recursive Green's function method. We drive a current using two leads connected on the same surface (top surface) and apply a magnetic field throughout the system, perpendicular to the surface. We find that this results in a current flowing on the other surface in the direction opposite to the direction of the current on the top surface and we comment on the viability of observing such an effect in experiment. The recursive Green's function method we employ is exact and provides us with the Green's functions of the two surfaces as well as their connecting elements, which can be applied also for other numerical simulations where the effect of surface to surface transport is important.

Figures

Equations

$$H(\mathbf{k})=$$

$$\dot{\mathbf{k}}\propto \mathbf{v}\times (B\hat{z}),$$

$$H=\sum _{ij}\sum _{\eta ,{\eta }^{\prime }}{t}_{i\eta ,j{\eta }^{\prime }}{c}_{i\eta }^{†}{c}_{j{\eta }^{\prime }}+h.c.,$$

$$H^{j\to i}$$

$$j_{j\to i}=(-e){\frac{\partial H^{j\to i}(\varphi )}{\partial \varphi }|}_{\varphi =0},$$

$$H^{j\to i}(\varphi )=t_{i\eta ,j{\eta }^{\prime }}{e}^{-i\varphi }{c}_{i\eta }^{†}{c}_{j{\eta }^{\prime }}+h.c..$$

$$⟨j_{j\to i}⟩=(-e)2\mathrm{Im}\left[\sum _{\eta ,{\eta }^{\prime }}{t}_{i\eta ,j{\eta }^{\prime }}⟨c_{i\eta }^{†}(\tau )c_{j{\eta }^{\prime }}(\tau )⟩\right],$$

$$⟨j_{j\to i}⟩=-2e\mathrm{Im}\sum _{\eta {\eta }^{\prime }}{t}_{i\eta ,j{\eta }^{\prime }}\int d\omega \sum _{\lambda }{f}^{\lambda }(\omega ){\chi }_{j{\eta }^{\prime },i\eta }(\omega ),$$

$${\chi }_{j{\eta }^{\prime },i\eta }(\omega )$$

$$=(2\pi {)}^2\sum _{\lambda }{f}^{\lambda }(\omega ){\left(G(\omega )V^{\lambda }{G}^{†}(\omega )\right)}_{i\eta ,j{\eta }^{\prime }},$$

$$G(\omega )={\left(\omega -H-i\Sigma (\omega )\right)}^{-1},$$

$${\Sigma }_{m\eta ,n{\eta }^{\prime }}(\omega )=-i\sum _{\lambda }{K}_{\alpha ,m\eta }^{\lambda \ast }[g^{\lambda }(\omega ){]}_{\alpha \beta }{K}_{\beta ,n{\eta }^{\prime }}^{\lambda }.$$

$$\left(\omega \mathbb{I}-\tilde{H}\right)G(\omega )=\mathbb{I},$$

$$\tilde{H}$$

$$(\omega -h)G_{n,\alpha }=A{G}_{n+1,\alpha }+B{G}_{n-1,\alpha },$$

$$(\omega -h_t)G_{0,\alpha }=\mathbb{I}{\delta }_{\alpha ,0}+A{G}_{1,\alpha },$$

$$(\omega -h_b)G_{N,\alpha }=\mathbb{I}{\delta }_{\alpha ,N}+B{G}_{N-1,\alpha },$$

$$(\omega -h^{(1)})G_{n,\alpha }=A^{(1)}{G}_{n+2,\alpha }+B^{(1)}{G}_{n-2,\alpha },$$

$$(\omega -h_t^{(1)})G_{0,\alpha }=\mathbb{I}{\delta }_{\alpha ,0}+A^{(1)}{G}_{2,\alpha },$$

$$(\omega -h_b^{(1)})G_{N,\alpha }=\mathbb{I}{\delta }_{\alpha ,N}+B^{(1)}{G}_{N-2,\alpha },$$

$$A^{(1)}=A^{(0)}(\omega -h^{(0)}{)}^{-1}{A}^{(0)},$$

$$B^{(1)}=B^{(0)}(\omega -h^{(0)}{)}^{-1}{B}^{(0)},$$

$$h^{(1)}=h^{(0)}+B^{(0)}(\omega -h^{(0)})A^{(0)},$$

$$h_t^{(1)}=h_t^{(0)}+A^{(0)}(\omega -h^{(0)}{)}^{-1}{B}^{(0)},$$

$$h_b^{(1)}=h_b^{(0)}+B^{(0)}(\omega -h^{(0)}{)}^{-1}{A}^{(0)}.$$

$$(\omega -h_t^{(k)})G_{0,0}=\mathbb{I}+A^{(k)}{G}_{N,0},$$

$$(\omega -h_b^{(k)})G_{N,0}=B^{(k)}{G}_{0,0},$$

$$(\omega -h_b^{(k)})G_{N,N}=\mathbb{I}+B^{(k)}{G}_{0,N},$$

$$(\omega -h_t^{(k)})G_{0,N}=A^{(k)}{G}_{N,N}.$$

$$\left(\omega \mathbb{I}-h^{\prime }\right)G^{\prime }=\mathbb{I},$$

$$\displaystyle h^{\prime}=\left(\begin{array}[]{cc}h_{t}^{(k)}&A^{(k)}\\ B^{(k)}&h_{b}^{(k)}\\ \end{array}\right),\quad G^{\prime}=\left(\begin{array}[]{cc}G_{00}&G_{0N}\\ G_{N0}&G_{NN}\\ \end{array}\right).$$

Full text

I Introduction

Weyl semimetals (WSM) Vishwanath2011 ; Burkov2011a ; Burkov2011b ; Zyuzin2012a ; Hosur2012 ; rmp are three dimensional topological systems with gap-less points in the spectrum near Fermi energy, called Weyl nodes. Quasi-particle excitations near these momentum points follow a relativistic Weyl equation and behaves like Weyl fermions. Such Weyl nodes always appear in pairs of opposite chirality Nielsen1981 . The low-energy excitations of the system are expected to show a number of novel behaviors: such as extremely high mobility and negative longitudinal magneto-resistance. Possibility of observing novel physical phenomena as well as the promise WSMs present for futuristic applications created an enormous excitement in the condensed matter and material physics community for the last few years, yielding a number of achievements in terms of theoretical studies Vazifeh2013 ; Turner2013 ; Son2013 ; Biswas2013 ; Hosur2013 ; Burkov2014 ; Gorbar2014 ; Uchida2014 ; Khanna2014 ; Ominato2014 ; Sbierski2014 ; Burkov2015a ; Burkov2015b ; Goswami2015 ; nonlocal1 ; Khanna2016 ; Behrends2016 ; Rao2016 ; Baireuther2016a ; Tao2016 ; Marra2016 ; Li2016 ; Baireuther2016b ; Madsen2016 ; Obrien2017 ; Khanna2017 ; Bovenzi2017 ; strainWSM and experimental verification WSMARPES ; Xu2015a ; Xu2015b ; Lv2015a ; Lv2015b ; Lu2015 ; Jia2016 ; aag .

A novel aspect of the physics of WSM is the non-local electro-magnetic response in bulk and surface transport nonlocal1 ; nonlocal2 ; nonlocal3 . The surface states of Weyl semimetals, often called Fermi arc, connect the two Weyl nodes of opposite chiralities FA . In presence of a magnetic field, Landau levels are formed with dispersion along the direction parallel to field. The low energy bulk Landau levels and Fermi arcs intermix to form closed orbits, which have been theoretically studied paralllel and experimentally verified expt . Potter et al. potter showed that these orbits are responsible for a ‘conveyor belt’ motion of electrons in presence of magnetic field and applied bias, imparting a non local nature to the currents at top and bottom surface. A similar ‘wormhole tunneling’ between Fermi arcs in momentum space was described by Wang et al. warmhole to justify the plateaus in quantum hall conductivity at the top surface, where the Weyl nodes act as wormholes in the momentum space, transporting electrons from one surface to another. In this work we aim to verify this phenomena using quantum transport calculations of surface currents, with a focus on modifying the recursive Green’s function formalism capable of taking into account the connection between surface states through the bulk. Our setup, as shown in Fig. 1, consists of a Weyl semimetal block and two leads, connected to the surface of the Weyl semimetal, which can drive a surface current $I$ on the top surface. Because the system is metallic current will flow throughout the bulk as well as the bottom surface. Now, when a magnetic field is applied in a direction as shown in the figure, we expect a reversal in the direction of the current on the bottom surface potter , which can serve as a signature of the non-locality of the Fermi arcs, allowing the electrons to tunnel from one surface to another through the bulk of the system. Our simulation confirms this prediction.

The paper is organized as following: in the section  II, we introduce the model for our set up and we discuss basic bulk and surface properties of the system, connecting to the origin of the expected non-local transport signature. Next, in section III and IV, we set up the numerical method, including the method of recursive Green’s function. We discuss our numerical results in section V followed by a discussion.

II Model of WSM

We consider the following Hamiltonian of a minimal model minmodel for a Weyl semimetal in a cubic lattice:

[TABLE]

Here $\sigma_{\alpha}$’s are Pauli matrices acting in the orbital space. For $\gamma<2t_{x}$ the bulk Hamiltonian is gap-less at two momentum points $(\pm k_{0},0,0)$. Near these two momentum points the low-energy Hamiltonian follows Weyl equation with opposite chirality. For simplicity, we keep $\gamma=0$, $t_{x}=t$ ($t$ serves as our unit of energy) and the lattice spacing is taken to be unity. The Hamiltonian is a minimal model of time-reversal broken Weyl semimetal where the time reversal operator is $T=K$, with $K$ being complex conjugation. The low-energy band-structure is shown in Fig. 2.

For a slab of finite thickness in $z$ direction, say of height $L_{z}$, the surface states of the Hamiltonian Eq. (1), i.e, the Fermi arc, connects the Weyl nodes of opposite chiralities, shown in Fig. 2. These surface states have predominant velocity component along the direction perpendicular to line joining weyl nodes and the Fermi arcs on opposite surfaces have velocities in opposite directions.

The interconnectedness of the Fermi arcs is due to their topological origin and can be checked in the following way, which follows Ref. nonlocal1, . Consider the system as shown in Fig. 1. The leads applies a finite bias along $y$ direction. This sets up a current along $y$ direction, with the charge carriers being predominantly from the Fermi arcs. In presence of a magnetic field along the $z$ direction, from semi-classical arguments, we can argue that these electrons will experience a Lorentz force, responsible for a change in lattice momentum of electrons-

[TABLE]

which is a vector tangential to the Fermi arc. Thus the electrons will slide along the Fermi arc, in the momentum space, from one chiral node to another (say, point $a$ to $b$ in Fig. 2(b)). But, due to the finite length of the arc, once the electrons travel to $b$, they must continue to the other surface (point $c$). Now, once they start sliding to point $d$, the change of momentum, following Eq. (2) will acquire the electrons with a velocity opposite to the direction of the current. This process will result in directions of the currents flowing in opposite surfaces to be opposite. We shall test this prediction in our numerical results.

III Formalism for bond-current

In this section we present briefly the Green’s function formalism we employ for the transport simulation. The geometry of our simulation is shown in Fig. 1. The translation symmetry is kept intact in $x$ direction by making a suitable Gauge choice for the vector potential, $\mathbf{A}=(-By,0,0)$. We use a modified form of Landauer-Buttiker formalism to calculate the current at the top and the bottom surface between site $i$ and $i$+1 for every site, where we also sum over the contributions of modes along the transverse direction, marked by quantum number $k_{x}$ (we also consider $\hbar=1$ and the charge of electron $-e=-1$ throughout the text).

For a generic non-interacting Hamiltonian

[TABLE]

where $i,j$ are site indices and $\eta,\eta^{\prime}$ are other orbital indices. For a hopping term in the Hamiltonian,

[TABLE]

the bond current operator is

[TABLE]

where

[TABLE]

On performing the lead average transport , we get the lead averaged bond current-

[TABLE]

which would be independent of the time $\tau$ for a system in steady state. Going to the Fourier space and using Fermi distributions ($f^{\lambda}(\omega)$ for $\lambda^{\rm th}$ lead) of the leads transport , we have

[TABLE]

which is the quantity of interest. The correlation functions are defined as

[TABLE]

where $\lambda$ denotes the left and right leads. $G$ is the full Green’s function of the system, which also includes the self-energy, $\Sigma(\omega)$, from the leads,

[TABLE]

[TABLE]

$K_{\alpha,m\eta}^{\lambda}$ is the hopping amplitude from $\lambda^{th}$ lead’s site $\alpha$ to the system-site $m$ and orbital $\eta$. For simplicity, we assume that the leads have a single orbital and the hopping is equally likely to either orbital of the system. Further, the density of state of leads is taken to be constant throughout the spectra (flat-band limit). Thus, the Green’s functions of leads can be represented by $g^{\lambda}=-i\pi\rho^{\lambda}$, which is independent of $\omega$ and $\rho^{\lambda}$ is the effective density of states at lead $\lambda$. Finally, we define $V^{\lambda}=K^{\lambda\dagger}\rho^{\lambda}K^{\lambda}$. The leads’ chemical potentials and temperatures are encoded in the Fermi functions $f^{\lambda}(\omega)$, and a bias can be applied by using a chemical potential difference between the leads.

IV Recursive Green’s function Method

The aforementioned Green’s function for our setup can be written in a matrix form as,

[TABLE]

with $\tilde{H}=H+i\Sigma_{l}(\omega)+i\Sigma_{r}(\omega)$ where, $\Sigma_{l}$ and $\Sigma_{r}$ are the self energies due to the left and the right leads. We write $\tilde{H}$ for given $k_{x}$ as

[TABLE]

where each of the $h$ are the Hamiltonian for a layer (given by constant $z$) and $A$ (and $B=A^{\dagger}$) is the connecting matrices from one layer to another. The self energy contribution due to leads attached at the top surface is included in $h_{t}$ and since no lead is attached to the bottom surface, $h_{b}=h$. Including the vector potential for the magnetic field by Peierls substitution, each layer’s Hamiltonian, $h$ is identical, except for the top layer ($0^{\rm th}$ layer), as the self-energies in our setup only appear to the sites which are on the top surface. We consider the number of layers in $z$ direction to be $L_{z}=N+1$.

There are many available ways to obtain surface ($G_{00}$) and bulk Green’s function recursively rec1 ; rec2 ; rec3 ; rec4 ; rec5 , especially, when $N\rightarrow\infty$. As our requirement is to find the bond current not only on the top surface, but as well as on the bottom surface, we require not only $G_{00}$, but also $G_{0N}$ and $G_{N0}$. So, we proceed to modify the recursive Green’s function method in the following way.

On multiplying the above matrices, as in Eq. (12), to have

[TABLE]

where $n=1,2...N-1$ and $\alpha=0,N$. Please note that the identity matrix appearing above is for each layer and not the same as appeared in Eq. (12). Now, on substituting the values of $G_{n+1,\alpha}$ and $G_{n-1,\alpha}$ in terms of $G_{n+2,\alpha}$, $G_{n-2,\alpha}$ and $G_{n,\alpha}$ in the above equations we get,

[TABLE]

where $A^{(1)}$, $B^{(1)}$, $h^{(1)}$, $H_{u}^{(1)}$, $H_{d}^{(1)}$ are identified as,

[TABLE]

where $X^{(i)}$ denotes modified matrices after $i^{\rm th}$ iteration where $0^{\rm th}$ iteration means unmodified matrices. Now, after $k$ iterations, where $L_{z}=N+1=2^{k}+1$ (number of layers along $z$ direction), we get,

[TABLE]

These equations can be written in simplified form

[TABLE]

with

[TABLE]

On solving the above equation, we can find the matrices $G_{0,0}$, $G_{N,N}$, $G_{N,0}$ and $G_{0,N}$, which we use in Eq. (8) to obtain the bond currents in top and the bottom surfaces.

The novelty of the above method is that, effectively, it reduces the computation time for a 2D system (for each $k_{x}$) to an effective 1D system, evaluating both the surfaces’ Green’s functions as well as their connections, where no approximations has been made. Moreover the distance between the surfaces can be controlled highly efficiently, where the separation of the surfaces grow exponentially with increasing iteration.

V Numerical results

Before we present our result with a finite magnetic field, we briefly discuss the $B=0$ case. Due to the time-reversal broken nature of the Hamiltonian, even without an applied bias, the individual bond currents are non-zero (but small) on the surfaces of the WSM for momentum $k_{x}$ between the Weyl nodes (values of $k_{x}$ when the Fermi arc exists), although the net current summing over $k_{x}$ and the layers of $z$ vanishes. This is inherent to the model we use and can be attributed to the anomalous Hall effect Vazifeh2013 ; aHall1 , giving rise to a circulating current along the edges. With an applied bias, such circulating current is lost and the direction of the net current is set by the bias. Further, as electrons from the Fermi arc are chiral on the surfaces (that is, they can only move with negative/positive velocity in top and bottom surface respectively), it is required to apply a bias accordingly to draw a surface current for small bias. This sets the sign of chemical potential of left and right leads. Since we are interested in differential conductance at $\omega=0$, a small chemical potential $\pm e\mu/2$ for left and right leads would suffice.

Now, once we apply the magnetic field, because of the mechanism of Eq. 2, the electrons from the Fermi arc of the top surface slide along the Fermi arc to reach the other surface. This results in a circular motion of the electrons involving the Fermi arcs of both surfaces and the bulk states. Consequently this results in opposite sign of the current flowing on the opposite surfaces. Writing the current as $J=\sigma(0)\delta V$, for a small bias $\delta V$, we show the plot of $\sigma(0)$ (in the unit of $e^{2}/h$) versus $B$ for both surfaces in Fig. 3 and in Fig. 4, where in the later figure we plot the transverse momentum resolution of the current. Fig. 3 clearly shows a reversal of the sign of the current in the bottom surface due to the applied magnetic field, confirming the non-local transport from one surface to another. A further interesting aspect to note that the current on the bottom surface seems to be proportional to the derivative of the top current, which is reminiscent of the relationship between the transverse and longitudinal conductance of quantum Hall systems Halperin , although the origin of such relation is not clear to us for our system. The positions of the initial Weyl nodes in absence of applied magnetic field are marked in Fig. 4, where we see that the Fermi arc states carry most of the currents on the top surface, whereas on the bottom surface the bond current is vanishingly small note1 . With increasing magnetic field, the effective positions of the Weyl nodes change in the momentum space and the net current on the top surface decreases, although the contributions continue to come from the Fermi arcs. Conversely, for the bottom surface, with increasing magnetic field, current continues to increase, in opposite direction to the bias. This is our main result.

We also show in the inset of Fig. 3, how the differential conductance on the two surfaces change with increasing the height of the sample, going as $2^{n}+1$ layers. Interestingly, the $\sigma(0)$ becomes independent of the height of the sample for large sample width, which is precisely the non-local character of the system. This is unrealistic from an experimental perspective and appears in our case because in our setup any momentum scattering mechanism between states of different $k_{x}$ is missing. In realistic situation, because of possible scattering mechanisms, the maximum width when the effect can be observed might be limited.

VI Dicussion

Our results establishes, numerically, that the non-local transport signatures are consistent with the predictions made. Here two issues deserve special attention. First, the model we use breaks time-reversal symmetry, which results in a circulating current through the edge of the system with attached leads of the same bias and without an applied magnetic field. The crucial difference between this circulating current and the non-local transport under a finite bias and an applied magnetic field is that in the later case a net amount of current flows through the bulk of the system, resulting from the mechanism of Eq. 2.

In a time-reversal symmetric Hamiltonian, such circulating current without magnetic field will be absent, but there will be multiple Fermi arcs present on either surfaces. The minimal model for a time-reversal invariant WSM contains two Fermi arcs, where the electrons can move in opposite directions. Similarly, a number of approximation we have made can be lifted for more realistic setup. Namely, we used a constant density of the states of the leads; we did not allow mixing of transverse momentum in presence of the leads and we did not take into account the presence impurities. Whereas the flat-band limit of the leads’ density of states is justified for small applied bias, the other constraints can be removed in our formalism, although we believe, for appropriate range of parameters, the main conclusion of our work to remain unaffected. Consideration of models of multiple Fermi arcs and more realistic leads are left for future studies.

In realistic materials, the expected length of the Fermi arcs are of the order of $0.1a^{-1}$ Huang2015b ; Yang2015b , $a$ being the lattice spacing. In our simulation, we used a larger separation of the Fermi arcs in Fig. 3 and in Fig. 4. For smaller Fermi arc and smaller applied magnetic field, it is required to simulate a much larger system which is numerically difficult. The dependence on the length of the Fermi arc for the current in opposite surfaces has been plotted in Fig. 5, for magnetic field strength of $B=0.005\hbar^{2}/ea^{2}$ (which is of the order of tens of Tesla, considering the lattice spacing is of the order of 5Å), showing that for considerably small length of the Fermi arcs, the differential conductance reduces, but still in the regime where it can be measured experimentally.

VII Summary

To summarize, we study the possibility of non-local transport signature of Weyl semimetal surface-states using a recursive Green’s function method. We find that in presence of an applied magnetic field the surface currents flow in a Weyl semimetal in opposite directions on opposite surfaces due to their connection through the bulk states. Our study agrees with the intuition one derives from semi-classical analysis and provides the ground for further analysis and experimental observation.