Field theory approach to the quantum transport in Weyl semimetals

TL;DR

This paper investigates how boundary conditions influence surface states, Fermi arcs, and transport properties in Weyl semimetals using a field theory approach, revealing tunable surface polarization and deviations from universal thermal conductance.

Contribution

It provides an analytical and numerical study of boundary condition effects on surface states, transport, and thermal properties in Weyl semimetals, including a correction to the thermal Hall conductivity.

Findings

Surface polarization and Fermi arc shape depend on boundary conditions.

Surface states' transport properties can be tuned by Zeeman terms.

Boundary effects cause deviations from universal thermal Hall conductivity.

Abstract

We analyze the structure of the surface states and Fermi arcs of Weyl semimetals as a function of the boundary conditions parameterizing the Hamiltonian self-adjoint extensions of a minimal model with two Weyl points. These boundary conditions determine both the pseudospin polarization of the system on the surface and the shape of the associated Fermi arcs. We analytically derive the expectation values of the density profile of the surface current, we evaluate the anomalous Hall conductivity as a function of temperature and chemical potential and we discuss the surface current correlation functions and their contribution to the thermal noise. Based on a lattice variant of the model, we numerically study the surface states at zero temperature and we show that their polarization and, consequently, their transport properties, can be varied by suitable Zeeman terms localized on the surface. We also provide an estimate of the bulk conductance of the system based on the Landauer-B\"uttiker approach. Finally, we analyze the surface anomalous thermal Hall conductivity and we show that the boundary properties lead to a correction of the expected universal thermal Hall conductivity, thus violating the Wiedemann-Franz law.

Equations

$$\begin{gathered} H_{\mathrm{lat}}(\mathbf{p})=\tilde{v}\left(\cos p_0-\cos p_x\right){\sigma }_x+ \\ v\left(2-\cos p_y-\cos p_z\right){\sigma }_x+v\sin p_y{\sigma }_y+v\sin p_z{\sigma }_z. \end{gathered}$$

$$H_{\mathrm{lat}}(-\mathbf{p})=U{H}_{\mathrm{lat}}(\mathbf{p})U^{-1},\text{with }U={\sigma }_x,$$

$$H(\mathbf{p})=\frac{v}{2p_0}\left(p_x^2-p_0^2\right){\sigma }_x+v{p}_y{\sigma }_y+v{p}_z{\sigma }_z,$$

$$H(-i\nabla )=-\frac{v}{2p_0}\left({\partial }_x^2+p_0^2\right){\sigma }_x-vi{\partial }_y{\sigma }_y-vi{\partial }_z{\sigma }_z.$$

$$⟨\psi H(-i\nabla )|\phi ⟩-⟨\psi |H(-i\nabla )\phi ⟩=0,$$

$$\varphi^{\dagger}(\bm{r})\sigma_{z}\psi(\bm{r})\,\bigr{|}_{z=0}=0\,,$$

$$\left(\sigma_{x}\cos\gamma+\sigma_{y}\sin\gamma\right)\psi(\bm{r})\,\bigr{|}_{z=0}=\psi(\bm{r})\,\bigr{|}_{z=0}\,.$$

$$\{{\zeta }_{\mathrm{s}}^{\pm }(\mathbf{r},p),{\zeta }_{\mathrm{b}}^{\pm }(\mathbf{r},\mathbf{p})\}$$

$$\{{\epsilon }_{\mathrm{s}}(p),\pm {\epsilon }_{\mathrm{b}}(\mathbf{p})\}$$

$$g(p_x)=\frac{p_x^2-p_0^2}{2p_0},{\tilde{p}}_z(p)=p_y\cos \gamma -g(p_x)\sin \gamma .$$

$${\epsilon }_{\mathrm{s}}(p)=v\left[g(p_x)\cos \gamma +p_y\sin \gamma \right],{\tilde{p}}_z(p)>0.$$

$${\epsilon }_{\mathrm{b}}(\mathbf{p})=v\sqrt{g(p_x{)}^2+p_y^2+p_z^2},p_z\ge 0.$$

$$p_y=-g(p_x)\cot \gamma ,{\tilde{p}}_z(p)>0.$$

$$L_{\mathrm{FA}}(\gamma )=p_0\left\{\frac{1}{\sin \gamma }+(\tan \gamma )\log \left[\frac{1+\cos \gamma }{\sin \gamma }\right]\right\}.$$

$$\Psi (t,\mathbf{r})={\Psi }_{\mathrm{s}}(t,\mathbf{r})+{\Psi }_{\mathrm{b}}(t,\mathbf{r}),$$

$$a_{\mathrm{s}}(p)$$

$$b_{\mathrm{s}}(-p)$$

$$a_{\mathrm{b}}(\mathbf{p})$$

$$b_{\mathrm{b}}(-\mathbf{p})$$

$${\Psi }_{\mathrm{s}}(t,\mathbf{r})$$

$${\Psi }_{\mathrm{b}}(t,\mathbf{r})$$

$$f(\epsilon ;\beta ,\mu )=\frac{\Theta (\epsilon )}{1+{\mathrm{e}}^{\beta (\epsilon -\mu )}},\beta =1/k_{B}T,$$

$$⟨a_{\mathrm{s}}^{†}(p)a_{\mathrm{s}}(q)⟩=f\left({\epsilon }_{\mathrm{s}}(p),{\beta }_{\mathrm{s}},{\mu }_{\mathrm{s}}\right)(2\pi {)}^2\delta (p-q),$$

$$⟨b_{\mathrm{s}}^{†}(-p)b_{\mathrm{s}}(-q)⟩=f\left(-{\epsilon }_{\mathrm{s}}(p),{\beta }_{\mathrm{s}},-{\mu }_{\mathrm{s}}\right)(2\pi {)}^2\delta (p-q).$$

$$⟨a_{\mathrm{b}}^{†}(\mathbf{p})a_{\mathrm{b}}(\mathbf{q})⟩=f\left({\epsilon }_{\mathrm{b}}(\mathbf{p}),{\beta }_{\mathrm{b}},{\mu }_{\mathrm{b}}\right)(2\pi {)}^3\delta (\mathbf{p}-\mathbf{q}),$$

$$⟨b_{\mathrm{b}}^{†}(-\mathbf{p})b_{\mathrm{b}}(-\mathbf{q})⟩=f\left(-{\epsilon }_{\mathrm{b}}(\mathbf{p}),{\beta }_{\mathrm{b}},-{\mu }_{\mathrm{b}}\right)(2\pi {)}^3\delta (\mathbf{p}-\mathbf{q}).$$

$$:{\Psi }^{†}\Psi :(t,\mathbf{r})=:{\Psi }_{\mathrm{s}}^{†}{\Psi }_{\mathrm{s}}:(t,\mathbf{r})+:{\Psi }_{\mathrm{b}}^{†}{\Psi }_{\mathrm{b}}:(t,\mathbf{r})$$

$$+:{\Psi }_{\mathrm{s}}^{†}{\Psi }_{\mathrm{b}}:(t,\mathbf{r})+:{\Psi }_{\mathrm{b}}^{†}{\Psi }_{\mathrm{s}}:(t,\mathbf{r}),$$

$$⟨:{\Psi }^{†}\Psi :(t,\mathbf{r})⟩=⟨:{\Psi }_{\mathrm{s}}^{†}{\Psi }_{\mathrm{s}}:(t,\mathbf{r})⟩+⟨:{\Psi }_{\mathrm{b}}^{†}{\Psi }_{\mathrm{b}}:(t,\mathbf{r})⟩,$$

$${\partial }_t{n}_{\mathrm{s}}(t,\mathbf{r})=-\nabla \mathbf{j}(t,\mathbf{r}),n_{\mathrm{s}}(t,\mathbf{r})=:{\Psi }_{\mathrm{s}}^{†}{\Psi }_{\mathrm{s}}:(t,\mathbf{r}).$$

$$j_x(t,\mathbf{r})$$

$$j_y(t,\mathbf{r})$$

$$j_z(t,\mathbf{r})$$

$$⟨n_{\mathrm{s}}(t,\mathbf{r})⟩$$

Full text

I Introduction

Weyl semimetals are the focus of intense theoretical (see, for example, the reviews hosur13, ; hasan17, ; yan17, ; burkov18, ; vishwanath18, ) and experimental lv2015 ; weng2015 ; xu2015 ; xu2015_5 ; xu2015_3 ; huang2015 ; zhang2016 studies. They are the prototypical example of three-dimensional systems that, despite being gapless in the bulk, display topologically protected properties wan11 ; balents11 . These properties stem from the chiral behavior of the Weyl band-touching points characterizing these materials, and they include peculiar transport phenomena such the chiral magnetic effect burkov2012 and the anomalous Hall response balents11 ; yang2011 , which are different manifestations of the quantum chiral anomaly that these systems displaynielsen1983 ; burkov2012 ; huang2015 ; zhang2016 ; goswami2013 ; parame14 ; huang17 ; lepori2018 .

