Quantum oscillations in Weyl semimetals - a surface theory approach

TL;DR

This paper develops a surface theory for Weyl semimetals to analyze quantum oscillations, revealing how bulk tunneling and Weyl point proximity influence phase offsets in oscillation patterns.

Contribution

It introduces an effective surface theory that captures surface-bulk interactions and explains phase offsets in quantum oscillations near Weyl points.

Findings

Surface quantum oscillations are affected by bulk tunneling.

Proximity of Weyl points causes phase offsets in oscillations.

Bulk chiral Landau levels are essential for orbit completion.

Abstract

We develop an effective surface theory for the surface states of a Weyl semimetal. This theory includes the peculiar Fermi arc states on the surface as well as leakage of the states from the surface to the bulk. Subjecting the model to a magnetic field perpendicular to the surface results in quantum oscillations. The oscillations are different from the usual ones since they do not involve a closed Fermi surface cross section. It has been shown previously that the Quantum oscillations can be understood semiclassically as resulting from motion of electrons on the surface Fermi arcs as well as tunneling through chiral Landau levels associated with the bulk. In this work we develop an effective surface theory and use it to analyze the quantum oscillation in the semiclassical regime and beyond. Specifically, we show that when a pair of Weyl points are close to each other the surface quantum oscillations acquire a phase offset which originates from the bulk. While the surface states are responsible for a large part of the electron motion, tunneling through the bulk is necessary for completing the orbit. This tunneling makes use of the bulk, zero energy, chiral Landau level in each Weyl node. When the nodes are close in momentum space their chiral levels overlap and a gap at zero energy is formed. This gap causes the phase offset in the surface quantum oscillations.

Equations

$$\begin{array}{r}{H}_W=\hslash \mathbf{v}\cdot \mathbf{k}{\sigma }_0+\sum _{i,j=x,y,z}\hslash h_{ij}{k}_i{\sigma }_j,\end{array}$$

$$\begin{array}{r}{ϵ}_n=\frac{\pi \hslash v}{k_a{\ell }_B^2+L}\left(n+\gamma \right),\end{array}$$

$$\begin{array}{rl}{H}_{3D}=& t_s\left(\sin k_x{\sigma }_x+\sin k_y{\sigma }_y+\sin k_z{\sigma }_z\right)\\ & +\left(m+t^{\prime }(2-\cos k_x-\cos k_y)\right){\sigma }_z,\end{array}$$

$$\begin{array}{rl}{E}_{\pm }=& \pm [t_s^2\left({\sin }^2{k}_x+{\sin }^2{k}_y\right)\\ & +{{\left(t_s\sin k_z+m+t^{\prime }(2-\cos k_x-\cos k_y)\right)}^2]}^{\frac{1}{2}}.\end{array}$$

$$\begin{array}{r}{H}_{3D}=\left(\begin{array}{ccccccc}{H}_0& R& 0& 0& \dots & 0& 0\\ R^{†}& H_0& R& 0& \dots & 0& 0\\ 0& R^{†}& H_0& R& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \dots & R^{†}& H_0\end{array}\right),\end{array}$$

$$\begin{array}{r}R=\left(\begin{array}{cc}-\frac{1}{2}& -\frac{t_s}{2}\\ \frac{t_s}{2}& \frac{1}{2}\end{array}\right).\end{array}$$

$$\begin{array}{rl}{H}_0& =t_s\sin k_x{\sigma }_x+\left(m+2-\cos k_x+t_s\sin k_z\right){\sigma }_z\\ & \equiv g_1(\vec{k}){\sigma }_x+g_3(\vec{k}){\sigma }_z.\end{array}$$

$$G_{\text{eff}}(i{\omega }_n)={\left[G_s^{-1}(i{\omega }_n)-T^{†}{G}_b(i{\omega }_n)T\right]}^{-1},$$

$$H_s$$

$$T=$$

$$\begin{array}{r}{G}_{\text{eff}}={\left[G_0^{-1}(i{\omega }_n)-R^{†}{G}_{\text{eff}}R\right]}^{-1},\end{array}$$

$$\begin{array}{rl}& R=\left(\begin{array}{cc}0& 0\\ t_s& 0\end{array}\right),\\ & H_0=g_1(\mathbf{k}){\sigma }_z-g_3(\mathbf{k}){\sigma }_x.\end{array}$$

$$\begin{array}{r}G(i{\omega }_n,k)=\left(\begin{array}{cc}{G}_{\text{eff}}^{(2)}& G_{\text{eff}}^{(3)}\\ G_{\text{eff}}^{(3)}& G_{\text{eff}}^{(1)}\end{array}\right),\end{array}$$

$$\begin{array}{rl}{G}_{\text{eff}}^{(1)}& =\frac{1}{2t_s^2(i\omega +g_1)}\left(t_s^2+(i\omega {)}^2-g_1^2-g_3^2\pm \sqrt{p}\right),\\ G_{\text{eff}}^{(2)}& =\frac{i\omega +g_1}{2t_s^2{g}_3^2}\left(-t_s^2+(i\omega {)}^2-g_1^2-g_3^2\pm \sqrt{p}\right),\\ G_{\text{eff}}^{(3)}& =\frac{1}{2t_s^2{g}_3}\left(t_s^2-(i\omega {)}^2+g_1^2+g_3^2\mp \sqrt{p}\right).\end{array}$$

