Quantized Fermi-arc-mediated transport in Weyl semimetal nanowires

TL;DR

This paper investigates how Fermi-arc surface states influence electrical conductance in Weyl semimetal nanowires, revealing two distinct transport regimes and demonstrating their robustness against disorder.

Contribution

It identifies and characterizes two regimes of transport in Weyl semimetal nanowires, highlighting the role of Fermi-arc surface states and their impact on conductance behavior.

Findings

Surface regime shows quantized conductance steps with flux oscillations.

Bulk-surface regime exhibits quadratic conductance dependence on chemical potential.

Both regimes are robust to disorder effects.

Abstract

We study longitudinal transport in Weyl semimetal nanowires, both in the absence and in the presence of a magnetic flux threading the nanowires. We identify two qualitatively different regimes of transport with respect to the chemical potential in the nanowires. In the "surface regime", for low doping, most of the conductance occurs through the Fermi-arc surface states, and it rises in steps of one quantum of conductance as a function of the chemical potential; furthermore, with varying flux the conductance changes in steps of one quantum of conductance with characteristic Fabry-P\'erot interference oscillations. In the "bulk-surface regime", for highly-doped samples, the dominant contribution to the conductance is quadratic in the chemical potential, and mostly conditioned by the bulk states; the flux dependence shows clearly that both the surface and the bulk states contribute. The two aforementioned regimes prove that the contribution of the Fermi-arc surface states is salient and, therefore, crucial for understanding transport properties of finite-size Weyl semimetal systems. Last but not least, we demonstrate that both regimes are robust to disorder.

Equations

$${\mathcal{H}}_{\mathrm{lat}}=v\left[\sin p_y{\tilde{\sigma }}_x-\sin p_x{\tilde{\sigma }}_y\right]{\tilde{\tau }}_z+v\sin p_z{\tilde{\tau }}_y+$$

$$t\sum _{i=x,y,z}(1-\cos p_i){\tilde{\tau }}_x+v{p}_0{\tilde{\sigma }}_z,$$

$$\mathcal{H}=v{p}_z{\sigma }_z{\tau }_z-v{p}_0{\sigma }_z{\tau }_0+v\left(p_x{\sigma }_x+p_y{\sigma }_y\right){\tau }_z,$$

$${\mathcal{H}}_{\lambda }=\lambda v\left[\left(p_z-\lambda p_0\right){\sigma }_z+p_x{\sigma }_x+p_y{\sigma }_y\right].$$

$${\Psi }_{\lambda }(r,\varphi ,z)=\left(\begin{array}{c}{\rho }_{+}^{\lambda }(r)e^{i(m-1)\varphi }\\ {\rho }_{-}^{\lambda }(r)e^{im\varphi }\end{array}\right)e^{i{p}_{z}z},$$

$${\rho }_{-}^{\lambda }(r)$$

$${\rho }_{+}^{\lambda }(r)$$

$$\displaystyle\left[\rho_{-}^{\lambda}(r)-i\lambda\rho_{+}^{\lambda}(r)\right]\Big{|}_{r=R}=0.$$