One of the most striking effects of the chirality of the Weyl points is the presence of gapless chiral states that are localized on the surfaces of Weyl semimetals and are protected against disorder. These chiral surface states appear along surfaces orthogonal to the separation of the Weyl points in momentum space, and they are responsible for an anomalous Hall conductivity proportional to the distance of the projection of the Weyl points on the surface Brillouin zonebalents11 ; yang2011 .

Several experiments exploited ARPES techniques to detect these surface states xu2015_3 (see also hasan17, ; yan17, and references therein) and the corresponding Fermi arcs wan11 ; balents11 ; the shape and spin texture of the Fermi arcs, in particular, are non-universal features that strongly depend on the surface properties of the investigated materials. Despite being non-universal, however, the surface features play a fundamental role in defining some of the transport properties of Weyl semimetals: Ref. cheianov18, , for instance, shows that the boundary characteristics are crucial in evaluating the surface contribution to the current induced by the chiral magnetic effect in an alternating magnetic field beenakker16 ; Pesin16 ; cheianov17 .

In this work, we present an analytical description of the surface states, currents and anomalous Hall conductivity of a minimal model of a Weyl semimetal.

Our results are based on the study of the self-adjoint extensions of the bulk Hamiltonian to the surface witten15 ; hashimoto2017 ; seradjeh18 ; faraei18 , which determine the set of physical boundary conditions describing the interface between the Weyl semimetal and the vacuum. We consider only the ballistic regime, thus neglecting any disorder or interparticle scattering. This implies that our analysis neglects the dissipation effects from the surface to the bulk states (see, for example, Refs. gorbar16, ; resta18, for an analysis of these effects on the surface transport), and additional disorder effects on the Fermi arcs juricic17 .

Our results provide an analytical description of the surface transport beyond linear response theory. In particular, we derive the profiles of the current density and conductivity as a function of the distance from the surface and we analytically evaluate the anomalous Hall conductance as a function of chemical potential, temperature and boundary conditions of the system. We calculate the thermal noise due to the surface states and we derive the anomalous Hall thermal conductivity of the system, finding that the Wiedemann-Franz law is fulfilled by the surface states only in the zero-temperature limit.

We also verify numerically that it is possible to vary the boundary condition through suitable surface Zeeman interactions and, consequently, to change the value of the anomalous Hall conductivity of the system at non-vanishing chemical potentials.

Overall, these results provide useful tools to interpret the surface transport properties of samples of Weyl semimetals in the ballistic regime at finite temperature and non-zero chemical potential, and to suitable estimate the effects of their boundary conditions.

This article is structured in the following way: in Sec. II we introduce a low-energy description for a Weyl semimetal with two band-touching point; in particular, we emphasize the role of its boundary conditions and we construct a suitable quantum field theory for its study. In Sec. III its surface properties are studied and we present our main results about the anomalous Hall conductivity. Sec. IV provides an analytical estimate of the bulk conductance of the system. In Sec. V we compare the zero-temperature results with the numerical study of its corresponding lattice model. Sec. VI is devoted to additional properties of the system at finite temperature, focusing in particular on the surface thermal noise and anomalous thermal Hall conductivity. Finally we present our conclusions in Sec. VII.

II The model and its boundary conditions

II.1 The Hamiltonian and its spectral properties

The Weyl points in a topological semimetal appear always in pairs with opposite chirality nielsen1981 . For this reason, a minimal model describing a realistic Weyl semimetal must include two band-touching points and break time-reversal invariance, as in the case of layered intermetallic materials with a trigonal crystal structure soh2019 . Our starting point is therefore a toy model of fermions with a suitable pseudospin-1/2 degree of freedom, that can represent orbital, sublattice or spin degrees of freedom. The fermions move on a cubic lattice with a dynamics dictated by the Hamiltonian $H_{\rm lat}=\sum_{\bf p}c^{\dagger}_{\bf p}H_{\rm lat}({\bf p})c_{\bf p}$, where $c$ and $c^{\dagger}$ are two-component spinors and:

[TABLE]

In this equation, the Pauli matrices $\sigma_{i}$ act on the pseudospin, and hereafter we adopt units such that the lattice spacing is unity.

This minimal model yang2011 ; okugawa2014 ; lopez2018 displays two Weyl points in $\left(\pm p_{0},0,0\right)$ at zero energy, and we choose $\tilde{v}=v/\sin p_{0}$ in order to obtain an isotropic energy dispersion around both. For $\tilde{v}\cos p_{0}\neq 0$, the Hamiltonian (II.1) corresponds to a stack of 2-dimensional topological insulators laying on the $yz$ planes and coupled by the tunneling term along the $\hat{x}$ direction. In particular, in the two-dimensional $yz$ limit described by $\cos{p_{x}}=0$, we obtain a topological insulator with chiral gapless edge modes for $\tilde{v}\cos{p_{0}}<0$. When introducing the coupling in the $\hat{x}$ direction, these gapless modes evolve into Fermi arcs localized on the surfaces $xy$ and $xz$.

The Hamiltonian (II.1) is invariant under space-inversion symmetry,

[TABLE]

in such a way that the two Weyl points appear at the same energy. Concerning its boundaries, we neglect band-bending potentials at the surface; as a result, the Fermi arcs in this system do not display a spiraling dispersion, differently from the setups analyzed in Ref. andreev15, (see also Ref. lepori16, for the effect of trapping potentials).

To analytically study the behavior of the surface states, for small values of $p_{0}$, we approximate the low-energy behavior of (II.1) with the Hamiltonian

[TABLE]

leading to the differential operator

[TABLE]

Let us assume that the Weyl semimetal is located in the half-space ${\mathbb{R}}^{3}_{+}$ with $z>0$. Then (II.4) defines a symmetric operator on a suitable set ${\cal D}({\mathbb{R}}^{3}_{+})$ of smooth functions, which is dense in the set of square integrable functions on ${\mathbb{R}}^{3}_{+}$. The general theory (see e.g. Ref. RS, ) of such operators implies in our case that $H(-i{\bf\nabla})$ has self-adjoint extensions, which involve one real parameter. The nature of this parameter can be deduced from the condition

[TABLE]

imposed for each pair of wave functions $\psi$ and $\varphi$ in the domain ${\cal D}({\mathbb{R}}^{3}_{+})$. The previous relation is equivalent to the surface condition

[TABLE]

where † stands for Hermitian conjugation. Analogously to the 2D case of graphene mccann2004 ; akhmerov2008 , Eq. (II.6) corresponds to the physical requirement of vanishing of the probability current flowing across the surface $z=0$. It is satisfied by the (maximal) set of wave functions with an arbitrary and translationally invariant polarization of the pseudospin parallel to the surface in all its points, namely

[TABLE]

The angle $0\leq\gamma<2\pi$ parametrizes all the self-adjoint extensions $H_{\gamma}$ of the Hamiltonian (II.3) and specifies the pseudospin polarization in the plane $xy$ of all the wavefunctions on the surface $z=0$. In the physical context the parameter $\gamma$ is expected to depend in general on the electric/magnetic properties of the Weyl material and the termination of its lattice that defines the surface, as experimentally verified in Refs. souma2016, ; morali2019, .

Once $\gamma\in[0,2\pi)$ is fixed, the spectral properties of $H_{\gamma}$ are uniquely determined and concisely described in Appendix A. Summarizing, the eigenfunctions

[TABLE]

are of two different types, called in what follows surface and bulk states. The surface eigenstates $\zeta_{\rm s}^{\pm}({\bf r},p)$, given by (A.1), depend only on the two momenta $p=(p_{x},p_{y})$ parallel to the $z$-plane and decay exponentially along the $z$-axis. The bulk eigenstates $\zeta_{\rm b}^{\pm}({\bf r,p})$, given by (A.2,A.3), depend instead on all the three momenta ${\bf p}=(p_{x},p_{y},p_{z})$ and oscillate along the $z$-axis. The associated eigenvalues

[TABLE]

are expressed in terms of the combinations

[TABLE]

In particular, the surface states $\zeta_{\rm s}^{\pm}({\bf r},p)$ correspond to the domains in momentum space with positive and negative energies $\varepsilon_{\rm s}(p)$ respectively. These surface eigenvalues of the Hamiltonian are defined by:

[TABLE]

The bulk energies result instead:

[TABLE]

The eigenvectors (II.8) form a complete orthogonal basis, which together with (II.11,II.12), defines uniquelyRS the self-adjoint extension $H_{\gamma}$.

The $\gamma$-dependence of the energy spectrum is a first indication that the physics of our system depends on the boundary conditions. It is instructive to consider in this respect the zero-energy eigenstates. There are four bulk states with this property, corresponding to the spin-degenerate states at the momenta $(p_{x}=\pm p_{0},\,p_{y}=0,\,p_{z}=0)$. In addition, there is the family of surface states with vanishing energy

[TABLE]

associated with the dispersion relation (II.11). The set of all these zero-energy eigenstates form the open Fermi arc on the surface Brillouin zone. In fact, equation (II.11) defines the shape of the arc which depends explicitly on the angle $\gamma$, which gives its orientation at the limiting points $(p_{x}=\pm p_{0},\,p_{y}=0,\,p_{z}=0)$ (measured from the upper vertical direction clockwise and counterclockwise respectively, see Fig. 1). In particular, for $p_{x}=p_{0}$, the Fermi arc is always oriented orthogonally to the pseudospin polarization given by the boundary conditions (II.7). In the range $0<\gamma<\pi$, the length $L_{\rm FA}(\gamma)$ of the Fermi arcs following from (II.13) is