$$\begin{array}{rl}p=& -4t_s^2(i\omega -g_1)(i\omega +g_1)+(-t_s^2-(i\omega {)}^2+g_1^2+g_3^2{)}^2.\end{array}$$

$$\begin{array}{r}{G}_{\text{eff,ret}}^{(1)}=\frac{t_s^2-g_3^2+|t_s^2-g_3^2|}{2t_s^2\left(\omega +\sin k_x+i\eta \right)},\end{array}$$

$$\begin{array}{r}{A}_{\text{eff}}^{(1)}\propto \{\begin{array}{ll}(1-\frac{g_3^2}{t_s^2})\delta (\omega +\sin k_x)& \text{for ts2g32<1}\\ 0& \text{otherwise}\end{array}\end{array}$$

$$\begin{array}{r}|m+1+\sin k_z|<1\Rightarrow m+\sin k_z<0.\end{array}$$

$$\begin{array}{r}H=v_F{k}_x{\sigma }_x+v_F{k}_y{\sigma }_y+\chi v_z{k}_z{\sigma }_z,\end{array}$$

$$\begin{array}{r}{v}_z{\pi }_z=-i\sqrt{\frac{v_F{v}_z}{2{\ell }_B^2}}(a^{†}-a),v_F{\pi }_x=\sqrt{\frac{v_F{v}_z}{2{\ell }_B^2}}(a^{†}+a).\end{array}$$

$$\begin{array}{r}{H}^2=\left(\frac{2v_F{v}_z}{{\ell }_B^2}(a^{†}a+\frac{1}{2})+v_F^2{k}_y^2\right){\sigma }_0-\chi \frac{v_F{v}_z}{{\ell }_B^2}{\sigma }_y,\end{array}$$

$$\begin{array}{r}{E}_n=\pm \sqrt{\frac{2v_F{v}_z}{{\ell }_B^2}n+v_F^2{k}_y^2}.\end{array}$$

$$\begin{array}{r}\Delta ϵ=\sqrt{2}C(m^{\ast }{v}_F^2{\hslash }^2{\omega }_c^2{)}^{\frac{1}{3}}\exp \left(-\frac{2}{3}{\left(\frac{k_{a}a}{2}\right)}^3\frac{m^{\ast }{v}_F^2}{\hslash {\omega }_c}\right),\end{array}$$

$$\begin{array}{r}C=\frac{1}{\sqrt{2}}{\left(2\sqrt{\frac{\pi }{2}}\frac{\Gamma (\frac{7}{4})}{\Gamma (\frac{1}{4})}\right)}^{\frac{2}{3}}\approx 0.523\end{array}$$

$$H_2=t_2\sin (k_z),$$

$$\begin{array}{rl}{H}_0=& \sum _{k_x,k_z}[\sum _{n=1}^q(m+2+\sin \left(kz-2\pi /q\cdot n\right)){\sigma }_z{c}_n^{†}{c}_n\\ & +\sum _{n=1}^{q-1}(\frac{i{\sigma }_x-{\sigma }_z}{2})c_n^{†}{c}_{n+1}\\ & -\frac{1}{2}{e}^{-i{k}_x}{\sigma }_z{c}_1^{†}{c}_q+\frac{1}{2i}{e}^{-i{k}_x}{\sigma }_x{c}_1^{†}{c}_q+\text{h.c.}],\end{array}$$

$$\begin{array}{r}R=\left(-\frac{1}{2}{\sigma }_z+\frac{1}{2i}{\sigma }_y\right)\otimes 1_q,\end{array}$$

$$\begin{array}{r}H=v_F{k}_x{\sigma }_x+v_F{k}_y{\sigma }_y+(\gamma +\frac{k_z^2}{2m^{\ast }}){\sigma }_z,\end{array}$$

$$\begin{array}{r}H=B{v}_F\tilde{z}{\sigma }_x+v_F{k}_y{\sigma }_y+(\gamma +\frac{k_z^2}{2m^{\ast }}){\sigma }_z,\end{array}$$

$$\begin{array}{r}H=(m^{\ast }{v}_F^2{\omega }_c^2{)}^{\frac{1}{3}}\left(Z{\sigma }_x+(\Gamma +\frac{K^2}{2}){\sigma }_z\right),\end{array}$$

Full text

I Introduction

Topological semimetalic states of matter have been theoretically predicted in 2007 by MurakamiMurakami (2007) and a simple lattice model of a Weyl semimetal (WSM) was proposed by Burkov and Balents in 2011Burkov, Hook, and Balents (2011). The early studies led to further theoretical workWan et al. (2011); Burkov and Balents (2011); Witczak-Krempa and Kim (2012); Chen and Hermele (2012); Heikkilä, Kopnin, and Volovik (2011) as well as the experimental realization of a type-I WSM in TaAsWeng et al. (2015); Huang et al. (2015); Lv et al. (2015a, b); Yang et al. (2015); Belopolski et al. (2016), TaNbXu et al. (2015a), NbPLiu, Ye, and Qi (2013) and TaPXu et al. (2015b). Type-II WSM behavior has been predicted in $\text{WTe}_{2}$Soluyanov et al. (2015), $\text{TaIrTe}_{4}$Koepernik et al. (2016) as well as $\text{MoTe}_{2}$Sun et al. (2015); Wang et al. (2016). Additionally, several theoretical predictions of Weyl semimetals have been made including in $\text{SrSi}_{2}$Huang et al. (2016), $\text{HgCr}_{2}\text{Se}_{4}$Xu et al. (2011) and $\text{Mo}_{x}\text{W}_{1-x}\text{Te}_{2}$Chang et al. (2016).

Weyl semimetals are three dimensional materials which are characterized by band crossing points in momentum space. Close to these points the energy disperses linearly with momentum and the system is described by the Weyl Hamiltonian. In three dimensions the Weyl nodes are robust against deformations and appear in pairs due to time reversal or inversion symmetry. The Weyl points are Berry curvature monopoles in momentum space and are responsible for the chiral anomalyNielsen and Ninomiya (1983); Adler (1969); Bell and Jackiw (1969) which is exhibited in a variety of transport effectsHosur and Qi (2013) such as negative magnetoresistanceNielsen and Ninomiya (1983), the anomalous quantum Hall effectBurkov and Balents (2011); Xu et al. (2011); Yang, Lu, and Ran (2011); Zyuzin and Burkov (2012); Vafek and Vishwanath (2014); Goswami and Tewari (2013); Chen, Wu, and Burkov (2013), the chiral magnetic effectJian-Hui et al. (2013); Vazifeh and Franz (2013); Chen, Wu, and Burkov (2013) as well as coupling between magnons and plasmonsLiu, Ye, and Qi (2013).

When the Fermi level of the system is close to the energy of the Weyl nodes the bulk low energy properties of the system are governed by these points. Expanding the Hamiltonian to linear order in momentum about the Weyl points leads to,

[TABLE]

where the indices $i$ and $j$ run over spatial directions, $\sigma_{0}$ is a two dimensional unit matrix while other $\sigma_{n}$s are Pauli matrices. The second term in the Hamiltonian is the usual Weyl Hamiltonian in which the matrix $h_{ij}$ determines the spinor direction. The determinant of $h_{ij}$ is $\chi=\pm 1$, the chirality of the Weyl point. The first term is unique to condensed matter systems, as it breaks Lorentz invariance. Unlike the second term, it is proportional to the unit matrix. This term, named ’tilt’ does not influence the spin direction but does have an effect on the spectrum as it tilts the Weyl cone in the energy-momentum space. In type-I WSM the tilt term is not strong enough to alter the nature of the Weyl point while in type-II WSM the tilt causes the system to have a finite density of states at the Weyl node energy. The density of states at zero energy alters the low temperature thermodynamic properties such that type I and type II WSMs can be distinguished experimentallySoluyanov et al. (2015).

Pairs of Weyl nodes with opposite chirality and at different momentum also lead to Fermi arcs. These are zero energy states localized on the sample surface and characterized by surface momentum along a line connecting the Weyl nodes projected to the surface momentum. Their existence can be understood in the context of the quantum anomalous Hall effect in two dimensions by slicing the system into many two dimensional systems. In each slice a topological invariant can be calculated and the presence of the Weyl nodes indicate a jump in this number. The Fermi arc states are therefore the collection of zero energy chiral edge states.Burkov, Hook, and Balents (2011) In this paper we concern ourselves with studying the Fermi arc states and their response to magnetic field.

In order to address the surface states we define the surface Brillouin zon (SBZ) which is made of the allowed momenta parallel to the surface. The behaviour of single particle wavefunctions perpendicular to the surface is described in real space. As mentioned above, the zero-energy surface states form Fermi arcsWan et al. (2011) which have been observed in photoemission measurement in TaAsXu et al. (2015c). Recently, Potter et al.Potter, Kimchi, and Vishwanath (2014); Zhang et al. (2016) proposed that the surface states on the arcs respond to a perpendicular magnetic field and produce magnetic quantum oscillations. Semi-classically the electrons on the Fermi arcs on the two surfaces connect though the bulk chiral bulk Landau level and can therefore move on a closed magnetic path. The semiclassical analysis predicts the following level quantization:Potter, Kimchi, and Vishwanath (2014)

[TABLE]

where, $k_{a}$ is the length of the Fermi arc, $l_{B}$ is the magnetic length, $L$ is the thickness of the slab and $v$ is the Fermi velocity. The phase shift $\gamma$ not accessible semiclassically. The period of the oscillations has two contributions. The first is given by the propagation of the electrons along the arc and is therefore proportional to the arc length $k_{a}$. The second contribution is proportional to $L$ since it is due to tunneling through the bulk. This dependence of the quantum oscillations on the geometry of the sample provides a signature of the mixed surface-bulk magnetic path, unique to WSM. Experimental evidence of this dependence was found in $\text{Cd}_{3}\text{As}_{2}$Moll et al. (2015).

This rest of this paper is structured as follows. In section II we introduce the model and develop the effective surface theory. In section III we apply a magnetic field to the system and use the effective surface theory to study the surface quantum oscillations and study the short arc length regime. In section IV we compare our findings to a full numerical treatment of the WSM slab.

II Effective Surface Theory

II.1 The model

We begin with a two-orbital tight binding modelChang and Yang (2015),

[TABLE]

where time-reversal symmetry $H(\mathbf{k})=\sigma_{y}H^{*}(-\mathbf{k})\sigma_{y}$ is broken by the second term. The choice to break time reversal symmetry is not limiting as one could devise a similar, time reversal invariant Hamiltonian, while breaking inversion symmetry. Our choice here is made in order to work with small matrices. The bulk energies are given by,

[TABLE]

Throughout the remainder of the paper, we will set $t^{\prime}=1$ and measure energy in units of $t^{\prime}$. The nodes of Eq. (4) show that the model has different phases depending on the parameters $m$ and $t_{s}$. For example, when keeping $t_{s}$ fixed to 1 and varying $m$, one finds the following phases:

• •

For $m>t_{s}$ the model is gapped and trivial.

• •

At $m=t_{s}$ a gap closure appears at $\mathbf{k}=(0,0,-\pi/2)$.

• •

For $-t_{s}<m<t_{s}$ the gap closure splits into two Weyl nodes which recombine for $m=-t_{s}$ at $\mathbf{k}=(0,0,\pi/2)$.

• •

For $-3t_{s}<m<-t_{s}$ there are two pairs of Weyl nodes which appear for $m=-t_{s}$ at $\mathbf{k}=(\pi,0,-\pi/2)$ and $\mathbf{k}=(0,\pi,-\pi/2)$.

• •

The Weyl points recombine again at $m=-3t_{s}$ at $\mathbf{k}=(\pi,0,\pi/2)$ and $\mathbf{k}=(0,\pi,\pi/2)$.

• •

For $-5t_{s}<m<-3t_{s}$ two Weyl nodes emerge at $\mathbf{k}=(\pi,\pi,-\pi/2)$ when $m=-3t_{s}$ and split when $m$ is decreased. The two points recombine at $\mathbf{k}=(\pi,\pi,\pi/2)$ when $m=-5t_{s}$.

• •

For $m<-5t_{s}$ the model is again gapped and trivial.

We choose to work in one of the Weyl semimetal regimes above, where $|m|<t_{s}$. The analysis can be easily extended to other regimes.

We work in a slab geometry with $(010)$ surfaces such that $k_{x}$ and $k_{z}$ remain good quantum numbers. We work partially Fourier transformed operators, $c_{i_{y},\vec{k}}$, where $\vec{k}=(k_{x},k_{z})$ and $i_{y}$ is a discrete coordinate in the $\hat{y}$ direction ranging between $1$ and the number of layers, $N_{y}$. With this definition the Hamiltonian can be written as a block matrix in $i_{y}$ and spin space. The diagonal blocks $H_{0}(\vec{k})$ represent hopping within the $x-z$ layer while the off diagonal ones represent inter-layer terms. The inter-layer terms connect neighbouring layers $i_{y}$ and $i_{y}\pm 1$ and are given by $2\times 2$ matrices in spin space, $R$. We can therefore write the system’s three dimensional Hamiltonian matrix as:

[TABLE]

with

[TABLE]

and

[TABLE]

In this geometry, surface states appear and are arranged in arcs in the surface Brillouin zone. The arcs extend between the projections of the Weyl points onto the surface Brillouin zone. In our lattice model and in the regime we choose to work, the arc connecting the two Weyl nodes is a straight line along the $k_{z}$-axis of length $k_{a}=2\arccos{(\frac{m}{t_{s}})}$.

II.2 Effective surface propagator

In order to derive an exact effective surface theory, we treat the surface degrees of freedom independently from the bulk and integrate out the bulk degrees of freedom. We denote the sites with $i_{y}=1,N_{y}$ as ’surface’ and the sites with $1<i_{y}<N_{y}$ as bulk. The matrix $H_{3D}$ is rearranged in this manner and one can identify 4 blocks: $H_{b}$ is an $2(N-2)\times 2(N-2)$ matrix which contain the bulk terms, $H_{s}$ is a $4\times 4$ matrix containing terms within the two surfaces and the off diagonal blocks $T$ and $T^{\dagger}$ couple the two. They are $2(N-2)\times 4$ and $4\times 2(N-2)$ matrices. Following Marchand and FranzMarchand and Franz (2012), the expression for the surface Green’s function,

[TABLE]

where $G_{b,s}(i\omega_{n})=-(i\omega_{n}-H_{b,s})^{-1}$ are the uncoupled bulk and surface Greens functions, respectively. It is important to note that the propagator in Eq. 8 is not directly related to an effective Hamiltonian since it contains a finite lifetime due to the decay of surface states into the bulk.

II.3 Numeric evaluation of the effective Green’s function in a slab geometry

The Hamiltonian of the full three dimensional system was given in Eq. 5. In a slab geometry the bulk Hamiltonian is the same as the full Hamiltonian with the first and last two rows and columns removed. The surface Hamiltonian reads:

[TABLE]

and the coupling matrix is given by:

[TABLE]

with $N-2$ such rows.

For realistic system sizes, the matrices involved are very large and the solution to Eq. (8) is only accessible numerically. We perform the numerical calculation and present the spectral function $A=-\frac{1}{\pi}\Im{[\text{Tr}(G_{\text{eff}})]}$ in Fig. 1, Panel 1 (a) shows the spectral function for the combined top and bottom surfaces for fixed $k_{z}$. It is plotted as a function of $k_{x}$ and the energy. One can see the top and bottom surface states in the gap. The positive slop line represents surface states on the top surface while the negative slope line contains states confined to the bottom surface. These states exist for any $k_{z}$ between the two Weyl points and at zero energy they form the Fermi arcs. The Fermi arc states are also seen in the spectral function cut at zero energy, as seen in Figs 1 (b-d) for different arc lengths, controlled by the parameter $m$.

II.4 Analytic Green’s function at low energy, semi-infinite sample

In order to advance analytically, we change the geometry of the system to a semi-infinite slab in $y$-direction by taking $N_{y}\rightarrow\infty$. Thus, we set $H_{s}=H_{0}$ and the only non-zero block of the matrix $T$ is equal to $R$. Looking at a semi-infinite slab has the advantage that the system with one layer removed is identical to the system before removing the layer. We therefore envision that we’re looking for the effective Green’s function for the $n$th layer when the effective Green’s function for the $(n+1)$th is known. Since the system is unchanged by removing a single layer, the two Green’s functions above are identical. This leads to the following recursive equation:

[TABLE]

where we have defined the uncoupled Green’s function, $G_{0}=(i\omega_{n}-H_{0})^{-1}$. The equation for $G_{\text{eff}}$ is now simply a $2\times 2$ matrix equation and can be solved analytically. Nonetheless, for general parameter values, the solution is quite complicated and not very insightful. It is therefore useful to simplify it by transforming the system via the unitary transformation $U=\exp{\left(-i\frac{\pi}{4}\sigma_{y}\right)}$, which corresponds to a rotation around the $y$-axis in orbital space. In addition, we set $t_{s}=t^{\prime}=1$ in order to simplify the result. This leads to,

[TABLE]

The solution for the Green’s function reads,

[TABLE]

where

[TABLE]

and

[TABLE]

Analytic continuation then yields the retarded/advanced Green’s functions. The low energy part of the spectrum is governed by the poles of $G^{(1)}$ at $\omega=-g_{1}({\bf k})$, where $G^{(2)}$ and $G^{(3)}$ approximately vanish for on-shell momenta. In this regime we find:

[TABLE]

which leads to the spectral function

[TABLE]

Therefore, at low energy the weight of the spectral function is concentrated in a limited part of the Brillouin zone. The zero energy states are obtained from the Green’s function by setting $\omega=0$ and therefore $k_{x}=0$ (we ignore the case of $k_{x}=\pi$ as there is vanishing spectral weight there). Since we have set $t_{s}=t^{\prime}=1$, this gives the condition:

[TABLE]

The left hand side of the above expression is zero when $k_{z}$ is at the Weyl points and is negative when $k_{z}$ is between them. Therefore the zero energy states reside on a straight line between the two Weyl point projections on the surface Brillouin zone.

III Application of a magnetic field

When put in magnetic field the density of states of metals oscillates as a function of inverse field. The oscillation frequency is proportional to the area enclosed by the Fermi surface cross section. This effect has been long utilized for characterization of materials. In two dimensions these quantum oscillations represent the Landau level quantization and in three dimensions the levels broaden and even overlap due to the dispersion along the field direction. Semiclassically the the oscillations can be viewed as follows. In the presence of magnetic field the quasiparticles encircle the Fermi surface due to the Lorentz force in momentum space. This closed orbit produces a maximum in the density of states when the quasiparticle phase accumulated during the motion is an integer times $2\pi$.

In a Weyl semimetal the Weyl points give rise to bulk quantum oscillations in any field direction. Importantly, the broadened Landau level include a zero energy level for each Weyl point which disperses linearly with the momentum along the field. We refer to this level as the chiral level and later comment about possible gapping of this level due to inter-nodal scattering. On the other hand, the surface states do not exhibit closed Fermi surfaces and quantum oscillations do not appear in the usual way. Instead, a path which includes the arcs on the top and bottom surfaces as well as tunneling through the low energy bulk states was proposed by Potter et al.Potter, Kimchi, and Vishwanath (2014) This semiclassical argument leads to surface level quantization of the form found in Eq. 2.

These semiclassical orbits require the use of low energy bulk modes near the Weyl points and therefore the bulk chiral Landau level is replaced by the fine quantization above. Bulk like levels are broader and begin at energies corresponding to $n>1$ bulk Landau level. Their quantization condition is different from that of the surface. We discuss both types of oscillations in this section.

III.1 Bulk Landau levels of a WSM

The bulk Landau levels can be obtained by considering the continuum low energy Weyl Hamiltonian, Eq. (1). In our case ${\bf v}=0$ and $h_{ij}$ is a diagonal matrix. Therefore,

[TABLE]

where $v_{z}=a\sqrt{t_{s}^{2}-m^{2}}$, $v_{F}=t_{s}a$ and $\chi$ is the chirality of the node. Here, $a$ is the lattice constant. We apply a magnetic field in the $y$-direction via substituting the canonical momentum $\mathbf{\pi}={\bf k}+{e\over\hbar c}\mathbf{A}$ and choosing the Landau gauge, $\mathbf{A}=-Bxe_{z}$. Defining raising and lowering operators gives:

[TABLE]

where $\ell_{B}=\sqrt{\hbar c/eB}$ is the magnetic length. The spectrum is found by squaring the Hamiltonian:

[TABLE]

where we set $\hbar=1$ from now on. Clearly the eigenstates of $H$ are eigenstates of $\sigma_{y}$. With the ansatz $\psi_{0}=(|0\rangle,i\chi|0\rangle)^{T}$, the zeroth Landau level dispersion is found to be $E_{0}=\chi v_{F}k_{y}$. Thus, we end up with a single chiral Landau level, independent of the magnetic field. Higher Landau levels are given by

[TABLE]

Note that each Weyl node exhibits only one zeroth Landau level. As explained by the Nielsen-Ninomiya theoremNielsen and Ninomiya (1981); Friedan (1982), Weyl nodes always come in pairs with opposite chiralities and in the full lattice model the chiral Landau levels are connected at high energy.

III.2 Surface quantum oscillations

III.2.1 Frequency

We now turn to the quantum oscillations associated with the surface Fermi arcs and analyze them using our effective surface theory. This allows us to test the semiclassical quantization condition of Eq. 2 and extend it beyond the decoupled Weyl points regime. The decoupled Weyl point approximation is valid for arc length $k_{a}$ such that $1/k_{a}\gg\ell_{B}$ and we therefore calculate the spectrum in and out of this regime. This is done by numerically solving the effective Green’s function equation 11.

In Fig. 2 (a) we show the surface density of states (DOS) calculated using the effective Greens function Eq. (8) for varying arc lengths in the long arc length regime for a fixed magnetic field and fixed slab thickness. In these graphs we represent the magnetic field by $q$, the number of unit cells in a magnetic unit cell. The magnetic unit cell is chosen such that the flux threading it is the flux quantum $\Phi_{0}$. This analysis produces Landau level separation which is in good agreement with Eq. 2. In order to further test the compatibility of the effective surface model with the semiclassical theory we repeat this analysis in Fig. 2(b) for varying slab thickness. In these cases we find that the observed oscillations and the effective surface theory coincide well with the semiclassical theory. For higher energies the clear oscillatory behaviour gets obscured by the fact that the bulk Landau levels overlap with the Fermi arc oscillations.

III.2.2 Phase offset

Another physical quantity that can be extracted from the plots is the phase offset $\gamma$. The graphs suggest that there is no Berry phase contribution to the surface Landau levels as $\gamma=\frac{1}{2}$. This can be interpreted as the cancellation of the Berry phase contributions of the two Weyl nodes. This is consistent with both chiral Landau levels participating in producing the surface Landau levels.

Moving from long arc lengths limit towards the small arc length regime, we expect hybridization between the two chiral Landau levels. This results in a gap and therefore a contribution to the phase offset. Indeed, when analyzing the energy offset in Fig. 3 (a) one can see that for short arc lengths a gap opens between the positive and negative surface Landau levels and the quantum oscillations vanish when approaching the point where the two Weyl nodes fuse at $m=t_{s}$. We find that the offset $\gamma$ in the short arc length regime originates from gapping of the chiral Landau level due to Weyl point mixing. We estimate the hybridization energy and compare it with the surface Landau level spectrum gap. This can be seen in Fig. 3.

In the short arc length limit the hybridization between two chiral Landau levels can be analyzed as tunneling in a double well potential in momentum space. In the appendix we analyze this problem using the WKB approximation and find the following energy splitting:

[TABLE]

where

[TABLE]

where we have explicitly included the lattice constant $a$. The Fermi velocity $v_{F}$, the effective mass $m^{*}$ for our model are defined in the appendix and $\omega_{c}$ is the cyclotron frequency. This offset is derived from the bulk model and in Fig. 3 (b) it is compared to the value of the zeroth Landau level at $k_{y}=0$ for the full bulk lattice model. One can see that the construction overestimates the energy gap by a small amount. This is due to the fact that the WKB approximation used in the derivation works better for higher Landau levels and the fact that we ignored a linear term in the potential. Nonetheless, the approximation captures the behaviour of the full system well.

Another interesting regime is a type II Weyl semimetal. In this regime, due to the tilt term there is no chiral Landau level while other bulk Landau levels are present. As suggested by the semiclassical analysis, the existence of surface Landau levels depends crucially on the bulk chiral Landau levelSoluyanov et al. (2015). In its absence we do not expect to see surface level quantization.

To test whether the surface states are quantized in Landau levels in a type II WSM we add a term to the lattice Hamiltonian in Eq. 5

[TABLE]

which turns the Weyl nodes into type-II Weyl nodes for $t_{2}>t_{s}$. In this regime $H_{2}$ completely dominates the spectrum and applying a magnetic field in a direction perpendicular to $z$ leads to a gapped spectrum. In Fig. 4 we show the results for various values of $t_{2}$, where one can see that when increasing $t_{2}$ the low energy regime which is dominated by the Fermi arc quantum oscillations shrinks until it completely vanishes for $t_{2}>t_{s}$. Our results therefore support the claim that the chiral Landau level is necessary for the formation of surface Landau levels.

IV Three Dimensional Lattice Model Analysis

In order to test the predictions of our effective surface theory, we use numerical diagonalization of a lattice model with an applied magnetic field via Peierls substitution. As in previous sections, the magnetic field $B$ is oriented in the $y$-direction and in the Landau gauge. In this gauge the hopping along $z$ acquires an $x$-dependent phase which breaks the translation invariance in the $x$-direction. We therefore define a magnetic unit cell, elongated along the $x$-direction. Choosing a cell of length $q$ lattice constants through which a flux quantum $\Phi_{0}$ is threaded amounts to a magnetic field $B=\Phi_{0}/qa^{2}$ where $a$ is the lattice constant. We vary $q$ to control the field strength. With this gauge the hopping along $z$ acquires a phase of $\exp\left(-i\frac{2\pi n_{x}}{q}\right)$, where $n_{x}$ is the index of the $n$th lattice site inside the magnetic unit cell. This increases the sizes of the matrices $H_{0}$ and $R$ to $2q$-by-$2q$ and the $y$-layer Hamiltonian reads

[TABLE]

where we have suppressed the $k$-indices on the creation/annihilation operators. The coupling between different $y$-layers is given by

[TABLE]

where ${\rm 1}_{q}$ is a $q$-by-$q$ unit matrix in the magnetic unit cell basis. The full $2qN_{y}$-by-$2qN_{y}$-Hamiltonian can be easily constructed. When diagonalizing the full system, we expect the low energy spectrum to be dominated by the surface quantum oscillations and by bulk contributions at higher energies.

In Fig. 5 (a) we show the low energy spectrum of a WSM slab in magnetic field. At low energy (below $0.6t_{s}$ in our model) one can see the Fermi-arc Landau level structure. For higher energies, one can observe bulk Landau levels mixed with the surface Landau levels. This is shown explicitly by varying the slab thickness. The frequency of the Fermi arc induced quantum oscillations varies with slab thickness, but the bulk Landau level spacing is approximately constant.

Analyzing the frequency of the the Landau levels enables us to test the predictions of the semiclassical theory and in addition further investigate the behaviour of the slab in regimes not accessible by the semiclassical theory. This is done in Fig. 5(b) where we confirm that the semiclassical approach accurately describes the frequency of the quantum oscillations when varying the slab width. Further, in Fig. 5(c) we vary the arc length $k_{a}=2\arccos{(\frac{m}{t_{s}})}$ by varying $m/t_{s}$ and compare the energy difference of the first and zeroth Landau level with the semiclassical theory Eq. (2). One can see that only for long arc lengths does the semiclassical theory describe the full quantum model well. This was done by taking into account the diabatic correctionPotter, Kimchi, and Vishwanath (2014) to the arc length $\propto\ell_{b}^{-1}$ which is due to the fact that an electron on the arc can tunnel through the bulk even before reaching a Weyl point. For small arc lengths the behaviour significantly deviates from the semiclassical theory and converges towards the behaviour of the merged Weyl points derived in appendix A.

Further, we can analyze the phase offset through the full bulk system and compare it to the findings of the effective surface theory. In Fig. 5(d) we show a comparison of the energy offset in the full bulk system with the semiclassical expectation with $\gamma=\frac{1}{2}$ as well as the tunnelling gap. One can see the expected crossover behaviour from the offset governed by the surface contribution to the bulk behaviour for $m\rightarrow t_{s}$.

V Conclusion

In this work we derived an effective surface theory for the surface of a Weyl semimetal. We used this theory as well as the full three dimensional lattice model to study the quantum oscillations. We find that the surface theory accurately predicts the quantum oscillations which are associated with the surface states. Our model allows us to find the regime of validity of the semiclassical analysis of the surface quantum oscillations. We find that it fails when the Weyl points hybridize. This hybridization gaps the bulk chiral Landau level and in turn gaps the surface quantum oscillations which are essentially a fine quantization inside the chiral Landau level. We estimate the hybridization gap using a double well potential model analyzed in the WKB approximation and find good agreement with the three dimensional model.

Appendix A Bulk Landau levels for overlapping Weyl nodes

For large Weyl node separations, one can model each Weyl node separately and arrive at the usual Weyl Landau levels. However, when the two Weyl nodes approach each other, they can not be treated separately. As a result the pair of chiral Landau levels hybridize and the spectrum developes a gap at zero energy.

In order to address this we consider a bulk model with an applied magnetic field in the $y$-direction, $\vec{B}=B\vec{e}_{y}$ and vector potential $\vec{A}=Bz\vec{e}_{x}$. As discussed in the introduction, the system experiences a gap closure at $m=t_{s}$. Lowering $m$ further, the gap closure splits into two Weyl nodes which traverse the Brillouin zone and recombine for $m=-t_{s}$. The point of the gap closure for $m=t_{s}$ is $\vec{k}_{0}=(0,0,-\pi/2)^{T}$. Close to this point one can expand the Hamiltonian as

[TABLE]

where $\gamma=m-t_{s}$ and we have omitted the quadratic terms in $k_{x}$ and $k_{y}$. We defined the Fermi velocity $v_{F}=ta$ and $m^{*}=\frac{1}{a^{2}t_{s}}$, where we explicitly wrote out the lattice constant $a$. Adding the magnetic field to the system, we arrive at

[TABLE]

where we define $\tilde{z}=z+k_{x}/B$. Normally one arrives at a zeroth Landau level which depends on the momentum along the field direction, in this case $k_{y}$, and is gapless for $k_{y}=0$. Here, however, the zeroth Landau level acquires a mass which is estimated below.

Defining the new variables $Z=\frac{\tilde{z}}{\alpha}$ and $K=\alpha k_{z}$ with $\alpha=(v_{F}m^{*}B)^{-\frac{1}{3}}$, we arrive at

[TABLE]

with the now dimensionless $\Gamma=\frac{\gamma}{(m^{*}v_{F}^{2}\omega_{c}^{2})^{\frac{1}{3}}}$ and the cyclotron frequency $\omega_{c}=\frac{B}{m^{*}}$. Squaring the Hamiltonian leads to

[TABLE]

where we have used $[Z,K^{2}]=2iK$. Ignoring the prefactor for now, one gets the eigenvalue equation,

[TABLE]

Using the ansatz $c_{2}=ic_{1}$, we get

[TABLE]

This has the general structure of the differential equation for an anharmonic oscillator when the operators are written in the momentum space basis as opposed to the real space basis. Using this analogy, we can solve this by defining the ’potential’ in momentum space, $V\equiv V(K)=(\Gamma+\frac{K^{2}}{2})^{2}+K$. In order to analyze the equation, we approximate it by ignoring the linear term. We want to investigate the $m=t_{s}$ limit, i.e., $\Gamma=0$. In this case the two Weyl nodes are combined to one single gap closure and one can solve the problem via the WKB approximation. The quantization condition reads,

[TABLE]

where $x_{\pm}=\pm\sqrt{2}E^{\frac{1}{2}}$ are the turning points. The left hand side can be transformed,

[TABLE]

with $t=\left(\frac{K}{x}\right)^{4}$. The integral is defined as the Euler beta function,

[TABLE]

at $x=\frac{1}{4}$ and $y=\frac{3}{2}$ and therefore,

[TABLE]

With $\mathcal{B}\left(\frac{1}{4},\frac{3}{2}\right)=\frac{\sqrt{\pi}\Gamma(\frac{1}{4})}{2\Gamma(\frac{7}{4})}$ we get

[TABLE]

and reinstating the prefactor leads to

[TABLE]

As usual, the WKB approximation works better for higher Landau levels and we have ignored the linear term. When comparing to the numerically calculated Landau levels of the full system, only the first couple of Landau levels acquire a correction $\zeta(n)$. We find numerically that the zeroth Landau level gets a factor of $\zeta(0)\approx 0.811$ and already for $n\geq 1$ we find $\zeta(n)\approx 1$, which is very close to what the authors in [Dietl, Piéchon, and Montambaux, 2008] have found.

Appendix B Chiral Landau level for small Weyl node separation

In order to estimate the energy shift of the chiral Landau level as a function of the Weyl node separation, we will start from equation (31), from which one can read off the potential in momentum space as $V(K)=(\Gamma+\frac{K^{2}}{2})^{2}+K$. This is the potential of an anharmonic oscillator in momentum space and for $\Gamma<0$ we have two distinct minima. Ignoring the linear term, we have a symmetric double well problemLandau and Lifshitz (1965), where $\Gamma=\frac{m-t_{s}}{(m^{*}v_{F}^{2}\omega_{c}^{2})^{\frac{1}{3}}}$ controls the separation of the two wells. We approach the problem by restricting our Hilbert space to that spanned by two wavefunctions, $\psi_{r}$ and $\psi_{l}$. These are the ground states of the left and right wells located at $\pm\sqrt{2\Gamma}$, when they are completely separated form each other. When the two wells approach each other, the two wavefunctions hybridize due to the finite tunnelling probability. The new eigenstates are the symmetric and anti-symmetric combinations,

[TABLE]

The wavefunctions obey the Schrödinger equations,

[TABLE]

where the first equation is valid due to the fact that the double well potential equals the single well potential in the regime of $\psi_{r}$. In addition, the amplitude of $\psi_{r}$ in the left well is vanishingly small. Without loss of generality, we pick $\psi_{+}$. Multiplying the first equation by $\psi_{+}$ and the second one by $\psi_{r}$, subtracting the first equation from the second and subsequently integrating from [math] to $\infty$ and using integration by parts, one arrives at,

[TABLE]

Restoring the prefactor from (31) and taking into account that we squared the Hamiltonian in order to derive Eq. (31), via $\Delta\epsilon=(m^{*}v_{s}^{2}\omega_{c}^{2})^{\frac{1}{3}}\sqrt{\Delta E}$ we arrive at,

[TABLE]

In order to evaluate this expression, we use the WKB approximation and we have,

[TABLE]

Evaluating the integral we arrive at,

[TABLE]

The distance of the Weyl node from the middle point around which we have expanded the Hamiltonian, is given by $k_{a}/2=\arccos(m/t_{s})\approx\sqrt{2\frac{(t_{s}-m)}{t_{s}}}$, where we have expanded around $m=1$. With this, we arrive at the final result,

[TABLE]

The constant $C$ can be inferred from the lowest Landau level at $m=t_{s}$ (see Appendix A) and we get

[TABLE]