$${\mathcal{H}}_{\lambda }=\{\begin{array}{ll}\left(\begin{array}{cc}\lambda (v{p}_z-\lambda \Delta )& \lambda v(p_x-i{p}_y)\\ \lambda v(p_x+i{p}_y)& -\lambda (v{p}_z-\lambda \Delta )\end{array}\right)& \text{for}r<R,\\ \left(\begin{array}{cc}\lambda (v{p}_z-\lambda {\Delta }_{\infty })& \lambda v(p_x-i{p}_y)\\ \lambda v(p_x+i{p}_y)& -\lambda (v{p}_z-\lambda {\Delta }_{\infty })\end{array}\right)& \text{for}r>R,\end{array}$$

$${\Psi }_{\lambda }(r,\varphi )\sim \left(\begin{array}{c}{\rho }_{+}^{\lambda }(r)e^{i(m-1)\varphi }\\ {\rho }_{-}^{\lambda }(r)e^{im\varphi }\end{array}\right)e^{-{\varkappa }_{\lambda }(r-R)},$$

$$\lambda (v{p}_z-\lambda {\Delta }_{\infty }){\rho }_{+}^{\lambda }-i\lambda v\left[{\partial }_r+\frac{m}{r}-\sqrt{(v{p}_z-\lambda \Delta {)}^2-E^2}/v\right]{\rho }_{-}^{\lambda }=E{\rho }_{+}^{\lambda }$$

$$-i\lambda v\left[{\partial }_r-\frac{m-1}{r}-\sqrt{(v{p}_z-\lambda {\Delta }_{\infty }{)}^2-E^2}/v\right]{\rho }_{+}^{\lambda }-\lambda (v{p}_z-\lambda {\Delta }_{\infty }){\rho }_{-}^{\lambda }=E{\rho }_{-}^{\lambda }$$

$${\rho }_{+}^{\lambda }+i\lambda {\rho }_{-}^{\lambda }=0$$

$$i\lambda {\rho }_{+}^{\lambda }-{\rho }_{-}^{\lambda }=0.$$

$${\rho }_{-}^{\lambda }(R)-i\lambda {\rho }_{+}^{\lambda }(R)=0$$

$$\mathcal{H}=\left(\begin{array}{cccc}v(p_z-p_0)& v({\pi }_x-i{\pi }_y)& 0& 0\\ v({\pi }_x+i{\pi }_y)& -v(p_z-p_0)& 0& 0\\ 0& 0& -v(p_z+p_0)& -v({\pi }_x-i{\pi }_y)\\ 0& 0& -v({\pi }_x+i{\pi }_y)& v(p_z+p_0)\end{array}\right)$$

$${\pi }_x\pm i{\pi }_y=-i\hslash e^{\pm i\varphi }\left[{\partial }_r\pm \frac{i}{r}{\partial }_{\varphi }\mp \frac{eB}{2\hslash }r\right],{\pi }_x^2+{\pi }_y^2=-{\hslash }^2\left[{\partial }_r^2+\frac{1}{r}{\partial }_r+\frac{1}{r^2}{\partial }_{\varphi }^2\right]-i\hslash eB{\partial }_{\varphi }+\frac{e^2{B}^2}{4}{r}^2$$

$$\left[{\pi }_x,{\pi }_y\right]=\left[p_x-eBy/2,p_y+eBx/2\right]=\frac{eB}{2}\left[p_x,x\right]-\frac{eB}{2}\left[y,p_y\right]=-i\hslash eB.$$

$${\mathcal{H}}_{\lambda }=\lambda v(p_z-\lambda p_0){\sigma }_z+\lambda v({\pi }_x{\sigma }_x+{\pi }_y{\sigma }_y).$$

$${\mathcal{H}}_{\lambda }=\left(\begin{array}{cc}\lambda v(p_z-\lambda p_0)& -i\lambda \hslash v{e}^{-i\varphi }\left[{\partial }_r-\frac{i}{r}{\partial }_{\varphi }+\frac{eB}{2\hslash }r\right]\\ -i\lambda \hslash v{e}^{+i\varphi }\left[{\partial }_r+\frac{i}{r}{\partial }_{\varphi }-\frac{eB}{2\hslash }r\right]& -\lambda v(p_z-\lambda p_0)\end{array}\right)$$

$$\left(\begin{array}{cc}{v}^2(p_z-\lambda p_0{)}^2+v^2({\pi }_x^2+{\pi }_y^2+\hslash eB)& 0\\ 0& v^2(p_z-\lambda p_0{)}^2+v^2({\pi }_x^2+{\pi }_y^2-\hslash eB)\end{array}\right){\Psi }_{\lambda }=E_{\lambda }^2{\Psi }_{\lambda },{\Psi }_{\lambda }=\left(\begin{array}{c}{\psi }_{+}^{\lambda }\\ {\psi }_{-}^{\lambda }\end{array}\right)$$

$$\left[-{\hslash }^2\left[{\partial }_r^2+\frac{1}{r}{\partial }_r+\frac{1}{r^2}{\partial }_{\varphi }^2\right]-i\hslash eB{\partial }_{\varphi }+\frac{e^2{B}^2}{4}{r}^2+\sigma \hslash eB+(p_z-\lambda p_0{)}^2-E^2/v^2\right]{\psi }_{\sigma }^{\lambda }=0.$$

$$\left[-{\hslash }^2\left[{\partial }_r^2+\frac{1}{r}{\partial }_r-\frac{(m-\Theta (\sigma ){)}^2}{r^2}\right]+\hslash eB(m-\Theta (\sigma ))+\frac{e^2{B}^2}{4}{r}^2+\sigma \hslash eB+(p_z-\lambda p_0{)}^2-E^2/v^2\right]{\rho }_{\sigma }^{\lambda }(r)=0.$$

$$\left[-{\hslash }^2\left[{\partial }_r^2+\frac{1}{r}{\partial }_r-\frac{(m-\Theta (\sigma ){)}^2}{r^2}\right]+(p_z-\lambda p_0{)}^2-E^2/v^2\right]{\rho }_{\sigma }^{\lambda }(r)=0.$$

$$\left(\begin{array}{cc}\lambda (p_z-\lambda p_0)& -i\lambda \hslash \left[{\partial }_r+\frac{m}{r}\right]\\ -i\lambda \hslash \left[{\partial }_r-\frac{m-1}{r}\right]& -\lambda (p_z-\lambda p_0)\end{array}\right)\left(\begin{array}{c}{\rho }_{+}^{\lambda }(r)\\ {\rho }_{-}^{\lambda }(r)\end{array}\right)=\frac{E}{v}\left(\begin{array}{c}{\rho }_{+}^{\lambda }(r)\\ {\rho }_{-}^{\lambda }(r)\end{array}\right).$$

$${\rho }_{-}^{\lambda }(r)$$

$${\rho }_{+}^{\lambda }(r)$$

$$\underset{0}{\overset{R}{\int }}r{J}_m^2\left(\alpha r\right)dr=\frac{R^2}{2}\left[J_m^2\left(\alpha R\right)+J_{m+1}^2\left(\alpha R\right)\right]-\frac{mR}{\alpha }{J}_m\left(\alpha R\right)J_{m+1}\left(\alpha R\right).$$

$$\left\{\xi {\partial }_{\xi }^2+{\partial }_{\xi }+\left[-\frac{\xi }{4}+\frac{E^2/v^2-(p_z-\lambda p_0{)}^2-\sigma \hslash eB}{2\hslash eB}-\frac{m-\Theta (\sigma )}{2}-\frac{{\left(m-\Theta (\sigma )\right)}^2}{4\xi }\right]\right\}{\tilde{\rho }}_{\sigma }^{\lambda }(\xi )=0.$$

$${\tilde{\rho }}_{\sigma }^{\lambda }(\xi )={\tilde{C}}_{m\sigma }{e}^{-\frac{\xi }{2}}{\xi }^{\frac{|m-\Theta (\sigma )|}{2}}{}_1{F}_1\left[\frac{1+m-\Theta (\sigma )+|m-\Theta (\sigma )|}{2}-\frac{E^2/v^2-(p_z-\lambda p_0{)}^2-\sigma eB}{2eB},1+|m-\Theta (\sigma )|,\xi \right],$$

$${\rho }_{\sigma }^{\lambda }(r)=C_{m\sigma }{e}^{-\frac{eB{r}^2}{4}}{r}^{|m-\Theta (\sigma )|}{}_1{F}_1\left[\frac{1+m-\Theta (\sigma )+|m-\Theta (\sigma )|}{2}-\frac{E^2/v^2-(p_z-\lambda p_0{)}^2-\sigma eB}{2eB},1+|m-\Theta (\sigma )|,\frac{eB{r}^2}{2}\right].$$

$${\rho }_{-}^{\lambda }(r)=C_{m-}{e}^{-\frac{eB{r}^2}{4}}{r}^{|m|}{}_1{F}_1\left[\frac{1+m+|m|}{2}-\frac{E^2/v^2-(p_z-\lambda p_0{)}^2+eB}{2eB},1+|m|,\frac{eB{r}^2}{2}\right].$$

$${\Psi }_{\lambda }=\left(\begin{array}{c}{\rho }_{+}^{\lambda }(r)e^{i(m-1)\varphi }\\ {\rho }_{-}^{\lambda }(r)e^{im\varphi }\end{array}\right)e^{i{p}_{z}z}$$

Full text

In 1929 the German mathematician and theoretical physicist Hermann Weyl proposed massless solutions of the Dirac equation Dirac (1928), the so-called “Weyl fermions” Weyl (1929). He demonstrated that in the absence of a mass term, the Dirac equation decoupled into two independent ones, also known as Weyl equations, each describing fermions of a given chirality, right or left. Despite numerous theoretical predictions Herring (1937); Murakami (2007); Wan et al. (2011); Burkov et al. (2011); Weng et al. (2015); Huang et al. (2015), it was not until 2015 that Weyl fermions were first observed as the low-energy excitations in TaAs Xu et al. (2015a) and NbAs Xu et al. (2015b).

Several decades ago Nielsen and Ninomiya showed that for continuous and periodic Hamiltonians with real spectra right and left chiralities are always bound to appear together Nielsen and Ninomiya (1981). Hence, a minimal low-energy model for Weyl fermions must embody both chiralities, and the numbers of particles with right and left chiralities must be equal. The low-energy Weyl fermions of a given chirality $\lambda=\pm 1$ disperse linearly with momentum, $E^{2}_{\lambda}=v^{2}\left[p_{x}^{2}+p_{y}^{2}+(p_{z}-\lambda p_{0})^{2}\right]$, with group velocity $v$ and band touching points $(0,0,\pm p_{0})$ referred to as “Weyl points” or “Weyl nodes”, chosen without loss of generality to be along the $p_{z}$ axis. In three dimensions the density of states of such quasiparticles grows quadratically with energy, while their group velocity is constant, and therefore, the bulk semiclassical conductance of Weyl semimetals $G\propto\mu^{2}$, where $\mu$ is the chemical potential of the sample, versus $G\propto\mu^{3/2}$ in ordinary metals.

Apart from this peculiar bulk property, Weyl semimetals are also known for their surface states — Fermi arcs. The contribution of these states to the transport properties of nanowires made of Weyl semimetals was considered both experimentally in Ref. [Nair and et al., 2018] and theoretically in Refs. [Baireuther et al., 2016; Gorbar et al., 2016; Baireuther et al., 2017; Igarashi and Koshino, 2017; Fu et al., 2018; Deng et al., 2019; Breitkreiz and Brouwer, 2019]. The experimental work mostly focuses on measuring the Shubnikov–de Haas effect, whereas theoretical papers contain semiclassical calculations of the conductance in different regimes.

In this paper, we show that despite being three-dimensional per se, Weyl semimetal nanowires may conduct only through the Fermi-arc surface states. In order to demonstrate the latter, we calculate both analytically and numerically the zero-bias conductance of Weyl semimetal nanowires (see Fig. 1), taking into account the contribution of the Fermi-arc surface states. We focus on studying transport properties of samples in which the transverse dimension of the system $W$ is, on one hand, much larger than the Fermi-arc localization length, but on the other hand, sufficiently small in order to resolve in energy the confinement gap appearing at the band touching points for the bulk states. The conditions above ensure spatial separation of the bulk and the surface states, as well as experimental accessibility of the proposed regime Li et al. (2015); Wang et al. (2016a, b, c, 2018). We show that depending on the chemical potential in the sample there exist two qualitatively different regimes of conductance: “surface regime” and “bulk-surface regime”. In the latter the conductance of the nanowire is conditioned both by the bulk and by the surface states, and it grows quadratically with the chemical potential, showing the expected hallmark of three-dimensional linearly-dispersed electrons. Surprisingly, in the former regime the nanowire shows effective one-dimensional behavior, and the conductance grows in steps of conductance quanta. We explain such a remarkable feature by showing that in finite-size Weyl semimetal systems there is always a window of energies, defined by the bulk confinement gap, where only surface states exist. This inherent feature may serve as a strong evidence of the presence of the Fermi-arc surface states. Furthermore, we study how the conductance is modified by magnetic flux penetrating the wire, and we demonstrate that in the surface regime it changes in steps of conductance quanta with characteristic Fabry-Pérot interference oscillations, whereas in the bulk-surface regime the changes are not quantized. Last but not least, we investigate the effects of weak and strong disorder on our results, and we find both regimes to be robust and our conclusions qualitatively unchanged.

Model.— In what follows we perform transport calculations both numerically and analytically. For the former, we use the following cubic-lattice Hamiltonian adopted from Ref. [Behrends et al., 2019]:

[TABLE]

where Pauli matrices $\bm{\tilde{\sigma}}=\{\tilde{\sigma}_{x},\,\tilde{\sigma}_{y},\,\tilde{\sigma}_{z}\}$ and $\bm{\tilde{\tau}}=\{\tilde{\tau}_{x},\,\tilde{\tau}_{y},\,\tilde{\tau}_{z}\}$ act in spin and orbital subspaces, respectively, $t$ denotes the hopping amplitude and $v$ parameterizes the low-energy velocity of Weyl fermions. We chose the $z$ direction to be the Weyl node separation axis with positions of the nodes given by $p_{z}=\pm p_{0}$. For the sake of brevity we set $\hbar$ and the lattice constant $a$ to unity, restoring them in what follows if needed. Note also that in all numerical simulations we chose $v=1$, $t=2/\sqrt{3}$. This choice of parameters for the Hamiltonian in Eq. (1) yields two Weyl cones of velocity $v$ in the band structure, localized at $\pm p_{0}$.

To perform analytical transport calculations we use a low-energy model with a block-diagonal form in the chirality subspace:

[TABLE]

where $\bm{\sigma}=\{\sigma_{x},\,\sigma_{y},\,\sigma_{z}\}$ and $\bm{\tau}=\{\tau_{x},\,\tau_{y},\,\tau_{z}\}$ denote Pauli matrices acting in orbital and chirality subspaces, respectively. Since the Hamiltonian in Eq. (2) is diagonal in the chirality subspace, for a given chirality $\lambda=\pm 1$ we can write a 2$\times$2 Hamiltonian as follows:

[TABLE]

In order to provide better understanding of transport properties of Weyl semimetal nanowires below we calculate the band structures of wires infinite in the $z$ direction with a finite cross-section in the $x$ and $y$ directions, in the absence of magnetic field. For the lattice Hamiltonian in Eq. (1) it is sufficient to impose zero boundary conditions (also known as open or hard-wall boundary conditions), whereas for the Hamiltonian in Eq. (2) it is necessary to derive boundary conditions, e.g., assuming a large-gap insulator outside of the wire Okugawa and Murakami (2014). We leave the detailed derivation of boundary conditions to Appendix A, presenting here the final result. We consider a cylindrical wire of radius $R$ defined by $x^{2}+y^{2}\leqslant R^{2}$ and we seek the solution of the Schrödinger equation for a given chirality $\mathcal{H}_{\lambda}\Psi_{\lambda}=E_{\lambda}\Psi_{\lambda}$ using the radial symmetry of the problem with the following ansatz

[TABLE]

where $m\in\mathbb{Z}$ denotes the angular momentum quantum number, and $p_{z}$ is the good momentum in the $z$ direction. Radial functions $\rho^{\lambda}_{\pm}$ are defined as follows

[TABLE]

where $\alpha\equiv\sqrt{E_{\lambda}^{2}/v^{2}-(p_{z}-\lambda p_{0})^{2}}$, $J_{m}(\dots)$ is the $m$-th Bessel function of the first kind. The boundary condition thus reads

[TABLE]

The equation above yields the allowed energies for given values of $m$ and $p_{z}$, hence defining the band structure of an infinite Weyl nanowire. Note that since chiralities are decoupled in Eq. (2), the boundary condition does not mix chiralities either.

In Fig. 2 we plot band structures obtained numerically for a square cross-section slab and analytically for a cylinder (see Fig. 1). Nearly flat bands on both panels represent Fermi-arc surface states with characteristic localization length $\ell_{0}\equiv\hbar/2p_{0}$ Okugawa and Murakami (2014). It is worth noting that in the window of energies $\left|E\right|\lesssim 0.175\,t$ only surface solutions exist, whereas bulk solutions are gapped out due to finite-size effects. The existence of such a window is conditioned by the relation between confinement gaps for the surface and the bulk states, namely, $\Delta_{\mathrm{surf}}<\Delta_{\mathrm{bulk}}$. Such a salient discrepancy stems from geometrical factors: indeed, the wave functions of the surface states are confined to a thin layer defined by the circumference of the cross-section of the wire, i.e., $4W$ for a slab and $2\pi R$ for a cylinder. The bulk states, however, have a smaller confinement length, namely, $W$ for a slab or $2R$ for a cylinder. Hence the bulk confinement gap is larger than the surface one: $\Delta_{\mathrm{bulk}}=\hbar v/W$ or $\Delta_{\mathrm{bulk}}=\hbar v/2R$ versus $\Delta_{\mathrm{surf}}=\hbar v/4W$ or $\Delta_{\mathrm{surf}}=\hbar v/2\pi R$, correspondingly. Therefore, we conclude that, remarkably, the window of energies with only surface solutions always exists in finite-size Weyl semimetal systems.

Conductance.— Below we compute the zero-bias conductance of the nanowire both analytically, via the transfer matrix approach Tworzydło et al. (2006); Bardarson et al. (2007); Xypakis and Bardarson (2017), and numerically, using the Kwant package Groth et al. (2014). Thus, we vary the chemical potential $\mu$ and the longitudinal magnetic field $B$, while calculating $G(\mu,B)=\lim_{V\to 0}I/V$, where $V$ is the bias, and $I$ is the current flowing through the wire. The analytical approach was described thoroughly in Refs. [Tworzydło et al., 2006; Bardarson et al., 2007; Xypakis and Bardarson, 2017], thus we leave the details of the calculation to Appendix C. We model the leads attached to the wire using the exact same Hamiltonians given in Eqs. (1) and (2), taken at the chemical potential $\mu_{\infty}$ large enough to emulate metallic electrodes. In practice, it means that in the analytical low-energy model $\mu_{\infty}$ is taken to be larger than any other energy scale, whereas in the lattice model it must be smaller than the bandwidth, and should be chosen to yield the largest possible number of scattering states at $E=\mu_{\infty}$ in the lead. Eventually $\mu_{\infty}$ drops out of all physically meaningful quantities such as, e.g., conductance. We note also that analytically computed curves coincide with those obtained numerically, thus to avoid redundancy we restrict ourselves to presenting here only the numerical data, while leaving the analytical data to Appendix G.

In Fig. 3 we plot the conductance of the wire in the units of the conductance quantum $G_{0}\equiv e^{2}/h$ as a function of the chemical potential $\mu$ in the sample, for fixed values of the magnetic flux penetrating the wire, namely for $\Phi/\Phi_{0}\in\{0,\,1/4,\,1/2,\,1\}$, where $\Phi_{0}\equiv h/e$ is the quantum of flux.

We start by considering the case of zero flux, $\Phi=0$. It is worth noting that there are two qualitatively different regimes of conductance: “surface regime” and “bulk-surface regime”. In the surface regime transport occurs mainly through Fermi-arc surface states, with the characteristic feature being the conductance rising in steps of $G_{0}$ (see the inset in Fig. 3). This peculiar property can be elucidated as follows. The dispersion of the surface states is effectively one-dimensional, i.e., their energy depends only on one of the two good momenta on the surface. Therefore, the transport properties of the surface bands are similar to those of one-dimensional quantum wires (cf. Landauer formula Landauer (1957)), and hence the conductance increases by $G_{0}$ every time the chemical potential crosses a new surface band. The range of chemical potentials for the surface regime is defined by the finite-size gap, which in its turn is $\propto\hbar v/R$, where $R$ is the radius of the wire Zhang and Liu (2018).

Below we turn to the case of $\Phi\neq 0$. In the absence of a Zeeman term, there are two main effects of the applied magnetic field on the band structure. First, the orbital effects lead to the formation of Landau levels Abrikosov (1998), and second, all bands are shifted either down or up, depending on whether the magnetic field is parallel to the axis of the wire or antiparallel, respectively. The reason for the latter is the fact that surface electrons have chiral dispersion, and therefore, their quasiclassical motion at the surface is clockwise or counterclockwise. Since the magnetic field is applied perpendicular to the plane of this motion, it either favors their motion or not, depending on the direction of the field. Thus, the more flux we apply the more we shift the band structure.

Several salient features of the conductance curves in Fig. 3 for nonzero values of flux are worth being discussed. First, it is clear that the shifts in the conductance curves for different values of flux are quantized in the surface regime (modulo interference oscillations) and irregular in the bulk regime. Such difference stems from the fact that, as already mentioned earlier, in the surface regime most of the transport is conditioned by the Fermi arc states, which are localized in a thin layer of width $\ell_{0}$. Such localization ensures that all surface electrons accumulate phases in a coherent manner. Contrary to that, in the bulk-surface regime both the surface and the bulk states are responsible for transport. The transverse parts of the wave functions of the latter are localized at different distances from the axis of the wire, thus making distinct bulk states be affected by different values of the flux. Such inhomogeneous influence of the magnetic field explains “arbitrary” shifts of the conductance curves for varying flux in the bulk-surface regime.

Finally, it is both of theoretical and experimental interest to study how the conductance changes with the applied flux at fixed values of the chemical potential in the wire. Previously, we have identified two qualitatively different regimes of transport depending on the chemical potential. Thus, on the left and right panels in Fig. 4 we plot the flux dependence of the conductance with the chemical potential fixed in the surface and bulk-surface regimes, correspondingly.

First, we analyze the surface regime. As expected from the previous subsection, changes in conductance in that regime occur in steps of the conductance quantum modulo Fabry-Pérot interference oscillations stemming from reflections from the leads. The origin of these oscillations is easy to corroborate: it is sufficient to reduce the length of the wire by a factor of two, and verify that their period doubles; we have checked that this is indeed the case. Second, we turn to the bulk-surface regime. Here, consistent with our antecedent findings, the conductance does not change in regular steps of $G_{0}$. However, since the surface states still contribute to the transport, we can still identify the aforementioned interference oscillations.

Last but not least, we have verified that our results hold in the presence of disorder by modeling the latter as a random uniform onsite variation of the chemical potential with amplitudes lying in $\left[-A_{\mathrm{dis}},\,A_{\mathrm{dis}}\right]$. In Fig. 5 we present conductance curves for disordered samples with disorder amplitudes ranging from $A_{\mathrm{dis}}=0.01$ to $A_{\mathrm{dis}}=0.5$. Average level broadening in the presence of such disorder can be estimated by $\Gamma\approx\frac{1}{3}\pi A^{2}_{\mathrm{dis}}$ Bruus and Flensberg (2004). Thus, for $A_{\mathrm{dis}}=0.25$ the level broadening $\Gamma\approx 0.065\,t$ becomes larger than the level spacing $\Delta E=0.05\,t$. Despite very strong disorder, both aforementioned regimes of conductance remain qualitatively unaffected. The surface regime is robust due to the fact that the dispersion of the Fermi-arc surface states is effectively one-dimensional and chiral, and thus bereft of backscattering.

Discussion and conclusions.— Above we have studied the longitudinal (magneto-)conductance of Weyl semimetal nanowires. First, we have found that depending on the chemical potential in the wire there exist two qualitatively different regimes of transport: surface regime and mixed bulk-surface regime. In the former only the Fermi-arc surface states conduct, giving rise to quantized conductance steps, characteristic for one-dimensional physics. Contrary to that, in the bulk-surface regime both the surface and the bulk states participate in transport, yielding the expected $G\propto\mu^{2}$ dependence. Furthermore, we have investigated how the conductance varies with the magnetic flux penetrating the wire at fixed values of the chemical potential; we have shown that if the chemical potential is tuned to be in the surface regime, then there are regular jumps of one conductance quantum with characteristic interference oscillations, whereas in the mixed bulk-surface regime the conductance changes irregularly with the increasing value of the magnetic flux.

Despite being obtained for a particular model of a Weyl semimetal, our results can be generalized and applied to a broader range of models, both for Weyl and Dirac semimetals (see, e.g., transport experiments in Ref. [Nishihaya and et al., 2019]). First, depending on the symmetries of a given material, realistic Weyl/Dirac semimetals may contain multiple pairs of Weyl cones at low energies, including spin-degenerate ones. Qualitatively, this may lead to higher conductance values, however, further investigations with more realistic models are required to corroborate this hypothesis. Moreover, for a Dirac semimetal with nonzero node separation in the presence of spin the height of the quantized conductance steps will double due to the spin degeneracy of the cones. It is worth discussing also the case of time-reversal invariant Weyl semimetal nanowires. In that case we have two pairs of cones in the Brillouin zone. If both pairs of cones are separated in the $z$ direction, they generate two sets of Fermi-arc surface states related by time-reversal symmetry. An applied magnetic flux shifts some of these states up in energy, whereas their time-reversal counterparts are shifted down. Thus, as a function of the flux the confinement gap closes and reopens periodically, and therefore, in this case the analog of the left panel of Fig. 4 will be periodic in flux, resembling the response of topological insulator nanowires Bardarson et al. (2010).

Our results can be tested experimentally with relatively thin nanowires made of, e.g., TaAs, NbAs, TaP, NbP, GdPtBi, Co3Sn2S2, etc., such that the confinement gap is sufficiently large to be resolved in energy Xu et al. (2015a, b); Lv et al. (2015); Liu et al. (2015); Suzuki et al. (2016); Liu et al. (2018). Note also that realistic wires made of topological semimetals are generally larger than those considered in our work. For instance, in Refs. [Li et al., 2015; Wang et al., 2016a, b, c, 2018] the radii of the nanowires lie in the range of $30$ to $200$ nm, and the Fermi-arc localization length is of the order of $1$ nm. To facilitate numerical simulations, our nanowires were taken to be approximately 10 times smaller than realistic ones, which means that in realistic systems the surface regime will occur in a smaller range of chemical potentials, namely $\mu\lesssim 10\,$meV. While such small energies require higher resolution to be observed in larger samples, our theoretical conclusions will remain qualitatively unchanged. Furthermore, in our geometry one quantum of flux through the wire is equivalent to having a magnetic field of $\approx 50\,$T. In realistic wires one quantum of flux is achieved at smaller values of the magnetic field due to a larger cross-section area.

Appendix A Derivation of boundary conditions

In order to derive boundary conditions for a cylindrical Weyl nanowire we follow Okugawa and Murakami Okugawa and Murakami (2014). The approach is based on modeling the vacuum outside of the wire as an insulator with a gap $\Delta_{\infty}$ larger than all the other energy scales in the system, or formally, $\Delta_{\infty}\to-\infty$. Thus, the system can be described by the following Hamiltonian

[TABLE]

with $R$ being the nanowire radius, and $\lambda$ the given chirality. We also denoted $\Delta\equiv vp_{0}$, where $p_{0}$ controls the positions of the Weyl nodes along the $p_{z}$ axis. Rewriting $p_{x}\pm ip_{y}$ in polar coordinates as $-i\hbar e^{\pm i\phi}\left[\partial_{r}\pm\frac{i}{r}\partial_{\phi}\right]$ and using the radial symmetry of the problem, we employ an exponentially decaying ansatz for the wave function of a given chirality outside of the nanowire:

[TABLE]

where $\varkappa_{\lambda}\equiv\sqrt{(vp_{z}-\lambda\Delta_{\infty})^{2}-E^{2}}/v>0$ since $\Delta_{\infty}\to-\infty$. Inserting the ansatz above into the Schrödinger equation defined by the Hamiltonian in Eq. (8), we get:

[TABLE]

We apply the limit of $\Delta_{\infty}\to-\infty$ to both equations above, and we end up with

[TABLE]

To have a non-trivial solution of the system above, we must satisfy $\rho^{\lambda}_{+}+i\lambda\rho^{\lambda}_{-}=0$, which in turn defines the boundary condition connecting radial parts of the wave function at $r=R$ in the limit of infinitely large $|\Delta_{\infty}|$:

[TABLE]

Appendix B Weyl nanowire Band structure

Below we consider a Weyl wire oriented along the $z$ axis. We use the Hamiltonian in Eq. (2) and we introduce the magnetic field along the $z$ direction (i.e., along the line connecting the Weyl nodes) using the symmetric gauge: $\bm{\mathcal{A}}=B\left(-y,\,x,\,0\right)/2$, therefore $\bm{B}=\operatorname{rot}\mathcal{A}=B\hat{n}_{z}$. Peierls substitution thus gives

[TABLE]

where $\pi_{x}=p_{x}-eBy/2,\pi_{y}=p_{y}+eBx/2$ with the electron charge given by $-e,e>0$. It is also of use to rewrite this Hamiltonian in cylindrical coordinates using

[TABLE]

as well as to compute the commutator

[TABLE]

To solve the Schrödinger equation we first note that problems for different chiralities are independent, therefore, we can rewrite it in a 2$\times$2 simplified form with chiralities $\lambda=\pm 1$:

[TABLE]

Equivalently in cylindrical coordinates:

[TABLE]

We square the Hamiltonian and we get

[TABLE]

For each component of the wave function we have:

[TABLE]

We use the following ansatz $\psi^{\lambda}_{\sigma}=\rho^{\lambda}_{\sigma}(r)e^{i(m-\Theta(\sigma))\phi}e^{ip_{z}z}$ (where $\Theta$ is the Heaviside step function), and we get:

[TABLE]

Below we solve Eq. (14) above in two different cases: $B=0$ and $B\neq 0$.

B.1 Zero magnetic field

At $B=0$ Eq. (14) simplifies to

[TABLE]

We should also keep in mind that $\rho^{\lambda}_{\sigma}(r)$ are coupled via

[TABLE]

As long as $E^{2}/v^{2}-(p_{z}-\lambda p_{0})^{2}\neq 0$, we can write a solution as

[TABLE]

where $J_{m}$ stands for the $m$-th Bessel function of the first kind. The corresponding normalisation constant can be found with the help of the integral

[TABLE]

B.2 Nonzero magnetic field

In Eq. (14) we perform a change of variable as follows $\xi=eBr^{2}/2\hbar\equiv r^{2}/2\ell_{B}^{2}$:

[TABLE]

Below we set $\hbar=1$ for the sake of simplicity. We use the two asymptotic limits of $\xi\to\infty$ and $\xi\to 0$ to build up a general solution, and we find

[TABLE]

where ${}_{1}F_{1}(a,b,z)$ is the confluent hypergeometric function of the first kind, and $\tilde{C}_{m\sigma}$ is a normalisation constant. Returning back to the original variable:

[TABLE]

It is worth noting that only one of the components of the wave function can be described by this expression. The other component should be found consistently using the original Schrödinger equation $\mathcal{H}_{\lambda}\Psi_{\lambda}=E\Psi_{\lambda}$. We choose here to express up to a normalisation constant the radial part of the lower component of the wave function:

[TABLE]

Next step is to find the upper component of the wave function using the lower one given by Eq. (19). To do so we use the initial Schrödinger equation $\mathcal{H}_{\lambda}\Psi_{\lambda}=E\Psi_{\lambda}$, where we insert

[TABLE]

and we get

[TABLE]

We use the equation above to express $\rho_{+}(r)$ in terms of found above $\rho_{-}(r)$:

[TABLE]

for $E/v\neq\lambda(p_{z}-\lambda p_{0})$. This very special case is never realized in finite-size samples. Thus, we get for $m>0$:

[TABLE]

and

[TABLE]

for $m\leqslant 0$. Using the boundary conditions in Eq. (9), we can calculate the bands in the presence of the magnetic field.

Appendix C Conductance of a cylindrical wire: transfer matrix approach

C.1 Defining scattering states

In order to compute the transfer matrix of the wire, we start by defining the scattering states of the problem. We model the leads by means of the same Hamiltonian as the wire

[TABLE]

taken at a value of the chemical potential $\mu_{\infty}$ larger than any other energy scale of the system, and we normalise the scattering states in such a way that they carry unit current in the $z$-direction. We note that regardless of whether we have or do not have a magnetic flux penetrating the wire, we set the magnetic field in the leads to zero, in order to simplify the calculation. Thus, to compute the scattering momenta we solve the following equation:

[TABLE]

which defines a set of momenta $p_{z}=\left\{p^{i}_{\infty},\,i\in\overline{1,2N}\right\}$, to which we will henceforth refer to as “scattering momenta”. We denoted the total number of scattering momenta as $2N$, because each chirality yields an equal number of scattering momenta. Indeed, Eq. (22) is symmetric under changing $p_{z}\to-p_{z}$ and simultaneously changing $\lambda\to-\lambda$. Knowing the scattering momenta, we can define a basis of scattering states:

[TABLE]

where $\rho_{\pm}^{\lambda}$ are defined in Eqs. (17) and (18). In what follows, we order and normalize the scattering states in such a way that the first and the last $N$ states are right- and left-movers, respectively, yielding all unit currents, $+1$ and $-1$, correspondingly. The current operator is defined as $\hat{j_{z}}=v\,\sigma_{z}\otimes\tau_{z}$ (hereinafter for the sake of brevity we will omit the tensor product $\otimes$). Chirality structure of the states in Eq. (23), as well as the angular parts and plane waves in the $z$-direction ensure that scattering states with different chiralities, scattering momenta and angular momenta are automatically orthogonal to each other. Thus, to ensure that the states carry unit current for a given chirality we need to satisfy:

[TABLE]

which yields:

[TABLE]

where we defined $\alpha_{i}\equiv\sqrt{\mu_{\infty}^{2}/v^{2}-(p^{i}_{\infty}-\lambda p_{0})^{2}}$, $\gamma_{i}\equiv\lambda(p^{i}_{\infty}-\lambda p_{0})-\mu_{\infty}/v$, and the function $I_{1}(m,\alpha,\beta,R)$ is defined in Appendix D.

Below we define the transverse part of the scattering states that will be consequently used to find the transfer matrix. Thus, uniting both chiralities to simplify notation, we define:

[TABLE]

C.2 Transfer matrix approach

We define a basis of scattering states

[TABLE]

and order the states in such a way that the first $N$ states are right-movers carrying unit currents from left to right, and the last $N$ states are left-movers carrying unit currents from right to left. We can write this formally in terms of the transverse parts as

[TABLE]

In the equation above $D_{ij}$ is a diagonal $2N\times 2N$ matrix defined as follows:

[TABLE]

We introduce the magnetic field parallel to the $z$ direction into the sample and we write down the Schrödinger equation for the wave function in the sample:

[TABLE]

We can represent the wave function in the basis of the scattering states $\Phi_{im}(r,\phi)$ defined above:

[TABLE]

We multiply our equation by $\Phi^{\dagger}_{im}(r,\phi)$ from the left, and we integrate it over all radii and angles:

[TABLE]

Thus we get:

[TABLE]

where we denoted

[TABLE]

The exact analytical expression for the matrix $U^{m}\equiv||U_{ij}^{m}||$ is calculated in Appendix E. The equations above show that there is no term mixing between different states with $m$ and $m^{\prime}$ while they propagate through the sample, therefore, we can compute the transfer matrix for channel $m$, and then sum up corresponding conductances over all relevant values of $m$. In order to simplify further calculations, we introduce a vector

[TABLE]

and, replacing $p_{z}\to-i\hbar\partial_{z}$, rewrite the equation above as

[TABLE]

where $\tilde{U}^{m}\equiv U^{m}\cdot D$. Assuming that $L$ is the length of the scattering region, we find the transfer matrix as matrix exponential in the following form:

[TABLE]

We must ensure that this expression is current-conserving, i.e., $T^{\dagger}_{m}\cdot D\cdot T_{m}=D$. The latter is easy to verify using the definition of the matrix exponential and mathematical induction, and we leave this proof to Appendix F.

Finally, to find the conductance we use the Landauer formula. The transfer matrix for a given channel $m$ is a $2N\times 2N$ matrix of the following form Mello et al. (1988)

[TABLE]

where $t_{m}$ and $t_{m}^{\prime}$ are the transmission matrices from left to right and right to left, respectively, and $r_{m}$ and $r_{m}^{\prime}$ are the corresponding reflection matrices. Thus, the transmission matrix $t_{m}$ is an $N\times N$ matrix, and can be defined as follows using the transfer matrix in Eq. (28):

[TABLE]

The conductance of the corresponding channel $m$ and the full conductance are, therefore, given by

[TABLE]

where $G_{0}\equiv e^{2}/h$ is the conductance quantum.

Appendix D Integrals for T-matrix calculation

In order to compute the overlap integrals we need to compute only the two following integrals:

$I_{1}(m,\alpha,\beta,R)=\int_{0}^{R}rJ_{m}\left(\alpha r\right)J_{m}\left(\beta r\right)dr$, 2. 2.

$I_{2}(m,\alpha,\beta,R)=\int_{0}^{R}r^{2}J_{m}\left(\alpha r\right)J_{m-1}\left(\beta r\right)dr$,

where $\alpha,\beta\neq 0$. We start with the first integral that has a closed analytical form. Integrating by parts we get:

[TABLE]

where $\delta_{ij}$ stands for the Kronecker’s delta symbol. The second integral does not have an analytical form for the most general case. So below we first restrict ourselves to cases where explicit analytical expressions can be derived. We start with the case where $\alpha^{2}=\beta^{2}$:

[TABLE]

Finally, in the case of $\alpha^{2}\neq\beta^{2}$ we can write down a recursive relation expressing the integral at arbitrary $m$ in terms of two integrals at $m=1$ and $m=0$ that can be calculated exactly due to symmetries and properties of Bessel functions:

[TABLE]

And finally, the recursive equation for $I_{2}(m,\alpha,\beta,R)$:

[TABLE]

Appendix E Calculating matrix elements $U_{ij}^{m}$

Below we calculate the matrix elements of the Hamiltonian in the scattering states:

[TABLE]

Since the Hamiltonian is diagonal in the chirality subspace (i.e., if we chose $\Phi^{\dagger}_{im}(r,\phi)$ and $\Phi_{jm}(r,\phi)$ to be of different chiralities, then we would get zero), we can simplify it to:

[TABLE]

The first two terms are constants, and therefore, we get:

[TABLE]

where $I_{1}\left(m,\alpha,\beta,R\right)$ is defined in Appendix D, and $\alpha_{i},\gamma_{i}$ are given in the subsection where we calculated the scattering states. We deal with the last term by rewriting it in polar coordinates:

[TABLE]

where $I_{2}\left(m,\alpha,\beta,R\right)$ is defined in Appendix D. Summing up Eqs. (34) and (35), we get the final result for $U^{m}_{ij}$.

Appendix F Current conservation proof

Below we prove that the expression for T-matrix in Eq. (28) is current-conserving. We must, therefore, prove that

[TABLE]

First, we note that

[TABLE]

Using the definition of the matrix exponential we have

[TABLE]

The inner sum for $m=0$ gives us $D$ which is cancelled by $-D$ outside of sums. It is easy to see that for $m=1$ and $m=2$ the inner sum equals [math], therefore, we assume that the inner sum is always zero. We prove this by mathematical induction. We assume that the statement is correct for $m=N$, and we prove that it implies that the statement is also correct for $m=N+1$ using the relation for binomial coefficients, namely $C_{N+1}^{k}=C_{N}^{k}+C_{N}^{k-1}$.

Appendix G Conductance as a function of the chemical potential

To demonstrate that the analytical approach to the problem yields qualitatively the same results as the numerical simulations with the Kwant package, we present in Fig. 6 the conductance curves at different values of the magnetic flux, penetrating the wire. Exactly as in the main text, here we add thermal broadening corresponding to $T=23\,$K. Since the electron charge $e$ and the reduced Planck constant $\hbar$ were set to unity in the analytical calculation, one quantum of flux $\Phi_{0}\equiv h/e=2\pi\hbar/e=2\pi$ corresponds to having $B=\Phi_{0}/\pi R^{2}=2/R^{2}=2/20^{2}=1/200$.