[TABLE]

We find Eq. (II.14) instructive because it relates the angle $\gamma$, which has an abstract mathematical origin, with the length of the Fermi arc $L_{\rm FA}$ and with $p_{0}$, which are experimentally measurable quantities in ARPES measurements xu2015 ; xu2015_3 ; souma2016 ; morali2019 .

The value $\gamma=\pi/2$ describes a straight Fermi arc connecting the projections of the Weyl points. In fact, at this point the function (II.14) reaches its minimum $L_{\rm FA}(\pi/2)=2p_{0}$. For $\gamma=0$ and $\gamma=\pi$ the orientation of the Fermi arc becomes orthogonal to the line connecting the Weyl points: in this case and, more in general, for $\pi\leq\gamma\leq 2\pi$, our choice of the function $g$ returns indeed two distinct and unbounded Fermi arcs. In tight-binding models of real materials, the two branches are always connected with each other, but our second-order approximation of the lattice Hamiltonian (II.1) describes only their behavior in proximity of the Weyl point projections and fails in depicting the global behavior of the Fermi arc across the Brillouin zone. The regime $\pi\leq\gamma\leq 2\pi$ may indeed correspond to systems and surfaces in which a Fermi arc connects Weyl points in adjacent Brillouin zones (see, for example, the recent experimental results in Ref. morali2019, ).

The above approach is very general and applies to other Hamiltonians in spaces with boundary as well. We stress in this respect that once the Hamiltonian $H$ is fixed, there is no further freedom for choosing neither the boundary condition nor the number of parameters, describing its self-adjoint extensions. In fact, the boundary condition is uniquely determined by requiring the Hermiticity (II.5) of $H$, whereas the number of parameters is fixed by the indices $n_{\pm}$ associated to $H$ (see e.g. Ref. RS, for details). If $n_{-}=n_{+}\equiv n$ according to the Von Neumann theoremRS $H$ admits a $n^{2}$-parameter family self-adjoint extensions. According to (II.7) a single angle $\gamma$ determines all self-adjoint extensions of the Hamiltonian (II.4), which implies that in our case $n=1$. This is related to our choice of a minimal 2-component model Hamiltonian for the description of the Weyl semimetal. We observe, however, that in the literature Dirac Hamiltonians involving 4-component spinors are often adopted to describe this kind of systems (see, for example, faraei18 ; faraei2019 and the analogous case for graphene akhmerov2008 ). In that case the self-adjoint extensions depend on a larger set of parameters. For the sake of simplicity we base our analysis on the model in Eq. (II.3), which allows us to capture all the main physical features of the surface transport of Weyl semimetals, without resorting to larger spinors.

Summarizing, the Hamiltonian $H_{\gamma}$ on the half space ${\mathbb{R}}^{3}_{+}$ has quite remarkable spectral properties, which represent the core of the quantum field description developed below. In particular, the main physical properties of the surface states can be derived from the Hamiltonian (II.3) for boundary conditions with $0<\gamma<\pi$. In the remainder of the paper, this regime will be assumed.

II.2 Quantum field approach

The strategy is well known and aims at the construction of a quantum field

[TABLE]

where $\Psi_{\rm s}$ and $\Psi_{\rm b}$ collect the bulk and surface contributions. The fundamental requirements are that the time evolution of $\Psi$ is generated by $H_{\gamma}$ and that $\Psi$ satisfies the canonical equal-time anti-commutation relations. In order to write the solution in explicit form we adopt for later convenience a Dirac type formulation and introduce the following notation: we label by $a_{\rm s}$ and $a_{\rm s}^{\dagger}$ the annihilation and creation operators for surface quasiparticles with $\varepsilon_{\rm s}>0$, and by $b_{\rm s}$ and $b_{\rm s}^{\dagger}$ the annihilation and creation operators for surface quasiholes, such that:

[TABLE]

where the operators $c_{\rm s}(p)$ are annihilation operators of the electrons on the surface and $\Theta$ is the Heaviside step function. Analogously, we set in the bulk

[TABLE]

where $c_{\rm b}({\bf p})$ are the annihilation operators of the electrons in the bulk. With these conventions $\Psi_{\rm s}$ and $\Psi_{\rm b}$ take the form

[TABLE]

in terms of the complete system (II.8) of eigenfunctions of $H_{\gamma}$ and the creation and annihilation operators (II.16-II.19), satisfying the canonical anti-commutation relations. We observe that the surface and the bulk components (II.20,II.21) obey separately the equation of motion. This is not he case for the equal-time canonical anti-commutation relations, which follow from the completeness of the energy eigenstates and hold therefore only for the total field (II.15).

At this point, the choice of representations of the oscillator algebras, generated by (II.20,II.21), is the only freedom we are left with. Since our goal is to study the Weyl semimetals at finite temperature and density, we will adopt below the Gibbs representation bratteli , keeping in general different (inverse) temperatures and chemical potentials $\{\beta_{\rm s},\mu_{\rm s}\}$ and $\{\beta_{\rm b},\mu_{\rm b}\}$. In terms of the Fermi distribution

[TABLE]

the non-vanishing two-point functions are:

[TABLE]

For the bulk excitations one has instead

[TABLE]

We observe that there is no interference between surface and bulk oscillators and that all higher point correlation functions can be expressed in terms of (II.23-II.26).

At this stage we can construct and investigate the basic physical observables of the system. Let us consider for instance the particle density operator

[TABLE]

where $:\cdots:$ stands for the normal ordering with respect to the creation and annihilation operators (II.16-II.19). Since the surface and the bulk operators have vanishing mixed two-point functions, one finds

[TABLE]

which implies that at the level of mean values there is no interplay between surface and bulk degrees of freedom. This is a general feature, which allows us to treat separately the mean values of the surface and bulk currents as well.

Let us stress finally that the above quantum field theory setting works directly in the thermodynamic limit of our system.

III Surface currents and Hall conductivity

III.1 Mean value of the surface current

The anomalous Hall conductivity of Weyl semimetals is determined by their Fermi arcs. In the minimal model (II.1), a simple decomposition of the 3D Hamiltonian (II.3) into a two-dimensional set of systems parameterized by $p_{x}$ shows that the number of surface states defining the the Fermi arc is proportional to the distance $2p_{0}$ between the projections of the Weyl points in the surface Brillouin zone. In particular, each 2D system defined by $H_{p_{x}}(p_{y},p_{z})$ constitutes a Chern insulator with chiral gapless edge modes for $-p_{0}<p_{x}<p_{0}$, and this implies that the contribution to the anomalous Hall conductivity of the surface states in the ballistic regime at zero temperature and half filling is given by $\sigma_{H}=e^{2}p_{0}/\pi h$ balents11 ; yang2011 . This value of the anomalous Hall conductivity is universal and does not depend on the boundary conditions of the system. Other physical quantities, as, for example, the behavior of $\sigma_{H}$ at finite chemical potential, depend instead on the boundary condition (II.7). In the following, we analytically study the physics of the surface modes focusing on several characteristics which are determined by the value black $0<\gamma<\pi$ of the surface polarization.

To investigate the transport properties of the system generated by the surface states we consider the surface current ${\bf j}(t,{\bf r})$, satisfying the continuity equation

[TABLE]

Using the equations of motion, one finds for ${\bf j}(t,{\bf r})$:

[TABLE]

Now, adopting the two-point functions (II.23,II.24) and the explicit form (A.1) of the surface eigenfunctions one obtains

[TABLE]

with $\tilde{p}_{z}(p)$ given by (II.10). As expected, the mean values (III.5-III.7) are time independent (invariance under time translations) and $(x,y)$-independent (invariance under space translations in the $(x,y)$-plane). Moreover, they manifestly satisfy the continuity equation (III.1) and the relation

[TABLE]

Equations (III.5-III.7) provide the distribution of the density of the electrons and the current generated by the surface states as a function of the distance $z$ from the boundary. We point out that the particle density in Eq. (III.5) is a different quantity compared to the density of the surface states, which, instead, can be easily derived from Eq. (II.11) (see Section V). The surface current is dissipationless in our model, since we neglect scattering effects from surface to bulk states, and it is responsible for the anomalous quantum Hall conductivity. Its dependence from $x$ and $y$ is trivial due to the translational invariance in these directions. The surface current along the $\hat{x}$ direction has a vanishing expectation value due to the contributions of states with positive and negative $p_{x}$ canceling each other.

III.2 The anomalous Hall conductivity

In the ballistic regime, the surface currents are dissipationless; this implies that, in a typical 2-terminal transport measurement with two external leads attached to a Weyl semimetal scatterer, the surface states can be considered in equilibrium with the leads they originate from. This is analogous to the standard quantum Hall devices in two dimensions and it allows us to consider the distribution of the surface states at a fixed chemical potential. The Hall conductivity can thus be obtained from the derivative of the expectation value of the current with respect to the chemical potential inherited by the source lead.

In order to investigate the anomalous Hall conductivity it is instructive to consider the $p_{x}$ components $\hat{j}_{y}(p_{x},z,\mu_{s})$ of the surface current density such that $j_{y}(z,\mu_{s})=\int dp_{x}\hat{j}_{y}(p_{x},z,\mu_{s})$. In particular, we can define the local contribution of each surface state labelled by $p_{x}$ to the total differential conductivity $\sigma_{H}$:

[TABLE]

such that:

[TABLE]

The momentum integrals in Eqs. (III.7) and (III.9) cannot be expressed in simple closed form in the general case. However, we can invert the order of the integrations, and determine $\sigma_{H}(\mu_{s})$ by integrating first the space coordinate $z$, and then the momenta. This allows for a general expression of the Hall conductivity as a function of the chemical potential $\mu_{s}$, the inverse temperature $\beta$ and the boundary polarization $\gamma$:

[TABLE]

where we reintroduced the Planck constant $\hbar$ for clarity. Here ${\rm Li}_{1/2}$ is a polylogarithm function mathbooks .

In the zero-temperature limit, $\beta_{s}\to\infty$, the anomalous Hall conductance (III.11) becomes:

[TABLE]

This expression is proportional to the density of surface states; it depends in general on the specific choice of the function $g(p_{x})$ and, in our case, it is valid for $vp_{0}+2\mu_{s}\cos(\gamma)>0$. For values of $\mu_{s}$ outside this regime at $T_{s}=0$, the differential anomalous Hall conductivity vanishes because the surface states are either completely empty (for $\cos\gamma>0$) or completely filled (for $\cos\gamma<0$); the dispersion relation (II.11) has indeed a minimum or a maximum for $\cos\gamma\gtrless 0$ respectively due to the constraint $\tilde{p}_{z}(p)>0$. Eq. (III.12) provides the known result $\sigma_{H}=e^{2}p_{0}/\pi h$ for $\mu_{s}\to 0$ and arbitrary boundary conditions $0<\gamma<\pi$, such that this value is universal (in the ballistic regime). Furthermore, we observe that for $\gamma=\pi/2$, the conductivity $\sigma_{H}$ is independent of $\mu_{s}$; for $0<\gamma<\pi/2$, it increases with $\mu_{s}$, whereas it decreases for $\pi/2<\gamma<\pi$ [see Fig. 2 (a)]. This shows the importance of the boundary conditions in the determination of the non-universal corrections to the anomalous Hall conductivity.

For finite temperature the system acquires a non-zero surface conductance also for values of the chemical potential such that $vp_{0}+2\mu_{s}\cos(\gamma)\leq 0$ [see Fig. 2 (b)]. The temperature dependence of the Hall conductivity (III.11) is non-trivial. For $\mu_{s}=0$, the conductivity is a function of the rescaled temperature $k_{B}T_{s}|\cos\gamma|/v$ and the parameter $p_{0}$. In particular, $\sigma_{H}$ is non monotonic in $T_{s}$ (see Fig. 3): by increasing $T_{s}$ from zero, the conductivity decreases from the universal value $e^{2}p_{0}/\pi h$ to a minimum which depends on $p_{0}$. Then it increases again and, asymptotically, it grows proportionally to $\sqrt{p_{0}k_{B}T_{s}|\cos\gamma|/v}$. This behavior reflects the particular density of states of the surface modes which is determined by the dispersion (II.11) combined with the constraint $\tilde{p}_{z}(p)>0$.

III.3 The case of straight Fermi arc $\gamma=\pi/2$

For the particular case of straight Fermi arcs, $\gamma=\pi/2$, the anomalous Hall conductivity acquires its universal value $\sigma_{H}=e^{2}p_{0}/\pi h$ and it is independent of both $T_{s}$ and $\mu_{s}$, as it can be verified by the limit of Eq. (III.11) [see, for instance, Fig. 2 (b)].

Besides the total Hall conductivity, the definition of the current density (III.7) allows us also to define the current as a function of the distance $z$ from the surface, thus the local contribution (III.9) to the Hall conductivity. To analyze the local transport properties of the system we begin by investigating the case of a straight Fermi arc, $\gamma=\pi/2$, for which $\varepsilon_{\rm s}$ depends on $p_{y}$ only. For this boundary conditions $\left\langle j_{y}(z)\right\rangle$ becomes independent of the chemical potential and temperature:

[TABLE]

hereafter the function $D_{+}$ labels the Dawson integral:

[TABLE]

The current (III.13) at $\gamma=\pi/2$ is consistent with the following local conductivity of the surface state labelled by $p_{x}$:

[TABLE]

independent of $\mu$. We plot $\hat{\sigma}(p_{x},z)$ for $\gamma=\pi/2$ in Fig. 4; its decay length in the bulk is given by $g(p_{x})^{-1}$: as expected, the states of the Fermi arc with $p_{x}$ approaching the Weyl point projections at $p_{x}=\pm p_{0}$ progressively penetrate deeper in the bulk and their contribution to the conductivity is weaker for $z=0$ but decays slower with the distance. For $p_{x}$ in the center of the Fermi arc, instead, the surface states are more localized and their contribution is stronger for $z=0$.

III.4 The zero-temperature limit

For generic values of the boundary condition parameter, the integrals in Eqs. (III.5) and (III.7) cannot be written in simple closed forms. The results simplify in the zero-temperature limit, that allows us to clearly evaluate the role of the boundary conditions and the chemical potential. For $T\to 0$ we can substitute the Fermi function (II.22) with:

[TABLE]

Concerning the operators $n_{\rm s}$ and $j_{y}$ in Eqs. (III.5) and (III.7), in the limit $\mu_{s}\to 0$, the difference of the $T_{s}\to 0$ limits of the Fermi distributions in the integral is meant to return the contribution of the zero-energy Fermi arc states. From Eq. (III.8) we observe that the angle $\gamma$ fixes the direction of the surface current density in our model, which gets inverted for $\pi<\gamma<2\pi$. This inversion can be obtained in the physical systems only for Fermi arcs that connect the two projections of the Weyl points by winding across the Brillouin zone: our approximation of the lattice Hamiltonian (II.1), however, is accurate only for momenta relatively close to the Weyl points (and small values of $p_{0}$); therefore, Eq. (II.3) does not describe correctly the physics of the whole surface Brillouin zone, resulting in unbounded Fermi arcs for $\pi<\gamma<2\pi$.

The momentum integral in Eq. (III.5) determines the decay of the total density of the surface states in the bulk; in particular, for large values of $z$ we may approximate the density with the asymptotic expansion:

[TABLE]

This relation is defined only for $p_{0}+2(\mu_{s}/v)\cos\gamma>0$. For values of $\mu_{s}$ outside this range of validity, the surface states are either completely filled (for $\cos\gamma<0$) or completely empty (for $\cos\gamma>0$) and $\left\langle n_{\rm s}\right\rangle\approx-p_{0}/\left(4\pi^{2}z^{2}\cos\gamma\right)$. The relation (III.17) implies that the expectation value of the density and current density of the surface state decays as $z^{-2}$ in the bulk. As emphasized by several works gorbar16 ; resta18 , this may yield in turn relevant scattering processes between surface and bulk states. The limit $z\to 0$, instead, provides an estimate of the surface density exactly at the boundary and it is given by:

[TABLE]

After defining $\hat{\sigma}(p_{x},z,\mu_{s})$ from Eq. (III.9), the integral in the $p_{x}$ momenta in Eq. (III.10) returns:

[TABLE]

where

[TABLE]

The conductance density (III.19) is represented in Fig. 5 as a function of the boundary condition parameter $\gamma$ and the distance from the surface $z$. It decays asymptotically as:

[TABLE]

All the previous results are based on the specific model (II.3) which is characterized by a single pair of Weyl cones and can be considered a good approximation for experimental systems with two Weyl cones only soh2019 . In the case of materials with more Weyl points, however, we expect that each Fermi arc will provide a contribution to the anomalous Hall conductance given by the previous equations (III.19-III.15). The additivity of these contributions can be verified in the limit of small $p_{0}$, in which several pairs of Weyl points are sufficiently distant from each other in momentum space.

IV Bulk conductance

The total conductance of a Weyl semimetal is given by the sum of its surface and bulk contribution. In the previous section, we evaluated the contribution of the conductivity due to the surface modes. The bulk conductance, instead, can be estimated through the Landauer approach by considering a Weyl semimetal connected to two external leads and evaluating the transmission amplitude associated to each of their modes. This approach has been applied in Ref. baireuther14, for the case of a pair of overlapping Weyl cones at the phase transition between topological and normal insulators and further detail can be found in Ref. beenakker2006, for the two-dimensional case of graphene.

Here we apply the same approach to approximate the conductance of the three-dimensional Weyl semimetal described by the Hamiltonian (II.3). To this purpose, we consider a system divided into three regions in the $\hat{y}$ directions, with width $W$ in the $\hat{x}$ and $\hat{z}$ directions. In order to estimate the bulk conductance, we consider periodic boundary conditions in the $\hat{x}$ and $\hat{z}$ directions, such that the momenta $p_{x}$ and $p_{z}$ are quantized in units of $2\pi/W$. The regions at $y<0$ and $y>L$ constitute two infinite leads, and we describe them with the Hamiltonian (II.3) and a chemical potential $\mu_{\rm lead}\to\infty$, in such a way that the two leads are effectively in a metallic phase with a large density of states. The central region $0\leq y\leq L$, instead, is characterized by $\mu=0$ and models a Weyl semimetal with chemical potential lying at the level of the two Weyl points.

We impose the continuity of the wavefunction at the interfaces $y=0$ and $y=L$ such that, for each value of $p_{x}$ and $p_{z}$ we obtain the transmission probability (see Appendix B for more detail):

[TABLE]

which generalizes in a straightforward way the result in Ref. baireuther14, for an arbitrary $g(p_{x})$. For our specific choice of $g(p_{x})$ and periodic boundary conditions along the $\hat{x}$ and $\hat{z}$ directions, the bulk conductance at vanishing chemical potential results:

[TABLE]

where we adopted the notation $p_{x,z}=2\pi n_{x,z}/W$. We observe that, in general, $G_{\rm b}$ is not scale invariant due to the parameter $p_{0}$, and it does not depend only on $L/W$. To evaluate the general behavior of the conductance $G_{\rm b}$ it is useful to distinguish two regimes: a fine-tuned regime, in which $p_{0}$ is an integer multiple of $2\pi/W$, and a standard regime in which $p_{0}$ is not a multiple of $2\pi/W$.

The fine-tuned regime is special because the Weyl points lie exactly on one of the momenta of the Brillouin zone, therefore there are two bulk zero-energy modes that contribute with a quantum of conductance to the bulk transport, independently on $L$, and the conductance decreases asymptotically to $2e^{2}/h$ for $L/W\to\infty$. This case is analogous to the result of Ref. baireuther14, for periodic boundary conditions, and the conductance of the Dirac semimetal in baireuther14, is recovered in the limit of large $p_{0}L$, where the contribution of the bulk modes to $G_{\rm b}$ is compatible with having two well-separated Weyl cones.

The most realistic scenario is the one with $p_{0}\neq n_{0}2\pi/W$ ($n_{0}\in\mathbb{N}$). Differently from the fine-tuned regime, the system has a vanishing conductance in the limit $L/W\to\infty$ for $p_{0}L\gg 1$. In order to estimate the conductance in this case, we may consider the behavior of the transmission probability in proximity of the two Weyl points, where $\mathcal{T}$ is maximized. Let us consider the case $n_{x}\approx 2\pi p_{0}W$. The maximum value of $\mathcal{T}$ can be approximated by observing that:

[TABLE]

such that $\mathcal{T}_{max}\approx\cosh^{-2}\left(2\pi{L}/{W}\right)$. By considering the decay of $\mathcal{T}$ with $n_{x}$ and $n_{z}$ away from the Weyl points, we conclude that the system has a vanishing conductance for $L/W\to\infty$.

Let us finally address the limit of close Weyl points, that describes a system approaching a phase transition in which the topological semimetal phase may be gapped. We model this regime by considering $p_{0}W\ll 1$ and we consider a wire geometry, thus $L\gg W$. In this situation we can estimate the behavior of the conductance in the following limits: for $p_{0}L\ll 2$, the conductance goes to $G_{\rm b}\to e^{2}/h$ because of the contribution of the term $n_{x}=n_{z}=0$; for $p_{0}L\to\infty$, instead, the conductance vanishes.

For a system with surfaces at $z=1$ and $z=W$, we assume that the surface at $z=W$ displays an opposite polarization with respect to the one in $z=1$ (as suggested by our numerical results based on the model (II.1), see also Ref. beenakker18, ). In this case:

[TABLE]

By considering the Hamiltonian (II.3), this implies that the values of $p_{z}$ must be taken as:

[TABLE]

analogously with the two-dimensional case of graphene beenakker2006 . This quantization of the momenta must be considered in calculating the bulk conductance $G_{\rm b}$ and the Weyl semimetal acquires a total conductance of the form:

[TABLE]

We conclude that the total conductance indirectly depends on the boundary conditions via the quantization of the momenta orthogonal to the surfaces: The bulk term in $G$ is indeed non-universal and depends, in general, on the choice of $g(p_{x})$ and the boundary conditions, whereas the surface term represents the universal anomalous Hall conductance for vanishing chemical potential and temperature. We finally observe that, for more general boundary conditions with independent polarizations $\gamma$ and $\gamma^{\prime}$ on the surfaces at $z=1$ and $z=W$, the determination of the correct set of momenta $p_{z}$ gives rise, in general, to non-analytical solutions.

Eq. (IV.1) describes the transmission probability in the system (II.3) with two Weyl points; we observe, however, that it can be extended also to materials characterized by well-separated dipoles of Weyl points, a common experimental situation hasan17 ; yan17 . Indeed, the transmission coefficient $\mathcal{T}(p_{x},p_{z})$ decays exponentially with the distance in momentum space from the pair of Weyl points, hence the value of $G_{\rm b}$ is dominated by the states in proximity of the Weyl pair. Therefore we expect that in a material with small $p_{0}$ and pairs of Weyl points sufficiently far from each other, the contribution of each pair to the bulk conductance will approximately add to each other. In this case, the resulting conductance $G_{\rm b}$ can be approximated with the sum over all the Weyl dipoles of the value (IV.2) of a single Weyl pair.

V Numerical comparisons

In this section we verify numerically the analytical results of the previous sections. For this task, we use the Kwant code kwant to simulate the following Hamiltonian of spin-1/2 fermions on the cubic lattice, corresponding to Eq. (II.1):

[TABLE]

with $b=v\left(2+\cot p_{0}\right)$. In the following, we will focus on finite size systems where the Weyl semimetal constitutes a ”scattering region” with size $W_{x}\times L\times W_{z}$ (in units of the lattice spacing $a\equiv 1$), and we will adopt different boundary conditions.

The results of the previous sections rely on the value of the surface polarization $\gamma$ defined in Eq. (II.7). Therefore, as a first step, we measure the value $\gamma(p_{0})$ for several values of $p_{0}\in[0,\frac{\pi}{2}]$. To this purpose we consider open boundary conditions along the $\hat{z}$ direction and periodic boundaries along $\hat{x}$ and $\hat{y}$ directions, in order to have only surfaces orthogonal to $\hat{z}$. In particular, we analyze systems with dimensions $L$ and $W\equiv W_{x}=W_{z}$ up to 150 and we estimate the parameter $\gamma$ by evaluating the surface polarization:

[TABLE]

here the expectation value is taken over the single-particle eigenstate of (V.1) corresponding to the lowest-energy eigenstate with positive energy. We verified that this state corresponds to a linear superposition of states localized on the two surfaces at $z\gtrsim 1$ and $z\lesssim W_{z}$ and it belongs to the (hybridized) Fermi arcs for the system with this geometry. In particular, we considered the site at ${\bf r_{0}}=\left(\frac{W}{2},\frac{L}{2},1\right)$ on the surface at $z=1$. We verified that the polarization $\gamma$ does not depend on the $x$ and $y$ coordinates when considering periodic boundary conditions in these directions and we checked that its dependence on $y$ is very weak also for open boundary conditions.

The relation (V.2) is easily obtained from the explicit form of the states localized on the surface at $z=0$ in the continuum model of the previous sections (see Eqs. (II.7) and (A.1)). In the entire range $p_{0}\in[0,\frac{\pi}{2}]$, we notably find only small deviations around the value $\gamma=\pi/2$, of the order of $10^{-2}$. This suggests that, in the thermodynamic limit, the lattice model (V.1) is indeed defined by the boundary polarization $\gamma=\pi/2$. Therefore, to be able of varying the parameter $\gamma$ and study its effect on the surface states, we introduce the following additional surface terms to the Hamiltonian:

[TABLE]

These surface interactions correspond to opposite Zeeman-like terms for the pseudospin of the system aligned along the $\hat{y}$ direction and localized on the two surfaces at $z=1$ and $z=W_{z}$. We adopted only Zeeman fields in the $\hat{y}$ directions because we verified that, for analogous values of the coupling constants, the Zeeman terms along $\hat{x}$ cause a much weaker effect on the surface eigenstates.

The engineering of such a surface term in physical systems strongly depends on the material properties. When the pseudospin is related to the occupation of specific orbitals within a unit cell, then, depending on the lattice termination, suitable weak voltage gates on the surface may favor or disfavor these occupations, thus giving rise to such effective Zeeman-like terms for the pseudospin on the surface. Effective surface Hamiltonians of the kind (V.3) can be obtained by the detailed analysis of the electrostatic properties of the interfaces between Weyl semimetals and other materials (or gates) which has been addressed in several works (see, for instance, Ref. goerbig17, ). Moreover, in the case of Weyl semimetals breaking time-reversal symmetry, it is also possible that the pseudospin is associated to the physical spin of the electrons; in such a situation, ferromagnets in contact with the surface may constitute useful tools for the engineering of such a surface Hamiltonian.

Our approach is dictated also by a more mathematical reason: we observe that any lattice Hamiltonian with a finite number of degrees of freedom, such as $H_{\rm lat}$, is self-adjoint without the necessity of specifying any extension through boundary conditions. Thus, we introduce the surface Hamiltonian $H_{\rm s}$ to compensate for the lack of this freedom in parameterizing the boundary conditions and we verify in the following that, indeed, the surface terms (V.3) modify the physical system in such a way that the analytical low-energy description presented in the previous sections provides an accurate approximation for the behavior of the system. In particular, we verify that, also for $B\neq 0$, the surface polarization is linked to the shape of the Fermi arcs and the density of surface states as a function of the energy, thus to $\sigma_{H}$.

The properties of the surface states for $B=0.2v$ are illustrated in Fig. 6. These results correspond to periodic boundary conditions along $\hat{x}$ and $\hat{y}$ and are obtained by diagonalizing the Hamiltonian in the subspaces defined by the conserved momenta $p_{x}$ and $p_{y}$. To select the surface states, we considered a threshold $A_{1/4}$ corresponding to the squared amplitude of the wavefunctions in the interval $z\in[0,W_{z}/4]$. The results in Fig. 6 correspond to all the states fulfilling $A_{1/4}>0.4$, thus sufficiently localized close to $z=1$. This threshold is however arbitrary and, in the comparison with the analytical low-energy model, the states we selected from the numerical simulation are only a subset of the corresponding surface states fulfilling the constraint $\tilde{p}_{z}>0$. Hence we expect some deviation of the numerical results from the analytical predictions for values of $\tilde{p}_{z}$ close to zero, thus close to the Weyl point projections and, more in general, close to the edges of the region where the analytical model predicts the existence of surface states, which are depicted as black lines on the surface Brillouin zone in all the panels of Fig. 6. Fig. 6(a) illustrate the probability $A_{1/4}$ for the eigenstates of the lattice model: in general, the amplitude $A_{1/4}$ decreases by approaching the predicted boundaries (black lines) and the Weyl point projections (green dots); small irregularities can be observed for $p_{y}=0$ due to the hybridization of the surface states in opposite surfaces. For increasing values of $|p_{y}|$, thus of $|\varepsilon_{s}|$, the domain of the selected numerical surface states is smaller than the analytical prediction; this is mostly due to the differences for energies comparable with $vp_{0}$ between the low-energy analytical model and the lattice model. The discrepancies in the domains is also partially due to the constraint $A_{1/4}>0.4$ which implies an underestimation of the surface states domain in the numerical data.

Fig. 6(b) displays the polarization (V.2) of the surface eigenstates. The value of $\gamma$ weakly varies as a function of the conserved momenta, and for $B=0.2v$, its average is $\gamma\approx 1.906$. We adopted this value in the analytical determination of the surface state domain and of the shape of the Fermi arc. The Fermi arc derived by the analytical prediction (II.13) matches very well the numerical results and, in general, for energies close to zero, the agreement between analytical and numerical models is very good, as shown by the comparison of the panels (c) and (d).

Fig. 7 displays the values of the surface parameter $\gamma$, as a function of the boundary field $B$. We compare two different ways of obtaining the estimate of $\gamma$: (i) we consider the average of the expectation value of the polarization (V.2) for the surface states, selected based on the constraint $A_{1/4}>0.4$ for the energy interval $\varepsilon_{\rm s}\in\left[-0.5v,0.5v\right]$; (ii) we consider the value of $\gamma$ obtained from a fit of the density of surface states.

Concerning the approach (ii), the number of surface states for an energy interval $d\varepsilon_{\rm s}$ is derived from Eqs. (II.11) and (II.10) and it results in:

[TABLE]

As mentioned in Section III, this quantity differs from the particle number on the surfaces, defined via Eq. (III.5).

We use Eq. (V.4) to perform a one-parameter fit of the numerical data concerning the number of states with $A_{1/4}>0.4$ of a system with dimension $300\times 300\times 90$ for the energy range $\varepsilon_{\rm s}\in\left[-0.3v,0.3v\right]$ with intervals $d\varepsilon_{\rm s}=0.02v$. Two examples of the fit result and density of surface states are shown in Fig. 8 for $B=0.3v,-0.05v$. This method is potentially affected by a larger systematic error because of the lattice and finite size effects, that determine considerable oscillations of the density for the surface states and cause a deviation from the analytical model for $\epsilon_{s}$ approaching $vp_{0}$.

The numerical results above are obtained for an isolated system with periodic boundary conditions in the $\hat{x}$ and $\hat{y}$ direction. In the following we focus on the transport properties of the Weyls semimetal in contact with two external leads.

In particular, to probe the surface conductivity, we include two semi-infinite semimetallic leads in the system, and we connect them to the Weyl scattering region on the facets at $y=1$ and $y=L_{y}$, thus reproducing the geometry of the previous section. These leads are described by the same Hamiltonian V.1 and they are approximately characterized by the same chemical potential of the scattering region.

This configuration is required to probe directly the conductivity of the scatterer, avoiding non-universal effects from the leads and the interfaces between leads and scattering region: such aspect is especially important for the dynamics of surface states. In more detail, to probe the low-energy physics around the Weyl nodes, we set the chemical potential of the leads $\mu_{\mathrm{lead}}=0$, and in the scattering region $\mu_{\mathrm{bulk}}$ at a slightly larger value (typically $\mu_{\mathrm{bulk}}/v\approx 10^{-3}-10^{-2}$). To simplify our numerical calculations, we maintained periodic boundary conditions in the $\hat{x}$ direction, thus diagonalizing the system in different Hilbert space sectors labeled by $p_{x}$, and we considered system sizes with $W_{x}=W_{z}\equiv W$.

The resulting conductance of the system corresponds, in the thermodynamical limit, to the anomalous Hall conductance $G_{H}=W\,\sigma_{H}$. This is a direct consequence of our choice of the leads, with the same Hamiltonian and vanishing chemical potential as the scatterer: the bulk density of states vanishes also in the leads, and the only contribution to the total conductance is given by the surface states. Under these conditions, even for small values of $L_{y}$ and $W$, we observe a clear quantization of $G_{H}$ in units of $e^{2}/h$ (see the inset of Fig. 9), analogously with the behavior of Weyl semimetal nanowires bardarson19 . Such quantization corresponds to the number of states in the Fermi arc at energy equal to $\mu\approx 0$; for the data of Fig. 9, $G_{H}$ assumes odd values of the conductance quantum. For the systems considered in the present Section, the total number of states in the Fermi arc for $\mu\approx 0$ and $0<\gamma<\pi$ results: N^{(tot)}_{s}=\Big{[}2p_{0}/(\frac{2\pi}{W_{x}})\Big{]}=\Big{[}p_{0}W_{x}/\pi\Big{]}, \big{[}\big{]} denoting the integer part. This value leads to the universal anomalous Hall conductance in Eq. (IV.5).

For sufficiently large $W$, thus with negligible hybridization of the surface states on the opposite surfaces, we find that the Hall conductivity $\sigma_{H}={G_{H}}/{W}$ tends to $1/6={p_{0}}/{\pi}$, the universal value predicted in Eq. (III.12) at $\mu=0$ (see Fig. 9). We observe, however, that $\sigma_{H}$ presents some discontinuities as a function of the width of the system, due to finite size effects. For the largest system size we probed, $W=140$, we obtained $\sigma_{H}(W=140)\approx 0.164e^{2}/h$.

We conclude by analyzing the transport properties of the bulk. To this purpose we consider systems with periodic boundary conditions in the $\hat{x}$ and $\hat{z}$ directions. By maintaining vanishing chemical potentials in both the leads and the scatterer, we calculated the bulk conductance $G_{\rm b}$ for $p_{0}=\pi/6$, in both the fine-tuned and standard regimes (the fine-tuned regime is given by $W=12n$ with $n\in\mathbb{N}$, such that $p_{0}W$ is a multiple of $2\pi$). For example, we consider $W=24$ and $W=25$. In the fine-tuned regime, already for $L_{y}=10$ and $W=24$, we measure $G_{\rm b}=2{e^{2}}/{h}$, up to an error of $10^{-13}$. This verifies the existence of the two expected zero-energy non-evanescent bulk states, which correspond to the Weyl band-touching points and match the continuum model description. However, differently from the previous section and the calculation in Eqs. (IV.1) and (IV.2), our choice of vanishing chemical potential in the leads implies a vanishing of their density of states, such that there is no other contribution to the transport in this regime apart from these two zero-energy bulk states.

A comparison with Eqs. (IV.1) and (IV.2) can be performed, instead, by adopting metallic leads, which may be simply modeled by a cubic lattice Hamiltonian of fermions with spin-independent nearest-neighbor hopping terms; in this case a high density of states occurs in the leads, similarly to the Landauer-Buttiker approach of Eq. (IV.2) (where the limit $\mu_{\mathrm{lead}}\to\infty$ is adopted).

The analytical estimate of the bulk conductance in Eq. (IV.2) is compared with the numerical result with metallic leads in Fig. 10 for $W=30$, $\mu=10^{-3}$, and $p_{0}=\pi/6$, as a function of $L$. The yellow dots represent the results from Eq. (IV.2), while the blue dots denote the numerical data. The analytical results underestimate the numerical conductance typically by a factor $\sim 3$, thus it provides only an estimate of the order of magnitude of the bulk conductance. This discrepancy is due to the different kinds of leads considered, Weyl semimetallic leads for the analytical estimate and metallic leads for the numerical results, and to additional interface effects that can stem from different boundary conditions at the interface between the wavefunctions of the leads and the scatterer than the ones considered in Appendix B.

VI The effects of temperature

VI.1 Thermal noise

The field theoretical approach described in Sec. III can be used to obtain an estimate of the current and conductivity behavior also as a function of temperature. We saw that, for $\gamma=\pi/2$, the surface current $\left\langle j_{y}\right\rangle$ does not depend on $T$ and $\mu$, whereas the average value of the currents along $x$ and $z$ always vanishes, as dictated by the symmetries of the system. In the ballistic regime, the chiral surface modes can be considered in thermal equilibrium with the leads they originate from; therefore, for the geometry considered in the previous section, we can introduce surface chemical potential $\mu_{\rm s}$ and inverse temperature $\beta_{\rm s}$ equal to the parameters of the lead at $y<0$ for $\sin\gamma>0$ (again $0<\gamma<\pi$ ). In the case $\sin\gamma<0$, which is in general not well-defined in our continuum model, the surface current changes direction and $\mu_{\rm s}$ and $\beta_{\rm s}$ would instead be derived from the lead at $y>L$.

From the definition of the current operator, we can derive the spectral density of the surface current noise at frequency $\nu$:

[TABLE]

For $\gamma=\pi/2$, the calculation can be explicitly done. In the limit $\nu\to 0$, for $x_{1}=x_{2}$ and $z\equiv z_{1}=z_{2}$, we obtain:

[TABLE]

which does not depend on $y_{1}$ and $y_{2}$ and decays asymptotically as $z^{-4}$. The shot noise vanishes,

[TABLE]

consistently with the chiral nature of the surface states. The pure thermal limit gives instead:

[TABLE]

which respects the Johnson-Nyquist law.

Eq. (VI.2) is the two-point correlation function of the surface current density only. A complete calculation of the noise spectrum, however, must take into account also the bulk-bulk and bulk-surface correlations at different points.

VI.2 Bulk currents

The expectation value (III.7) accounts exclusively for the surface current; the total current density is given by

[TABLE]

and includes both surface and bulk contributions. Employing the two-point functions (II.23-II.26) and the explicit form (A.1,A.2) of the surface and bulk eigenfunctions one obtains the following mean values:

[TABLE]

[TABLE]

The surface state contribution gives the first term in the curly brackets. The original contribution of the bulk states (A.2) involves an integration over $d^{3}p$. Using the Cauchy integral formula and the pole structure of the scattering matrices (A.4,A.5), one can perform the $p_{z}$ integral, which leads to the second term. In this operation one uses the fact that restricting $p_{z}$ in $\varepsilon_{\rm b}({\bf p})$ to the poles of (A.4,A.5) one gets precisely $\varepsilon_{\rm s}({p})$. For this reason the surface energy $\varepsilon_{\rm s}$ enters both the bulk and surface contributions.

In general, the assumption of having different temperatures and chemical potentials for the surface and bulk states is justified in the ballistic regime: our model neglects scattering terms and interactions between bulk and surface states. Such terms naturally appear in a physical system due to disorder gorbar16 or electron-phonon coupling resta18 and determine a relaxation time $\tau$ beyond which bulk and surface states equilibrate. However, if $\tau\ll L/(v\sin\gamma)$, we may assume that a quasiparticle does not equilibrate during the transport between the two leads in the geometry discussed in the previous section. In this situation, as a first approximation, it is legitimate to choose different bulk and surface Fermi distributions.

If instead we consider the system with $\beta_{\rm s}=\beta_{\rm b}$ and $\mu_{\rm s}=\mu_{\rm b}$, from (VI.9) one immediately infers that the mean value of total current $\left\langle{\bm{J}}(t,{\bm{r}})\right\rangle$ vanishes everywhere. In this case we are dealing in fact with the thermodynamic limit of a system at equilibrium. This result may not hold in finite systems: by including two parallel surfaces in our model with a finite separation, thus two different scattering matrices, the expectation value of the local current density at equilibrium may in general be different from zero and depend on position, consistently with previous numerical results displaying persistent currents in small Weyl semimetals lopez2018 ; zhang2018 .

Also concerning the spectral density of the noise, at equilibrium, the correlations between bulk operators provide a contribution identical to the surface states, thus doubling the result (VI.4). Additionally, one should consider also correlations between surface and bulk currents whose computation go beyond the scope of this work.

VI.3 Thermal Hall conductivity

Heat currents are typically more difficult to measure than electric currents. In experimental systems the heat transport is determined not only by the electrons, but also by the phonons propagating in the material. When we restrict our attention to the heat transport along the Weyl semimetal surfaces, the electronic contribution to the heat current typically scales with $T^{2}$, due to the chiral dispersion of the Fermi arcs, whereas the phonons are free of propagating in any direction, leading to a typical Stefan-Boltzmann behavior proportional to $T^{4}$ (or $T^{5}$ for more refined models wellstood1994 ). Therefore, analogously to the heat transport in quantum Hall setups, the electronic contribution dominates for low temperatures; we conclude that the calculation of the electronic heat conductivity of our model provides an approximate description of the heat transport for low temperatures.

Similarly to the charge transport umansky2017 , Weyl semimetals display an anomalous Hall effect also for the heat transportgoswami2013 ; Lundgren2014 . Former analysis based on the bulk properties of Weyl semimetals suggests the existence of a universal value of the thermal Hall conductivity, $\kappa_{H}=p_{0}{\pi k_{B}^{2}T_{\rm s}}/{3h}$ at vanishing chemical potential, such a value fulfills the Wiedemann-Franz law and is consistent with the energy transport of the chiral states on the surface kane1997 ; cappelli2002 . In the following we analyze the energy transport on the surface of our model and we show that, in general, the boundary conditions yield non-universal corrections of the thermal transport, determining, in turn, a violation of the Wiedemann-Franz law.

Based on the total field $\Psi$ and the surface field $\Psi_{\rm s}$, the total and surface energy currents of our system, in the $\hat{y}$ direction orthogonal to the Weyl point separation, are defined by:

[TABLE]

The electric and energy currents (III.3,VI.6,VI.10,VI.11) generateCallen the bulk and surface heat currents,

[TABLE]

Analogously to the case of the electric transport, one can extract the thermal Hall conductance $\kappa_{H}$ from the surface heat current (VI.13) as a function of the surface temperature $T_{\rm s}$. We obtain:

[TABLE]

In the case $\gamma=\pi/2$ of a straight Fermi arc, the Fermi arc dispersion depends only on $p_{y}$ and the previous expression simplifies. In particular, in the limit of vanishing chemical potential, we recover the predicted universal result goswami2013 :

[TABLE]

This result is consistent with the Wiedeman-Franz law, expected for non-interacting fermions in the absence of boundaries. For $\gamma=\pi/2$, indeed the density of surface states is constant, thus causing no corrections to the universal value (VI.15) as a function of temperature or chemical potential. We verify the Wiedeman-Franz law by comparing the heat and electric currents. The two current densities are proportional to each other and, in particular, we get:

[TABLE]

corresponding to a Lorentz number:

[TABLE]

for both bulk and surface states.

It is instructive to analyze the heat current (VI.17): its dependence from surface and bulk temperature is given by $T_{\rm s}^{2}-T_{\rm b}^{2}$. In a more realistic description we may consider a temperature that varies smoothly from the surface to the bulk, such that $T_{\rm s}^{2}-T_{\rm b}^{2}\to-T\partial_{z}{T}/2$. Such substitution determines indeed a Hall heat current orthogonal to the gradient of the temperature.

In the general case $\gamma\neq\pi/2$, the Lorentz ratio $\mathcal{L}$ of the Hall transport has a non-trivial dependence from the boundary polarization $\gamma$ (see Fig. 13) and the behavior of the anomalous thermal Hall conductivity can be estimated numerically from the integral in Eq. (VI.14). At $\mu_{s}=0$ we numerically observe that the Lorentz ratio is a non-trivial function of the parameter $k_{B}T_{s}|\cos\gamma|/v$ (see Fig. 11). In the limit $T_{s}\cos\gamma\to 0$, the universal ratio (VI.18) is recovered, but, in the case $\gamma\neq\pi/2$, the Wiedemann-Franz law is in general violated for $T_{s}>0$ due to the boundary conditions and the behavior of the density of surface states.

For $\mu_{s}\neq 0$, $\kappa_{H}$ depends separately on $T_{s}$ and $\gamma$ and the Lorentz ratio displays a rich behavior (see Figures 13 and 14). From Eqs. (III.11) and (VI.14) it is possible to derive that the Lorentz ratio depends solely on the parameter:

[TABLE]

For $\cos\gamma>0$, $\alpha=0$ when the chemical potential coincides with the minimum of the energy band of the surface states; for $\cos\gamma<0$ and $\alpha=0$, $\mu$ lies on the surface energy maximum instead. The values $\alpha>0$ thus correspond to the chemical potential lying within the energy band of the surface states; whereas for $\alpha<0$, the chemical potential lies outside.

From the integral in Eq. (VI.14) we derive:

[TABLE]

where ${\rm Li}$ labels polylogarithm functions mathbooks .

In Fig. 12 we illustrate the general behavior of $\mathcal{L}$ as a function of $\alpha$. The sign of the square bracket in (VI.19) is particularly important because it determines the low-temperature behavior of the Lorentz ratio. In the limit $\alpha\to+\infty$, the Lorentz ratio converges to the standard value (VI.18) and the Wiedemann-Franz low is fulfilled for $T_{s}\to 0$; for $\alpha>0$, indeed, the curve in Fig. 12 reproduces the same results as Fig. (11). On the contrary, for negative values of $\alpha$, the Lorentz ratio increases and diverges as $k_{B}^{2}\alpha^{2}/e^{2}$ for $\alpha\to-\infty$. For negative values of (VI.19), the surface states acquire indeed an insulating behavior in the low-temperature limit, since the chemical potential lies outside their energy band. The value at $\alpha=0$, instead, corresponds to the high-temperature limit for both the regimes. From Eq. (VI.20) we derive $\mathcal{L}(0)\approx 1.634$, which defines the high-temperature limit for all the values of $\gamma\neq\pi/2$ and all the values of $\mu_{s}$, compatibly with the regime of validity of the low-energy Hamiltonian (II.3). In this high-temperature limit, the Wiedemann-Franz law for the anomalous Hall conductivity is always violated for $\gamma\neq\pi/2$, despite the non-interacting nature of our model.

The dependence of the Lorentz ratio from the parameter $\alpha$ is reflected in the behaviors depicted in Figures 13 and 14 when considering the thermal transport as a function of the boundary conditions and the chemical potential respectively.

In Fig. 13 we consider the behavior of $\mathcal{L}$ for $\mu_{s}>0$ as a function of $\gamma$ for several values of the temperature: we immediately observe that there is a point, for $\gamma>\pi/2$, in which all the curves cross and the Lorentz ratio is independent of the temperature. This point coincides with $\alpha=0$ in Eq. (VI.19), thus for the chemical potential lying on the extremum of the surface energy band. For all the values of gamma on the left of this crossing point, we observe a non-monotonic behavior of $\mathcal{L}$ with the temperature, similar to the case of $\mu_{s}=0$ (Fig. 11). These values of $\gamma$ span the range $\alpha>0$ in Fig. 12 and, in this regime, the Wiedemann-Franz law is recovered for $T\to 0$. The values of $\gamma$ larger than the crossing point, instead, span the regime $\alpha<0$ due to $\mu_{s}>0$. For decreasing temperature and for each value of $\gamma$ in this range, the Lorentz ratio diverges.

A corresponding behavior is obtained also as a function of $\mu_{s}$. Fig. 14 illustrates the Lorentz ratio as a function of $\mu_{s}/v$ for different values of the temperature at $\gamma=\pi/3$. There is a crossing point in which $\mathcal{L}$ does not depend on the temperature. This is the value of $\mu_{s}$ such that $\alpha=0$ in Eq. (VI.19). On the left of this point the Lorentz ratio increases by lowering the temperature. On the right of the crossing point, instead, $\alpha>0$ and we recover the Wiedemann-Franz limit for $T_{s}\to 0$.

VII Conclusions

In this work we analyzed in detail the surface transport properties of a toy model of Weyl semimetal with two band-touching points, thus breaking time-reversal symmetry, which can be adopted, for example, to describe layered intermetallic materials with magnetically induced Weyl points soh2019 . Additionally, we expect that our results can be easily extended to account for multiple pairs of Weyl points, for materials in which such pairs are well-isolated in the Brillouin zone.

Our analysis stems from the self-adjoint extensions of the bulk Weyl Hamiltonian and allows us to focus on the non-universal properties of the transport that depend on the boundary conditions. In particular, the set of boundary conditions we considered is defined by a single angle $\gamma$ that determines both the pseudospin polarization at the surface of the system and the shape of the Fermi arcs.

From experimental data souma2016 ; morali2019 and ab initio simulations sun2015 of Weyl semimetals, it is known that the lattice termination of the system (related to the pseudospin polarization) strongly affects the properties of the Fermi arcs. In full generality, the boundary spin polarization may depend on the conserved momenta of the system, thus it can rotate along the Fermi arcs in the surface Brillouin zone sun2015 . Our model considers instead a simpler case with a constant pseudospin polarization $\gamma$ independent of the momenta, as obtained by the self-adjoint Hamiltonian extensions with local boundary conditions. The polarization $\gamma$ we consider is not necessarily equivalent to the physical spin of the system, thus our model is compatible with former results. Despite its simplicity, our work demonstrates that it is possible to consistently include the boundary conditions in a field theoretical determination of the surface transport properties.

We studied in detail the Fermi arcs of the system and their contribution to the anomalous Hall conductivity by calculating the corresponding current density, which typically decays with the square of the distance from the surface.

We derived a general formula for the anomalous Hall conductivity as a function of temperature, chemical potential and boundary parameter for both the electric and thermal transport. This allowed us to verify that, in general, the Wiedemann-Franz law for the anomalous Hall transport is fulfilled only in the very limit $T\to 0$, whereas the Lorentz ratio presents a non-trivial behavior as a function of chemical potential and boundary conditions for finite temperatures.

We additionally estimated the bulk conductance of our two-Weyl-point model based on a Landauer-Büttiker approach, and we verified with numerical calculations our predictions at zero temperature based on a suitable lattice model. In particular, we show that, also in lattice systems, it is possible to vary the boundary conditions through the introduction of suitable surface Zeeman interactions. Also in this case, the boundary polarization, the shape of the Fermi arcs and the density of surface states are linked to the same parameter $\gamma$.

Acknowledgements

We thank Ajit C. Balram, Ion Cosma Fulga, Pietro Novelli, Marco Polini, and Fabio Taddei for useful discussions and the developers of Kwant. M. B. was supported by a research grant (Project nr. 25310) from Villum Fonden.

Appendix A The eigenvectors of the self-adjoint extensions $H_{\gamma}$

We describe here the explicit form of the $H_{\gamma}$-eigenvectors and summarize their basic properties. The surface eigenstates are given by

[TABLE]

where $\varepsilon_{\rm s}(p)$ and $\tilde{p}_{z}(p)$ are defined by (II.11) and (II.10). The exponential decay along the $z$-axis is worth mentioning.

For the bulk eigenstates one has

[TABLE]

[TABLE]

where $\varepsilon_{\rm b}({\bf p})$ is defined by (II.12) and

[TABLE]

[TABLE]

The factors (A.4,A.5) satisfy

[TABLE]

and have a simple physical interpretation: they represent the scattering matrices for particles and antiparticles with incoming momenta ${\bf p}_{\rm r}$, which are reflected from the boundary $z=0$ and have final momenta $\bf p$.

The fundamental property of the surface and bulk states is that they form a complete system

[TABLE]

which is the main ingredient for constructing the canonical quantum field $\Psi$ defined by (II.15,II.20,II.21). For proving (A.7) one can proceed as follows. One starts by considering the three-dimensional integral in the square brackets, which gives the right hand side plus a rest. The latter can be reduced to a two-dimensional integral by integrating over $p_{z}$, using the Cauchy integral formula, $z>0$ and the analytic properties of the scattering matrices $S_{\pm}(\bf p)$ in the complex $p_{z}$-plane. At this point the rest precisely cancels the second term in the square brackets. This computation is very instructive because it shows that the surface states are generated by the bound states (poles in the upper half complex $p_{z}$-plane) of the scattering matrices (A.4,A.5).

Appendix B Estimation of the bulk conductance

For the evaluation of Eq. (IV.1) we consider wavefunctions at $E=0$ of the kind:

[TABLE]

Here, $S$ and $D$ label the source and drain external leads; $r$ and $t$ are the reflection and transmission amplitudes, such that $|t|^{2}+|r|^{2}=1$, and they depend on $p_{z}$ and $g(p_{x})$. The spinors in the leads are polarized along the $\hat{y}$ axis, consistently with the approximation $\mu_{\rm lead}\to\infty$ and $|p_{y}|\gg|p_{x}|,|p_{z}|$. The chemical potential in the central region is $\mu=0$ and

[TABLE]

due to the vanishing energy. To estimate the transmission probability $\mathcal{T}=|t|^{2}$ we impose the specific boundary conditions:

[TABLE]

we obtain:

[TABLE]

from which we derive Eq. (IV.1). A more rigorous estimate of the transmission coefficient should take into account more general boundary conditions at $y=0$ and $y=L$, which, also in this case, could be classified through the self-adjoint extensions of the full Hamiltonian of the system.