Topologically distinct Weyl-fermion pairs: on the effect of magnetic tunnelling

TL;DR

This paper explores how different topological configurations of Weyl-fermion pairs affect their electronic properties under magnetic fields, revealing that vorticity, not chirality, governs zero-energy modes and other phenomena.

Contribution

It identifies the impact of Weyl node vorticity on Landau levels and topological features, distinguishing between different connection topologies of Weyl nodes.

Findings

Zero-energy Landau levels are linked to vorticity, not chirality.

Magnetic tunneling does not always remove zero-energy modes.

Disorder effects and surface states vary with topological configuration.

Abstract

A Weyl semimetal has Weyl nodes that always come in pairs with opposite chiralities. Notably, different ways of connection between nodes are possible and would lead to distinct topologies. Here we identify their differences in many respects from two proposed models with different vorticities. One prominent feature is the behaviour of zeroth Landau levels (LLs) under magnetic field. We demonstrate that the magnetic tunnelling does not always expel LLs from zero energy because the number of zero-energy modes is linked to the vorticity of the Weyl nodes, instead of the chirality. Other respects in disorder effects for weak (anti-)localization, surface Fermi arcs, and Weyl-node annihilation, are interesting consequences that await future exploration.

Tables

Table 1. Table 1: Average inter-valley scattering rates ⟨ 1 τ I ⟩ delimited-⟨⟩ 1 subscript 𝜏 𝐼 \langle\frac{1}{\tau_{I}}\rangle for anisotropic impurities potentials. Here p 𝑝 p -wave ( p x subscript 𝑝 𝑥 p_{x} , p y subscript 𝑝 𝑦 p_{y} , and p z subscript 𝑝 𝑧 p_{z} ) impurities potentials are considered. The p x subscript 𝑝 𝑥 p_{x} -impurity potential (second row) does not differentiate the two models, but the form of ∼ p ∥ = p y similar-to absent subscript 𝑝 parallel-to subscript 𝑝 𝑦 \sim p_{\parallel}=p_{y} or p z subscript 𝑝 𝑧 p_{z} can tell the difference (last row). N F subscript 𝑁 𝐹 N_{F} is the density of states at the Fermi energy and u I subscript 𝑢 𝐼 u_{I} characterizes the strength of the impurity potential.

impurity potential	Model A	Model B
$U_{q,q^{\prime }}^{+-}$	$(+,+,+)\to (-,+,+)$	$(+,+,+)\to (-,-,-)$
$u_I\frac{\left(q_x-q_x^{\prime }\right)}{k_F}$	$⟨\frac{1}{{\tau }_I}⟩=\frac{2\pi }{\hslash }{N}_F\frac{8}{9}{u}_I^2$	$⟨\frac{1}{{\tau }_I}⟩=\frac{2\pi }{\hslash }{N}_F\frac{8}{9}{u}_I^2$
$u_I\frac{\left(q_{\parallel }-q_{\parallel }^{\prime }\right)}{k_F}$	$⟨\frac{1}{{\tau }_I}⟩=\frac{2\pi }{\hslash }{N}_F\frac{4}{9}{u}_I^2$	$⟨\frac{1}{{\tau }_I}⟩=\frac{2\pi }{\hslash }{N}_F\frac{8}{9}{u}_I^2$

Equations

$${\mathcal{M}}_{x}H(k_x,k_y,k_z){\mathcal{M}}_x^{-1}=H(-k_x,k_y,k_z),$$

$$H_{\mathrm{A}}=\frac{1}{2m}\left(k_x^2-k_{\mathrm{W}}^2-\alpha k_{\parallel }^2\right){\sigma }_x+v_{\parallel }{k}_y{\sigma }_y+v_{\parallel }{k}_z{\sigma }_z,$$

$$H_{\mathrm{B}}=\frac{1}{2m}\left(k_x^2-k_{\mathrm{W}}^2-\alpha k_{\parallel }^2\right){\sigma }_x+\frac{v_{\parallel }}{k_{\mathrm{W}}}{k}_x\left(k_y{\sigma }_y+k_z{\sigma }_z\right),$$

$$H_{\mathrm{A}}$$

$$H_{\mathrm{B}}$$

$$k_x{k}_y\to \frac{1}{2}\left({\Pi }_x{\Pi }_y+{\Pi }_y{\Pi }_x\right)=l_B^{-2}\left(k_x\bar{x}+\frac{i}{2}\right),$$

$$H_{\mathrm{B}}^{\prime}=\frac{\omega_{c}^{2}}{4E_{\mathrm{VH}}}\left\{\begin{array}[]{l}\left[\frac{1}{g}(q^{2}-1)+g\frac{\tilde{\alpha}}{4}\frac{\partial^{2}}{\partial q^{2}}-g\frac{\tilde{\alpha}}{4}q_{z}^{2}\right]\sigma_{x}\\ +i\left(q\frac{\partial}{\partial q}+\frac{1}{2}\right)\sigma_{y}+q_{z}q\sigma_{z}\end{array}\right\},$$

$$⟨\frac{1}{{\tau }_I}⟩=\frac{2\pi }{\hslash }\sum _{q,q^{\prime }}⟨|U_{q{q}^{\prime }}^{+-}{|}^2⟩\delta (E_F-{\xi }_{q^{\prime }})\delta (E_F-{\xi }_q)$$

$$\sqrt{\frac{2}{g}}q=\frac{a+a^{†}}{\sqrt{2}},\text{ and }\sqrt{\frac{g}{2}}\frac{\partial }{\partial q}=\frac{a-a^{†}}{\sqrt{2}}$$

$$a^{\dagger}=\left(\begin{array}[]{ccccc}0&0&0&\dots&0\\ 1&0&0&\dots&0\\ 0&\sqrt{2}&0&\dots&0\\ \vdots&\vdots&\ddots&\ddots&0\\ 0&0&\dots&\sqrt{L-1}&0\\ \end{array}\right)$$

$$\psi=\left(\begin{array}[]{c}\tilde{\phi}\\ 0\\ \vdots\\ 0\\ \tilde{\chi}\\ 0\\ \vdots\\ 0\end{array}\right),$$

$$H=\left(\begin{array}[]{ccccc|ccccc}&&&\times&\times&&&&\times&\times\\ \lx@intercol\hfil\scalebox{1.7}{$\tilde{H}_{11}$}\hfil\lx@intercol&\vdots&\vdots&\lx@intercol\hfil\scalebox{1.7}{$\tilde{H}_{12}$}\hfil\lx@intercol&\vdots&\vdots\\ &&&\times&\times&&&&\times&\times\\ \hline\cr&&&\times&\times&&&&\times&\times\\ \lx@intercol\hfil\scalebox{1.7}{$\tilde{H}_{21}$}\hfil\lx@intercol&\vdots&\vdots&\lx@intercol\hfil\scalebox{1.7}{$\tilde{H}_{22}$}\hfil\lx@intercol&\vdots&\vdots\\ &&&\times&\times&&&&\times&\times\\ \end{array}\right),$$

$$\tilde{H}=\left(\begin{array}[]{c|c}\tilde{H}_{11}&\tilde{H}_{12}\\ \hline\cr\tilde{H}_{21}&\tilde{H}_{22}\\ \end{array}\right),$$

$$\tilde{\psi}=\left(\begin{array}[]{c}\tilde{\phi}\\ \tilde{\chi}\end{array}\right)$$

$$H^{\dagger}=\left(\begin{array}[]{cccc|cccc}&&&&&&&\\ \lx@intercol\hfil\raisebox{6.0pt}{ \scalebox{1.7}{$\tilde{H}_{11}^{\dagger}$} }\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\raisebox{6.0pt}{ \scalebox{1.7}{$\tilde{H}_{12}^{\dagger}$} }\hfil\lx@intercol\\ \times&\cdots&\cdots&\times&\times&\cdots&\cdots&\times\\ \times&\cdots&\cdots&\times&\times&\cdots&\cdots&\times\\ \hline\cr&&&&&&&\\ \lx@intercol\hfil\raisebox{6.0pt}{ \scalebox{1.7}{$\tilde{H}_{21}^{\dagger}$} }\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\raisebox{6.0pt}{ \scalebox{1.7}{$\tilde{H}_{22}^{\dagger}$} }\hfil\lx@intercol\\ \times&\cdots&\cdots&\times&\times&\cdots&\cdots&\times\\ \times&\cdots&\cdots&\times&\times&\cdots&\cdots&\times\\ \end{array}\right),$$

$$H^{\dagger}H=\left(\begin{array}[]{ccccc|ccccc}&&&\times&\times&&&&\times&\times\\ \lx@intercol\hfil\raisebox{6.0pt}{ \scalebox{1.7}{$B_{11}$} }\hfil\lx@intercol&\vdots&\vdots&\lx@intercol\hfil\raisebox{6.0pt}{ \scalebox{1.7}{$B_{12}$} }\hfil\lx@intercol&\vdots&\vdots\\ \times&\cdots&\times&\times&\times&\times&\cdots&\times&\times&\times\\ \times&\cdots&\times&\times&\times&\times&\cdots&\times&\times&\times\\ \hline\cr&&&\times&\times&&&&\times&\times\\ \lx@intercol\hfil\raisebox{6.0pt}{ \scalebox{1.7}{$B_{21}$} }\hfil\lx@intercol&\vdots&\vdots&\lx@intercol\hfil\raisebox{6.0pt}{ \scalebox{1.7}{$B_{22}$} }\hfil\lx@intercol&\vdots&\vdots\\ \times&\cdots&\times&\times&\times&\times&\cdots&\times&\times&\times\\ \times&\cdots&\times&\times&\times&\times&\cdots&\times&\times&\times\\ \end{array}\right),$$

$$B\equiv\tilde{H}^{\dagger}\tilde{H}=\left(\begin{array}[]{cc}B_{11}&B_{12}\\ B_{21}&B_{22}\\ \end{array}\right).$$

$$\begin{array}{r}{H}_{\mathrm{B}}^{\prime }=\frac{{\omega }_c^2}{4E_{\mathrm{VH}}}\{\left[\frac{1}{g}(q^2-1)+g\frac{\tilde{\alpha }}{4}\frac{{\partial }^2}{\partial q^2}-g\frac{\tilde{\alpha }}{4}{q}_z^2\right]{\sigma }_x\\ +i\left(q\frac{\partial }{\partial q}+\frac{1}{2}\right){\sigma }_y+q_{z}q{\sigma }_z\},\end{array}$$

$$\begin{array}{r}\left\{-\left(q\frac{\partial }{\partial q}+\frac{1}{2}\right)\left(q\frac{\partial }{\partial q}+\frac{1}{2}\right)+\frac{2}{g}{q}^2+\frac{1}{g^2}(q^2-1{)}^2\right\}\chi (q)\\ ={\epsilon }^2\chi (q)\\ \left\{-\left(q\frac{\partial }{\partial q}+\frac{1}{2}\right)\left(q\frac{\partial }{\partial q}+\frac{1}{2}\right)-\frac{2}{g}{q}^2+\frac{1}{g^2}(q^2-1{)}^2\right\}\Phi (q)\\ ={\epsilon }^2\Phi (q)\end{array}$$

$$\begin{array}{rl}\left(q\frac{\partial }{\partial q}+\frac{1}{2}\right)\frac{\chi }{\sqrt{q}}=\frac{1}{\sqrt{q}}\left(q\frac{\partial }{\partial q}\right)\chi ,& \\ & \text{ and thus }\\ {\left(q\frac{\partial }{\partial q}+\frac{1}{2}\right)}^2\frac{\chi }{\sqrt{q}}=\frac{1}{\sqrt{q}}{\left(q\frac{\partial }{\partial q}\right)}^2\chi ,& \end{array}$$

$$\begin{array}[]{l}\left[q^{2}\frac{\partial^{2}}{\partial q^{2}}+\left(1-\frac{2}{g}q^{2}\right)q\frac{\partial}{\partial q}+\frac{2}{g^{2}}(1-2g)q^{2}\right]\tilde{\chi}(q)=\lambda^{2}\tilde{\chi}(q),\\ \left[q^{2}\frac{\partial^{2}}{\partial q^{2}}+\left(1-\frac{2}{g}q^{2}\right)q\frac{\partial}{\partial q}+\frac{2}{g^{2}}q^{2}\right]\tilde{\Phi}(q)=\lambda^{2}\tilde{\Phi}(q),\\ \end{array}$$

$$\begin{array}[]{l}q^{2}\frac{\partial^{2}f_{1}}{\partial q^{2}}+\left(1+2\lambda-\frac{2}{g}q^{2}\right)q\frac{\partial f_{1}}{\partial q}+\frac{2}{g^{2}}(1-2g-g\lambda)q^{2}f_{1}=0,\\ q^{2}\frac{\partial^{2}f_{2}}{\partial q^{2}}+\left(1+2\lambda-\frac{2}{g}q^{2}\right)q\frac{\partial f_{2}}{\partial q}+\frac{2}{g^{2}}(1-g\lambda)q^{2}f_{2}=0.\\ \end{array}$$

$$\begin{array}[]{l}\bar{q}^{2}\frac{\partial^{2}f_{1}}{\partial\bar{q}^{2}}+\left(1+2\lambda-2\bar{q}^{2}\right)\bar{q}\frac{\partial f_{1}}{\partial\bar{q}}+\frac{2}{g}(1-2g-g\lambda)\bar{q}^{2}f_{1}=0,\\ \bar{q}^{2}\frac{\partial^{2}f_{2}}{\partial\bar{q}^{2}}+\left(1+2\lambda-2\bar{q}^{2}\right)\bar{q}\frac{\partial f_{2}}{\partial\bar{q}}+\frac{2}{g}(1-g\lambda)\bar{q}^{2}f_{2}=0.\\ \end{array}$$

$$\begin{array}[]{l}4\rho^{2}\frac{\partial^{2}f_{1}}{\partial\rho^{2}}+4(1+\lambda-\rho)\rho\frac{\partial f_{1}}{\partial\rho}+\frac{2}{g}(1-2g-g\lambda)\rho f_{1}=0\\ 4\rho^{2}\frac{\partial^{2}f_{2}}{\partial\rho^{2}}+4(1+\lambda-\rho)\rho\frac{\partial f_{2}}{\partial\rho}+\frac{2}{g}(1-g\lambda)\rho f_{2}=0\\ \end{array}$$

$$\begin{array}[]{ccl}\bar{q}\frac{df}{d\bar{q}}&=&2\rho\frac{df}{d\rho},\\ \bar{q}^{2}\frac{d^{2}f}{d\bar{q}^{2}}&=&4\rho^{2}\frac{d^{2}f}{d\rho^{2}}+2\rho\frac{df}{d\rho}.\end{array}$$

$$\begin{array}[]{l}\rho\frac{\partial^{2}f_{1}}{\partial\rho^{2}}+(\lambda+1-\rho)\frac{\partial f_{1}}{\partial\rho}+(n-1)f_{1}=0,\\ \rho\frac{\partial^{2}f_{2}}{\partial\rho^{2}}+(\lambda+1-\rho)\frac{\partial f_{2}}{\partial\rho}+nf_{2}=0,\\ \end{array}$$

$$E_n=2\frac{{\omega }_c^2}{4E_{\mathrm{VH}}}{\left[n\left(\frac{1}{g}-n\right)\right]}^{\frac{1}{2}}.$$

$$\Phi (q)=\frac{1}{\sqrt{q}}{e}^{-\frac{q^2}{2g}}{q}^{\lambda }{L}_n^{\lambda }\left(\frac{q^2}{g}\right)=q^{\lambda -\frac{1}{2}}{L}_n^{\lambda }\left(\frac{q^2}{g}\right)e^{-\frac{q^2}{2g}},$$

$$\frac{1}{g}(q^2-1)+g\frac{\tilde{\alpha }}{4}\frac{{\partial }^2}{\partial q^2}\pm \left(q\frac{\partial }{\partial q}+\frac{1}{2}\right).$$

$$\left[g\frac{\tilde{\alpha }}{4}\frac{{\partial }^2}{\partial q^2}+\left(q\frac{\partial }{\partial q}+\frac{1}{2}\right)+\frac{1}{g}(q^2-1)\right]\psi (q)=0.$$

$$+\frac{g\tilde{\alpha }}{4}\frac{{\partial }^2\Phi }{\partial q^2}-\sqrt{1-\tilde{\alpha }}q\frac{\partial \Phi }{\partial q}-[\frac{1}{g}+\frac{\sqrt{1-\tilde{\alpha }}}{2}]\Phi =0.$$

Full text

Introduction.— The prediction of Weyl semimetals (WSMs) Wan et al. (2011); Burkov and Balents (2011); Halász and Balents (2012); Xu et al. (2011); Liu and Vanderbilt (2014); M. Turner and Vishwanath (2013) and its realization in real materials Huang et al. (2015a); Xu et al. (2015a); Weng et al. (2015); Lv et al. (2015a); Yang et al. (2015); Xu et al. (2015b); Lv et al. (2015b) make the relativistic chiral fermions find their counterpart in three-dimensional (3D) condensed matter systems. These chiral fermions reside in the nodes of the electronic structure around which the energy dispersion is linear in momentum, the so-called Weyl nodes (WNs). The WNs always come in pairs with opposite chiralities Nielsen and Ninomiya (1981, 1983a) which act as the source and the sink of the Berry curvature Yan and Felser (2017). Due to topological reasons, various unusual behaviours were found in the WSM, such as the chiral anomaly Adler (1969); Nielsen and Ninomiya (1983b); Zyuzin and Burkov (2012); Parameswaran et al. (2014); Jia et al. (2016), negative magnetoresistance Son and Spivak (2013); Huang et al. (2015b), chiral magnetic effect Landsteiner (2014), novel quantum oscillation Hosur (2012); Potter et al. (2014); Zhang et al. (2016); Bulmash and Qi (2016), Fermi arc from the surface states Wan et al. (2011); Lv et al. (2015c); Jia et al. (2016); Batabyal et al. (2016), and so on.

One prominent and important phenomena of the WN is the behaviour of Landau levels (LLs) under magnetic field, which were mapped out by magneto-optical study recently in NbAs Yuan et al. (2018). The zeroth LL ($n=0$), also called chiral Landau band, has linear dispersion along the field direction, say $\hat{z}$, as $E=\chi v_{z}k_{z}$. Interestingly, the chirality $\chi=\pm 1$ determines the sign of the band slope. In particular, the connection of opposite chiral bands by field applied along the WNs provides the platform of chiral anomaly in which the charge pumping breaks the chiral symmetry. This actually relies on the existence of the zero-energy modes at $k_{z}=0$ which are topologically protected, or otherwise the system becomes insulating and the charge pumping may be forbidden.

Semiclassically the LLs are formed through quantization conditions of cyclotron orbits. When cyclotron orbits encounter each other, go across density discontinuities Kane et al. (1988), or are blocked by some boundaries Goran et al. (2013), different quantization conditions may be formed and hence the LLs are changed. Besides, there exists tunnelling between cyclotron orbits, generally known as magnetic breakdown Shoenberg (1984). For example it has been discovered in adjecent quantum wells Goran et al. (2013). Particularly, inter-level tunnelling between levels from separate chiralities produce new features in transport in graphene Greenaway et al. (2015). Since level mixing is common, the question of interest now is whether the phenomenon happens in WSMs. As WNs are connected by bands, the magnetic tunnelling between LLs is expected when the field is applied perpendicular to the connection of WNs, possibly gapping the zeroth LL. If it is the case, this would also lead to failure of chiral anomaly, as indicated as the possible explanation for the increased magneto-resistance Zhang et al. (2017a); Ramshaw et al. (2018) or sharp sign reversal of Hall resistivity Zhang et al. (2017b).

However, the reverse inference might not be true, meaning that the magneto-resistance changes are not necessarily attributed to the gapping of zeroth LLs. They may also be caused by the gap opening in the system through multiple Weyl carrier interaction Ramshaw et al. (2018) or with the help of other non-Weyl singularities Kim et al. (2017). Besides, recently there are also concerns about interpreting the measured negative longitudinal magnetoresistance as direct evidence of chiral anomaly dos Reis et al. (2016); Li et al. (2017); Liang et al. (2018). Therefore, not only helping identify reasons of resistance changes affirmatively but also with theoretical significance, it is important and interesting to see if the gapping of zeroth LLs is an inevitable result. Actually, the above inference of repelling zero modes only considers the simple connection between WNs, while symmetry constraints could make situations change. For example, the mirror symmetry is the commonly seen constraint connecting the WNs Huang et al. (2015a); Xu et al. (2015a); Weng et al. (2015); Lv et al. (2015b); Gao et al. (2016).

To consider consequences from different ways of connection between WNs, we studied two models with different symmetry constraints imposed. Magnetic tunnelling was found to be different and the zero-energy modes are still robust in one model against the tunnelling. We attribute the findings to distinct topological invariants, suggesting that chirality alone is not sufficient to characterize a WN. Moreover, respective unique phenomena in impurity scattering, surface Fermi arcs, and WN annihilation are also studied and can differentiate the two topological distinct models. The mutual annihilation between WNs does not necessarily become gapped and a nodal ring is possible to be formed, consistent with the newly discovered conversion rule Sun et al. (2018). Therefore, the detailed ways of connection between WNs await future more investigation.

Models.— We consider a pair of WNs with opposite chiralities sitting on two sides of the mirror plane $M_{x}$. Their separation $2k_{\mathrm{W}}$ is relatively small compared to the size of the Brillouin zone (BZ) such that the magnetic field has a chance to couple them. Other WNs, if exist, in other regions of the BZ can be ignored as they are much far away. To have

[TABLE]

we find that there are two possible choice of $\mathcal{M}_{x}$ and then models. We dub them Model A and Model B: Model A is for $\mathcal{M}_{x}=\sigma_{0}$ that the two bands have equal mirror parity in the mirror plane, while Model B is for $\mathcal{M}_{x}=\sigma_{x}$ that two bands have opposite mirror parities. Specifically, the Hamiltonian are written as

[TABLE]

[TABLE]

where $k_{\parallel}^{2}=k_{y}^{2}+k_{z}^{2}$. In general, coefficient for $k_{y}$ and $k_{z}$ can be different, but the physics are the same. The Planck constant $\hbar$ is set as unity throughout the paper. We note that the origin is meaninglessly specified and might not be at a time-reversal-invariant momentum.

Both two models contain two WNs located at $(\pm k_{\mathrm{W}},0,0)$. Expanding around the WNs, they approximate, to linear order, as

[TABLE]

Here $\chi$ labels the chirality and also position of the WN, and $v_{x}=k_{\mathrm{W}}/m$ and $v_{y}=v_{z}=v_{\parallel}$. The values of $k_{\mathrm{W}}$, $m$, $v_{y}$ and $v_{z}$ are all assumed to be positive without loss of generosity. $\alpha$ ($>0$) in $H_{\mathrm{B}}$ is required in order to have a saddle point at energy $E_{\mathrm{VH}}=k_{\mathrm{W}}^{2}/2m$, above which close energy contours are assured. Although a $k_{\parallel}$-linear term is allowed to appear in the $\sigma_{x}$ term, it is omitted for an elegance reason. We have confirmed that its presence once being small does not change qualitative conclusions.

Without special regard to the symmetry or phase, Model A was usually adopted to study the effect of a pair of nodes Zhang et al. (2017b); Chan and Lee (2017). Applying the magnetic field along the perpendicular $z$ direction to this system, we solve the LL spectrum by substituting $\mathbf{k}$ in Eq. (2) by $\mathbf{\Pi}=\mathbf{k}+e\mathbf{A}$. We choose the Landau gauge $\mathbf{A}=Bx\hat{y}$, so we make $k_{x}\rightarrow\Pi_{x}=k_{x}$, $k_{z}\rightarrow\Pi_{z}=k_{z}$, and $k_{y}\rightarrow\Pi_{y}=l_{B}^{-2}\bar{x}$, where the magnetic length $l_{B}=\sqrt{\frac{1}{eB}}$ and $\bar{x}$ is the coordinate relative to the guiding center $x_{0}=-l_{B}^{2}k_{y}$. $\bar{x}$ is conjugate to $k_{x}$ by quantizing $\bar{x}\rightarrow i\frac{\partial}{\partial k_{x}}$. The WN separation is used as a measure to define the dimensionless momentum scale as $q=k_{x}/k_{\mathrm{W}}$. Since the magnetic field breaks the inplane translation symmetry, two WNs are expected to couple via the field. The coupling way can be revealed from the missing terms to Eq. (4). We use a dimensionless parameter $g$ to describe the degree of the coupling. The coupling will increase with the cyclotron energy $\omega_{c}=\sqrt{2v_{x}v_{y}}l_{B}^{-1}$ and decrease with the energy barrier $E_{\mathrm{VH}}$, defined as $g=\frac{\omega_{c}^{2}}{4E_{\mathrm{VH}}^{2}}$. $g$ is proportional to the magnetic field $B$, and the appreciable coupling $g\approx 1$ is achieved when the magnetic length $l_{B}$ is comparable to the scale defined by $k_{\mathrm{W}}^{-1}$ 111Take Weyl semimetal TaAs for example. The nodes separation of $W_{1}$ is $k_{\mathrm{W}}=0.0072\left(\frac{2\pi}{a}\right)$, and the Fermi velocity in the conduction band is $(v_{x},v_{y},v_{z})=(2.5,1.2,0.2)\times 10^{5}$ m/s Lee et al. (2015). The coupling $g=1$ under the field strength $B\approx 11.88$ T having magnetic length $l_{B}=7.44$ nm, which is close to the length 7.60 nm defined by $k_{\mathrm{W}}^{-1}$. By defining the remaining variables into dimensionless quantities, $q_{z}=v_{\parallel}k_{z}\left(\frac{\omega_{c}^{2}}{4E_{\mathrm{VH}}}\right)^{-1}=(\frac{2}{g})(\frac{v_{\parallel}}{v_{x}})(\frac{k_{z}}{k_{\mathrm{W}}})$, and $\tilde{\alpha}=\left(\frac{v_{x}}{v_{\parallel}}\right)^{2}\alpha$, we can study the Hamiltonian under magnetic field as a function of $q$ and $q_{z}$. We numerically solve this system with raising and lowering operators. Special treatment is developed to solve it more efficiently and the details are shown in the Supplemental Material Sup .

The LL spectrum at $k_{z}=0$, i.e. $q_{z}=0$, with respect to $g$ is shown in Fig. 1(a). In the limit of $g\rightarrow 0$, two WNs have independent and identical LLs, so each LL is doubly degenerate. As $g$ is turned on, the degeneracy is lifted off and band splits are visible at $g\approx 0.3$ (decrease with levels) in Fig. 1(a). We have ascribed the band splits to the magnetic tunnelling in Ref. Zhang et al. (2017b) as the cyclotron wave functions in $k$ space broaden with the $B$ field and hybridize with others when overlaps occur. Notice that the chiral symmetry is present for $\left\{H,\sigma_{z}\right\}=0$ at $k_{z}=0$, so the spectrum is symmetric about zero energy. In order to obey the chiral symmetry, the zeroth LL has to split into one with positive energy and one with negative energy.

The results shown in Model A are reasonable but not conclusive. We solve Model B as follows. In Model B, the Peierls substitution should be carefully treated in $k_{x}k_{y}$ and $k_{x}k_{z}$. To make the Hamiltonian Hermitian, we do the symmetrization

[TABLE]

and $k_{x}k_{z}\rightarrow(\Pi_{x}\Pi_{z}+\Pi_{z}\Pi_{x})/2=k_{x}k_{z}$. In terms of dimesionless parameters defined above, the Hamiltonian under magnetic field becomes

[TABLE]

where the prime stands for a system under a magnetic field.

The LL spectrum with respect to $g$ for Model B at $k_{z}=0$ is shown in Fig. 1 (b). For small $g$, the effect of $\tilde{\alpha}$ is small, so we provide the results for $\tilde{\alpha}=0$ as shown by red dots in Fig. 1 (b). The solutions are analytical to be $E_{n}=\frac{\omega_{c}^{2}}{2E_{\mathrm{VH}}}\sqrt{n\left(\frac{1}{g}-n\right)}$ with $n\in\mathbb{Z}_{\geq}0$ and each energy has twofold degeneracy. The LL energies emerge from the typical LL spectrum as $E_{n}=\sqrt{n}\omega_{c}$ in the limit of $g\rightarrow 0$ and deform with finite $g$. However, the analytical solutions for $\tilde{\alpha}=0$ are only applicable in a low-energy region. At high energies, with large $g$ or $n$, the LL quantization makes no sense when $E>E_{\mathrm{VH}}$ due to the flat dispersion along the $k_{y}$ axis and thus unbounded equal-energy contours. Specifically, the analytical solutions are applicable when $n\leq\frac{1}{2}\left(\frac{1}{g}-\frac{1}{2}\right)$.

Since the model is invariant under inversion in $q$, the eigenstates will be either even (e) or odd (o) in $q$ as denoted in Fig. 1 (b). The even- and odd states appear alternatively in energy, showing that they evolve from a degenerate spectrum for small $g$. We were unable to prove whether the degeneracy for $n>0$ is exact at small and finite $g$ but found it seeming to be in the applicable region of $\tilde{\alpha}=0$. The band splitting is reasonable as seen in a symmetric double-well with finite tunnelling probability where even and odd states have different energies. By contrast, Model A does not have this symmetry and therefore its eigenstates do not respect this symmetry in $q$.

Nevertheless, two zero-energy LL states persist for all values of $g$ and $\tilde{\alpha}$. They can be directly verified by solving the zero-eigenvalue problem. The eigenfunctions are found to be $\left(0,e^{-\kappa q^{2}}\Phi(q)\right)^{T}$, where $\Phi(q)$ can be either ${{}_{1}F_{1}}(-\frac{1}{4}\lambda;\frac{1}{2};\xi^{2}q^{2})$ or $H_{\lambda/2}(\xi q)$, the former being the Kummer confluent hypergeometric function and the latter the Hermite polynomial, indicating a double degeneracy. Here $\kappa=\frac{1+\sqrt{1-\tilde{\alpha}}}{g\tilde{\alpha}}$, $\xi=\sqrt{\frac{2\sqrt{1-\tilde{\alpha}}}{g\tilde{\alpha}}}$ and $\lambda=\frac{2+g\sqrt{1-\tilde{\alpha}}}{-g\sqrt{1-\tilde{\alpha}}}$. The derivation of the analytical form can be found in the Supplemental Material Sup .

The zero-energy LL is topologically guaranteed once the topological charge is nonzero. Therefore it was regarded as legitimate that the zeroth LL gaps for two WNs of zero net chiarlity under a strong magnetic field, as what we see in Model A, Fig. 1 (a). The interpretation has to be corrected when the persistent zero-energy LLs is demonstrated in Model B which retains zero net chirality as Model A. Therefore, chirality will not be responsible for the zero-energy LLs. Still, the zero-energy modes should be dictated by topology. One can understand that the chirality is a high dimensional topological invariant and hence is not suitable for explaining the LL system, since the systems for $k_{z}=0$ is restricted to a 2D system perpendicular to the magnetic field. In a 2D manifold pierced by holes (Weyl nodes), a 1D topological number is the end. Because of the presence of chiral symmetry, the systems belong to the AIII symmetry class and are classified by the winding number, a $\mathbb{Z}$-type topological invariant Chiu et al. (2016). In the paper, we dub it vorticity. In the chiral basis, the phase $\phi$ in the off-diagonal entry of the Hamiltonian characterizes the vorticity defined to be $\nu=\frac{1}{2\pi}\oint_{S^{1}}d\mathbf{k}_{\parallel}\cdot\nabla\phi$, where $S^{1}$ is a loop enclosing a (or multiple) WN(s) on the $k_{z}=0$ plane. Referring to Eqs. (4) and (5), two WNs in Model A take opposite vorticities, but equal vorticity in Model B. (The sign of vorticity might change by changing the basis, but the relative sign between two is invariant.) We illustrate chiralities and vorticities of WNs for the two models as the conclusion in Figs. 1 (c) and (d). Therefore, the net vorticity in Model A is $\nu_{\mathrm{A}}=0$ and is $\nu_{\mathrm{B}}=2$ in Model B. According to the index theorem Ansourian (1977); Niemi and Semenoff (1983); Jackiw (1984); Niemi and Semenoff (1984) and generalizing it, the absolute value of vorticity can be witnessed by the number of zero-energy LLs in a magnetic field. We have also examined a tilted field in the $y$-$z$ plane to strengthen this proof Sup .

Dispersion along $k_{z}$.— The chiral anomaly is a phenomenon that in parallel magnetic and electric fields electric charges are transmitted from one WN to the other. To investigate this effect, we consider two pairs of WNs separated in $k_{z}$ with each pair as studied before. We modify our models by $k_{z}\rightarrow\frac{1}{2k_{\mathrm{V}}}(k_{z}^{2}-k_{\mathrm{V}}^{2})$, and dub the modified models as Model $\tilde{\mathrm{A}}$ and $\tilde{\mathrm{B}}$ respectively; in this way the separation is $2k_{\mathrm{V}}$ and the absolute velocity component in $z$ is still $v_{z}$.

We show the calculation results for Landau bands along $k_{z}$ in Fig. 2. Since the magnetic field is along $z$, two pairs of WNs, at different $k_{z}$ connected in dispersion, are independent. In Model $\tilde{\mathrm{A}}$, as the zeroth LLs at both $k_{z}=\pm k_{\mathrm{V}}$ are gapped in large magnetic fields, a 3D insulating phase is present, as shown in Fig. 2(a). In Model $\tilde{\mathrm{B}}$, by contrast, the protected zero-energy states extend into the chiral Landau bands and result in a 3D metallic phase, Fig. 2(b). Moreover, the fact that the two chiral Landau bands crossing at either $k_{z}=k_{\mathrm{V}}$ or $-k_{\mathrm{V}}$ have opposite slopes is the proof of two WNs taking opposite chiralities. This feature reveals that to characterize the zeroth Landau bands in a WSM unequivocally, two topological invariants, chirality and vorticity, are necessary.

Impurity scattering.— We then discuss other phenomena that may distinguish the two models which are characterized by the same chiralities but different vorticities. In weak magnetic fields, the conductivity is highly influenced by disorders. Contrary to normal materials, the topological semimetals undergo weak anti-localization in the absence of magnetic field due to the $\pi$-Berry phase from the WN to suppress backscattering Lu et al. (2015). The anti-localization phenomenon will fade away when inter-valley scattering is taken into consideration for the lack of topological protection. In chiral anomaly, the latter determines the scale of transport time. Therefore, inter-valley scattering will influence transport properties the most.

We emphasize that two models have inherent distinction in inter-valley scatterings. Set $|v_{x}|=|v_{y}|=|v_{z}|\equiv v$ in Eqs. (4) and (5) for simplicity. Since each WN looks similar itself, the intra-valley scattering makes no difference between the two models. But for inter-valley scattering, whether the Fermi velocity changes sign or not from one valley to the other will affect the scattering probability. We denote the inter-valley scattering potential by $U^{+,-}_{q,q^{\prime}}$ for a scattering from $q$ to $q^{\prime}$ (momentum relative to WNs). Under Born approximation, the average scattering rate is given by

[TABLE]

where $\xi_{q}=\hbar v_{F}q$. We realize that when the impurity is anisotropic as $p$-wave, the differences in two models will be identified. Take $p_{y}$-wave impurity for instance that the scattering potential $U^{+-}_{k,k^{\prime}}\sim(q_{y}-q^{\prime}_{y})$ changes sign in $y$. As the Fermi velocity $v_{y}$ have opposite signs at two valleys in Model B, which indicates a $\pi$-phase difference between electrons at two valleys, inter-valley scattering will be enhanced by a $p_{y}$-wave impurity. In contrast, inter-valley scattering is weaker in Model A owing to equal sign of $v_{y}$. We conclude the results in Table 1.

Surface states.— By solving WSM slabs with semi-infinity in the $z$ direction and a hard wall potential for $z>0$, we can have the surface states and the corresponding energies as a function of $(k_{x},k_{y})$. For Model A, the energy $E=-v_{\parallel}k_{y}$, and for Model B the energy is $E=\frac{v_{\parallel}}{k_{\mathrm{W}}}|k_{x}|k_{y}$. Therefore, the energy contours can be shown in Fig. 3, from which we found the Fermi arcs do not connect to each other in Model B, since the solved surface state wavefunction is not continous in $k_{x}=0$. (See details in Ref. Sup ).

Weyl-node annihilation.— The two models differ in mirror parities of the two bands, so they give different results after the pair WNs collide and annihilate [by tuning $k_{\mathrm{W}}^{2}$ in Eqs. (2) and (3) to negative values]. After collision, WNs in Model A will gap the system while they evolve into a gapless nodal ring on the mirror plane in Model B Sup . These are simply consequences of symmetry-guaranteed anti-band crossing and band crossing. However, we point out that these are also consistent with the topological conditions. Since the annihilation process does not break chiral symmetry and hence vorticity is conserved all the way. For $k_{\mathrm{W}}^{2}<0$, the gapped phase in Model A assures $\nu_{\mathrm{A}}=0$, and the nodal ring in Model B piercing through the plane perpendicular to the mirror plane accounts for $\nu_{\mathrm{B}}=2$. We remark that the annihilation of WNs into a nodal ring or not by collisions is consistent with the newly discovered conversion rule Sun et al. (2018). Based on this idea, we can predict that when the mirror symmetry is broken by perturbation and if chiral symmetry is preserved, two WNs in Model A can still annihilate into an insulating phase, while WNs in Model B can collide but not annihilate.

Summary.— We have clarified that the chiral LLs need the topological protection from vorticity instead of chirality. With a strong magnetic field, two neighbor WNs coalesce into Landau orbits that possess finite- or zero-energy zeroth LLs depending on the net vorticity being 0 or 2 in the plane perpendicular to the field. The persistence of zero energies with vorticity 2 is robust and has conceptual appeal in searching topologically protected states. This finding reveals that to characterize topology of a WN chirality as well as vorticity are required in certain cases. Moreover, we also demonstrate that the net vorticity of a pair of WNs affects impurity scattering, surface states, and the way of WN annihilation.

I Algorithm for finding Landau spectrum

To solve the Hamiltonian with variable $q$ and its derivative $\partial/\partial q$, we can use the language of raising and lowering operators. The replacement is

[TABLE]

which guarantees the $[a,a^{\dagger}]=1$. By imposing $a|0\rangle=0$, all the necessary relations are then found as $a|n\rangle=\sqrt{n}|n-1\rangle$, $a^{\dagger}|n\rangle=\sqrt{n+1}|n+1\rangle$ and $a^{\dagger}a|n\rangle=n\rangle$, where $n=0,\cdots,L,\cdots$ labels the basis with well defined particle numbers. The eigen-differential problem is then converted into a matrix problem, and we can numerically solve the Landau spectrum of the system by matrix diagonalization. Since numerically we always solve it with a finite matrix of size $L$, the truncation of the operators $a$ and $a^{\dagger}$ always break canonical commutation relations from the highest few levels. For operators of order $k$, such as $a^{\dagger}a^{\dagger}a^{\dagger}$ of order 3, the levels which break the commutation relations happen at the highest $k$ levels. Especially, when written in the basis of $|0\rangle,\dots,|L-1\rangle$, the form of $a^{\dagger}$ is of the $L\times L$ matrix as

[TABLE]

The highest $L$ basis breaks the $[a,a^{\dagger}]=1$ and produces pseudo zero eigenvalues from the state of $[0,0,\dots,0,1]^{T}$. Therefore, no matter how large the matrix is used, there always would be pseudo zero or near zero states from the highest few particle number basis. This causes serious problems since what we concern is the low level states near zero energies. Some algorithm may use the regularization by adding some big numbers at the highest fewest levels to reduce their contribution. However, this is inconvenient for our case since we do not know how many zero energy states exist a priori. For two or multiple zero eigen energy solutions, linear combination from these eigenstates is always allowed and we do not have good rules to rule out pseudo solutions without ruining the true solutions.

Since what we concern is only the low-lying Landau spectrum near the Fermi level 0, we know their contributions all come from the low particle number bases. We then develop an efficient algorithm to exclude the pseudo solutions. For operators of order $k$, the pseudo states come from the highest $k$ basis, and we can restrict the solutions to be in the basis of $|0\rangle,\dots,|L-k\rangle$ of the Hilbert space. This can be effectively achieved by truncating the operator of $L\times L$ matrix into matrix of $L\times(L-k)$. The full Hamiltonian of size of $2L\times 2L$ then becomes matrix of size of $2L\times 2(L-k)$. Using the singular value decomposition (SVD) factorization, we can find the eigen-energies of system without contamination from high lying states. In our case of Hamiltonian which is at most of order 3, we always drop the last 4 bases, namely choosing $k=4$.

For the $m\times n$ matrix $M$, there exists the SVD factorization to be of the form of $M=U\Sigma V^{\dagger}$, where $U$ is an $m\times m$ unitary matrix whose columns are called the left-singular vectors of $M$, $V$ is an $n\times n$ unitary matrix with columns called right-singular vectors of $M$, and $\Sigma$ is a diagonal $m\times n$ matrix with non-negative real numbers on the diagonal. The right-singular vectors of $M$ are a set of orthonormal eigenvectors of $M^{\dagger}M$. For our purpose, the right-singular matrix $V$ serves to find the eigenstates and thus determines the eigen-energies. Below we will demonstrate how to find the low energy spectrum of interest.

Suppose the low energy eigenstates for the Hamiltonian $H$ has the support at most up to $L_{k}$, namely the mixture components from Landau levels than $L_{k}$ are zero. In the following SVD approach to get rid of pseudo solutions, we must guarantee that $L-k>L_{k}$. This can be always be achieved since $L_{k}$ is usually not very large and we can choose large enough $L$ to guarantee this requirement. The number of truncated columns $k$ can be chosen to be small, say $k=4$ for the system of order 3. Then we can write down the eigenstate of interests to be in the form of

[TABLE]

where the $(0,\dots,0)^{T}$ are located at the last $k$ Landau levels to be truncated at the up and down spin space separately. Both $\tilde{\phi}$ and $\tilde{\chi}$ are columns of size $(L-k)\times 1$ and have support up to $L_{k}$. Therefore, the weighting components of the $\{L_{k}+1,\cdots,L-k\}$ levels for them are actually zero. Assume that $k=2$ in the following sketch of proof, and then the finite size Hamiltonian $H$ of $2L\times 2L$ matrix can be written as

[TABLE]

where the $2k$ columns filled with $\times$ denotes the columns to be truncated. For the eigenstate $\psi$, we can have $H\psi=E\psi$, in which we are interested in low energy $E$ regime. Since components of the last $\{L-k+1,\dots,L\}$ levels for the eigenstates of interest are zero, The columns with $\times$ for the Hamiltonian actually have no effect. We can then drop them and collect the truncated Hamiltonian denoted as $\tilde{H}$ to be

[TABLE]

which is $2L\times 2(L-k)$ matrix. We then do the SVD factorization to have $\tilde{H}=U\Sigma V^{\dagger}$, with each column vector of $V$ denoted as

[TABLE]

which is a column of size $2(L-k)\times 1$. Padding with zeros in the form of Eq. (8) for $\tilde{\psi}$ to become $\psi$, we can have $\psi$ as the eigenstate of $H^{2}$ with eigenvalue of $E^{2}$. Doing some linear combination of eigenstates with the same eigenvalue $\lambda=E^{2}$, we obtain the eigenstate $\psi^{\prime}$ satisfying the eigenequation $H\psi^{\prime}=E\psi^{\prime}$. If such eigenstate cannot be found, it means that the eigenstate with eigenvalue $E$ for the system has components mixing from Landau levels no smaller than $L$ such that the form of Eq. (8) with the chosen size $L$ cannot the eigensolutions of $H$. In such case, we increase the numerical system size $L$ until the eigensolutions can be found for the low energy regime.

The reason that the $\psi$ constructed from $\tilde{\psi}$ from the SVD can be eigenstate of $H^{2}$ is simple. The Hermitian conjugate of the Hamiltonian $H$ is

[TABLE]

such that

[TABLE]

where

[TABLE]

The $\tilde{\psi}$ obtained from SVD of $\tilde{H}$ is just the eigenstates of $\tilde{H}^{\dagger}\tilde{H}$. As long as the eigenstate $\psi$ satisfying $H\psi=E\psi$ has the support $L_{k}<L-k$, namely components from highest $k$ levels are zero, the form of Eq. (8) constructed from $\tilde{\psi}$ would fall into eigen-solutions of $H^{2}\psi=H^{\dagger}H\psi=E^{2}\psi$.

II The analytical solutions for Model B with $\tilde{\alpha}=0$

The spectrum of Model B with $\tilde{\alpha}=0$ for small $g$ can also be found. These solutions are also plotted as red dots in Fig. 1 (b) in the main text for comparison. Model B under field defined in the main text and is rewritten here

[TABLE]

where the prime stands for a system under a magnetic field.

In the following we present the derivations for the solutions when $\tilde{\alpha}=0$ and $k_{z}=0$ for Model B. The Hamiltonian to solve is $H^{\prime}_{\mathrm{B}}$ in Eq. (9) with $\tilde{\alpha}=0$ and $q_{z}=0$, and we assume the eigenstate to be $(\chi(q),\Phi(q))^{T}$. By squaring $H_{\mathrm{B}}^{\prime}$, we can decouple $\chi$ and $\Phi$ as

[TABLE]

where the energy is rescaled to a dimensionless quantity defined as $\varepsilon=E/\frac{\omega_{c}^{2}}{4E_{\mathrm{VH}}}$. Owing to

[TABLE]

we firstly take $\chi(q)=\frac{1}{\sqrt{q}}e^{-\frac{q^{2}}{2g}}\tilde{\chi}(q)$ and $\Phi(q)=\frac{1}{\sqrt{q}}e^{-\frac{q^{2}}{2g}}\tilde{\Phi}(q)$ into the differential equations and obtain

[TABLE]

where $\lambda^{2}=\frac{1}{g^{2}}-\varepsilon^{2}$. Note that the exponential factor $e^{-\frac{q^{2}}{2g}}$ is to eliminate the quartic term $q^{4}$. Assuming $\lambda>0$, we continue to take $\tilde{\chi}(q)=q^{\lambda}f_{1}(q)$ and $\tilde{\Phi}(q)=q^{\lambda}f_{2}(q)$ and have

[TABLE]

Then we rescale the length by $\bar{q}=\frac{q}{\sqrt{g}}$, and obtain

[TABLE]

The final step is to change the variable from $\bar{q}$ to $\rho=\bar{q}^{2}$, which results in

[TABLE]

In the above, we already used

[TABLE]

Reformulating the last two equations, we have the so called associated (generalized) Laguerre equations

[TABLE]

where $n=\frac{1}{2g}-\frac{\lambda}{2}=0,1,2,\cdots$. The solutions $f_{1}$ and $f_{2}$ will be the associated Laguerre polynomials: $f_{1}(\rho)=L^{\lambda}_{n-1}(\rho)$ and $f_{2}(\rho)=L^{\lambda}_{n}(\rho)$ with $n=0,1,2,\cdots$. For $n=0$, we have to take $f_{1}=0$. As $n=\frac{1}{2g}-\frac{1}{2}\sqrt{\frac{1}{g^{2}}-\varepsilon^{2}}\geq 0$ is a non-negative integer, the energy eigenvalues are

[TABLE]

Here $E_{n}\leq\frac{1}{g}$ as $\lambda^{2}=\frac{1}{g^{2}}-\varepsilon^{2}>0$. The relevant state is the lowest Landau level $n=0$ which gives the zero-energy state $E_{n=0}=0$. It’s eigenstate are thus $(0,\Phi)^{T}$ with

[TABLE]

where $L^{\lambda}_{n}$ is a polynomial function of $q^{2}$ to degree $n$. Therefore, the normalizability demands $\lambda\geq\frac{1}{2}$, that is, $0\leq n\leq\frac{1}{2}\left(\frac{1}{g}-\frac{1}{2}\right)$. The missing states for large $n$ not satisfying the constraint become extended states whose spectra are continuous. This is an artefact of this model with $\alpha=0$ which omit the $k_{y}$ and $k_{z}$ dependence in the $\sigma_{x}$ term, since it produces an open equienergy contour for $E\geq E_{\mathrm{VH}}$ and hence the cyclotron orbit is not confined.

III The analytical form of zero energy solutions for Model $H_{\mathrm{B}}$

Here we demonstrate the analytical form of the zero energy solutions for Model B with $k_{z}=0$. The zero-eigenvalue problem is $H_{\mathrm{B}}^{\prime}\Psi=0$, where $H_{\mathrm{B}}^{\prime}$ is an $2\times 2$ off-diagonal matrix with elements

[TABLE]

The prime means the system in a magnetic field. We could try solutions either as $\Psi=\left(0,\psi\right)^{T}$ or $\Psi=\left(\psi,0\right)^{T}$, but it turns out that the second choice is not normalizable. Then we are going to solve the differential equation

[TABLE]

Its large-$q$ limit can be conquered by setting $\psi(q)=e^{-\kappa q^{2}}\Phi(q)$ and taking it into Eq. (21) to obtain $\kappa$. As $\kappa>0$ for normalizability, we have $\kappa=\frac{1+\sqrt{1-\tilde{\alpha}}}{g\tilde{\alpha}}$. The other choice is $\kappa=\frac{1-\sqrt{1-\tilde{\alpha}}}{g\tilde{\alpha}}$, but this would lead to additional exponential term from modified Hermite differential equation later. After absorbing the additional exponential term, this results in the same result from $\kappa=\frac{1+\sqrt{1-\tilde{\alpha}}}{g\tilde{\alpha}}$ and therefore we focus on the plus sign choice. It follows the differential equation for $\Phi(q)$, which is

[TABLE]

By defining $\xi=\sqrt{\frac{2\sqrt{1-\tilde{\alpha}}}{g\tilde{\alpha}}}$, $\lambda=\frac{2+g\sqrt{1-\tilde{\alpha}}}{-g\sqrt{1-\tilde{\alpha}}}$ and $\theta=\xi q$, the equation is then transformed into the Hermite differential equation

[TABLE]

The solution for $\Phi$ is ${{}_{1}F_{1}}(-\frac{1}{4}\lambda;\frac{1}{2};\theta^{2})$ and $H_{\lambda/2}(\theta))$, where ${}_{1}F_{1}(a;b;x)$ and $H_{\mu}(x)$ are the Kummer confluent hypergeometric function and the Hermite polynomial respectively.

It is known that $H_{\lambda/2}(\xi q)$ is not purely even or odd with respect to $q$. As one can find that the even part of $H_{\lambda/2}(\xi q)$ is actually ${}_{1}F_{1}(-\frac{1}{4}\lambda;\frac{1}{2};\xi^{2}q^{2})$, it is consistent that after reconstruction the zero-energy eigenfunctions are either of even or odd parity in $q$.

The zero-energy eigen-functions are plotted in Fig. 4. The $g$ value proportional to the field strength controls the coupling between the two nodes. In small $g$, the Landau orbits are well separated and have each center around the Weyl nodes. The orbits start to overlap for larger $g$ so that they tend to move toward the mirror plane. The trend of tuning $g$ for both $q$-even and $q$-odd solutions are the same.

IV Models including pairs of Weyl nodes separated in $k_{z}$

To study chiral anomaly of pairs of close Weyl nodes (here separated in $k_{x}$ direction) under the effect of magnetic field along $z$ direction, we must include the other pair of nodes separated in $z$ direction. By shifting the Weyl nodes in $H_{\mathrm{A}}$ and $H_{\mathrm{B}}$ to $k_{z}=k_{\mathrm{V}}$ and include the other pair of Weyl nodes at $k_{z}=-k_{\mathrm{V}}$, we can have the modified models, denoted as $H_{\tilde{\mathrm{A}}}$ and $H_{\tilde{\mathrm{B}}}$, in the form of

[TABLE]

The Weyl nodes of the first pair are then located at $(\pm k_{\mathrm{W}},0,k_{\mathrm{V}})$, while the other pair is located at $(\pm k_{\mathrm{W}},0,-k_{\mathrm{V}})$. Since we are looking at physics near the Weyl nodes and their low energy spectrum, the $\alpha k_{\parallel}^{2}$ term does not play much role. In most of the time they can be even dropped. The value of $\alpha$ under discussion is therefore small, and the $\alpha k_{\parallel}^{2}$ term affects mainly the dispersion in higher energy and does not influence the low energy of interest much.

Here Model $\tilde{A}$ and Model $\tilde{B}$ still have the mirror plane $k_{x}=0$, and we do not put in additional symmetry relation between the first pair of WNs and the second pair for simplicity. In this way, the effect of different choices of $M_{x}=\sigma_{0}$ or $\sigma_{x}$ can be clearly seen. Besides, usually in real materials, additional WNs are far away from the pair of interest such that their effect can be discarded since they are far from reach of the magnetic length scale under reasonable field strength. Therefore, for simplicity we only compare one mirror plane with different operator choices in order to elucidate the symmetry impacts.

In the usual Weyl semimetal, the separation of WNs in the $k_{z}$ direction is larger than the $k_{\mathrm{W}}$, i.e. $k_{\mathrm{V}}>k_{\mathrm{W}}$. Due to the other pair of nodes, the term of $q_{z}\sigma_{z}$ in Model A Hamiltonian under magnetic field $H_{\mathrm{A}}^{\prime}$ is replaced by $\frac{1}{g}A_{z}(q_{z}^{2}-1)$ in $H_{\tilde{\mathrm{A}}}^{\prime}$ under field while the term $q_{z}q\sigma_{z}$ in $H_{\mathrm{B}}^{\prime}$ is replaced by $\frac{1}{g}A_{z}q(q_{z}^{2}-1)$ in $H_{\tilde{\mathrm{B}}}^{\prime}$, where $A_{z}=g\frac{v_{\parallel}k_{\mathrm{V}}}{2}\left(\frac{\omega_{c}^{2}}{4E_{\mathrm{VH}}}\right)^{-1}=(\frac{v_{\parallel}}{v_{x}})(\frac{k_{\mathrm{V}}}{k_{\mathrm{W}}})$ with the new definition of the dimensionless $q_{z}=k_{z}/k_{\mathrm{V}}$. The $-\frac{\alpha}{2m}k_{\parallel}^{2}$ term in both models $H_{\mathrm{A}}^{\prime}$ and $H_{\mathrm{B}}^{\prime}$ would correspondingly become $\frac{\omega_{c}^{2}}{4E_{\mathrm{VH}}}\left[g\frac{\tilde{\alpha}_{y}}{4}\frac{\partial^{2}}{\partial q^{2}}-\frac{1}{g}\tilde{\alpha}_{z}q^{2}_{z}\right]$, where $\tilde{\alpha}_{y}=\alpha(\frac{v_{x}}{v_{\parallel}})^{2}$, $\tilde{\alpha}_{z}=\alpha(\frac{k_{\mathrm{V}}}{k_{\mathrm{W}}})^{2}$. Since $q$ and $q_{z}$ are independent to each other, $q_{z}$ as a good quantum number can be treated as a parameter and the Hamiltonian is solved at fixed $q_{z}$ each time.

Explicitly, the Hamiltonian under field to solve for Model $\tilde{A}$ is then

[TABLE]

while the Hamiltonian for Model $\tilde{B}$ is then

[TABLE]

V Solutions for rotated magnetic field in the yz plane

Here we consider the magnetic field rotated in the yz plane which is still perpendicular to the two nodes separated in the $k_{x}$ direction. Since we mainly concern about whether the chiral Landau levels and the zero energy levels can maintain, we focus on the case following discussions of model $H_{\mathrm{B}}$ and $H_{\tilde{\mathrm{B}}}$. The coefficients for $k_{y}$ and $k_{z}$ in the Hamiltonian can be different but the physics is the same. For general purpose, we can always rescale $k_{y}$ and $k_{z}$ such that $\alpha k_{\parallel}^{2}$ have the same coefficients for $k_{y}$ and $k_{z}$ while linear terms $\sim v_{y}k_{x}k_{y}\sigma_{y}$ and $\sim v_{z}k_{x}k_{z}\sigma_{z}$ have different parallel velocities $v_{y}$ and $v_{z}$ for $k_{y}$ and $k_{z}$ respectively. Note that the $k_{x}=0$ plane still need to be dispersional for all $(k_{y},k_{z})$ to ensure close energy contours such that the magnetic orbits can be formed under all rotated field directions. In the Hamiltonian, this means that $\alpha$ can be small but cannot be zero. For single pair of Weyl nodes, the Hamiltonian is

[TABLE]

which is already defined in the main text. The coefficient $\alpha$ mainly influence dispersion in higher energies. Usually $\alpha$ is small, and therefore does not change the low energy spectrum much, including the zero energy levels of concern. The simpler way to deal with rotated field is to define new momentum coordinates $k^{\prime}_{y}$ and $k^{\prime}_{z}$ where the new $\hat{z}$ direction is along the field. Suppose the field $\mathbf{B}=B(\sin\theta\hat{y}+\cos\theta\hat{z})$, then the new momentum coordinates are defined as $k^{\prime}_{y}=\cos\theta k_{y}-\sin\theta k_{z}$ and $k^{\prime}_{z}=\sin\theta k_{y}+\cos\theta k_{z}$, and still $k_{\parallel}^{2}=k_{y}^{2}+k_{z}^{2}={k^{\prime}_{y}}^{2}+{k^{\prime}_{z}}^{2}$. In such definition, we take the advantage of $k^{\prime}_{z}$ still being a good quantum number and $k_{x}k^{\prime}_{y}\rightarrow l_{B}^{-2}(k_{x}\bar{x}+\frac{i}{2})$ similar to model $H_{\mathrm{B}}$ with modified guiding center $x_{0}=l^{2}_{B}k^{\prime}_{y}$. With definition of $\sigma^{\prime}_{y}=\cos\theta\sigma_{y}-\sin\theta\sigma_{z}$ and $\sigma_{z}^{\prime}=\sin\theta\sigma_{z}+\cos\theta\sigma_{z}$, the Hamiltonian in new coordinates can be written as

[TABLE]

The dimensionless momentum in the field direction is defined as $q^{\prime}_{z}=v_{z}k^{\prime}_{z}\left(\frac{\omega_{c}^{2}}{4E_{\mathrm{VH}}}\right)^{-1}=(\frac{2}{g})(\frac{v_{z}}{v_{x}})(\frac{k^{\prime}_{z}}{k_{\mathrm{W}}})$. The dimensionless quantities now are $\tilde{\alpha}_{y}=\alpha(\frac{v_{x}}{v_{y}})^{2}$ and $\tilde{\alpha}_{z}=\alpha(\frac{v_{x}}{v_{z}})^{2}$. Therefore, the Hamiltonian under field in rotated coordinate is then

[TABLE]

where the prime on the left side refers to the Hamiltonian under magnetic field. The independent parameters are $\tilde{\alpha}$, $g$, and rotated angle $\theta$. Among them, $\tilde{\alpha}$ and $g$ are related to materials properties, i.e. Weyl nodes quantities, and $g$ and $\theta$ are related to field amplitude and direction respectively. As an example, we present the case of $\tilde{\alpha}=0.05$ and with fixed value of $g=0.8$, we rotate the angle $\theta$ of field with respect to $\hat{z}$ axis from $0^{\circ}$ to $180^{\circ}$ and present the result in Fig. 5 (a). It can be found that the zero energies persist in all angles $\theta$ in the plane of $k^{\prime}_{z}=0$ perpendicular to the rotated field.

To see if the chiral anomaly can remain for rotated field when two pairs of Weyl nodes are located at $k_{z}=\pm k_{\mathrm{V}}$, the Hamiltonian for Model A and Model B are Eq. (24), and we focus on discussing Model B.

Similarly, when written in the new coordinate $(k_{x},k^{\prime}_{y},k^{\prime}_{z})$ of the rotated frame, Model B Hamiltonian is then

[TABLE]

Different to $H_{\mathrm{B}}$, the dimensionless momentum along the field direction is defined as $q_{z}=k^{\prime}_{z}/k_{\mathrm{V}}$. Therefore, the dimensionless parameters involving $k^{\prime}_{z}$ would change the dependence from $v_{x}/v_{\parallel}$ to $k_{\mathrm{V}}/k_{\mathrm{W}}$. Explicitly, the dimensionless parameters are defined as follows. $\tilde{\alpha}_{y}=\alpha(\frac{v_{x}}{v_{y}})^{2}$, $\tilde{\alpha}_{z}=\alpha(\frac{k_{\mathrm{V}}}{k_{\mathrm{W}}})^{2}$, $A_{z}=\left(\frac{v_{z}}{v_{x}}\right)(\frac{k_{\mathrm{V}}}{k_{\mathrm{W}}})$, $A_{y}=(\frac{v_{x}v_{z}}{4v_{y}^{2}})(\frac{k_{\mathrm{W}}}{k_{\mathrm{V}}})$ and $A_{yz}=\frac{v_{z}}{2v_{y}}$, where the definition of $A_{z}$ is the same as that in single pair of Wely nodes. Not all the defined dimensionless parameters are independent. Among them, the independent parameters are chosen to be $\alpha$, $k_{\mathrm{W}}/k_{\mathrm{V}}$, and $\frac{v_{x}}{v_{y}}$. The Hamiltonian under magnetic field to solve is then

[TABLE]

With the reasonable choice of $\alpha=0.05$, $v_{y}/v_{x}=v_{z}/v_{x}=0.5$, and $k_{\mathrm{W}}/k_{\mathrm{V}}=0.2$, we present the case of $g=0.8$ and $\theta=30^{\circ}$ in Fig. 5 (b). The corresponding values are $\tilde{\alpha}_{y}=0.2$, $\tilde{\alpha}_{z}=1.25$, $A_{y}=0.1$, $A_{z}=2.5$, and $A_{yz}=0.5$. It is found that the chiral anomaly still remains since the chiral Landau levels near each pair of nodes are robust.

VI surface states

We are going to solve Weyl semimetal slabs for Model A and Model B. In the $x$ and $y$ directions, the sizes are infinity, while it is semi-infinity in the $z$ direction. Assume that the Weyl semimetal systems are built for $z<0$ adjacent to vacuum for $z>0$. We will analyze Model B first and then Model A since the former model is new to us.

VI.1 Model B

To model a vacuum-semimetal interface, we introduce a mass term in the Hamiltonian as

[TABLE]

where $M(z)=0$ for $z<0$ and $M(z)=M\rightarrow\infty$ for $z>0$. Here we simplify the model by dropping some constants which will be restored later.

For the localized surface states, we take the ansatz:

[TABLE]

where $\mathrm{Re}\lambda_{<}$ and $\mathrm{Re}\lambda_{>}$ are positive. For $z>0$, in the limit of $M\rightarrow\infty$, we have, by taking the ansatz into $H_{\mathrm{B}}\psi=0$ and neglecting small numbers,

[TABLE]

[TABLE]

leading to

[TABLE]

Since $\lambda_{>}>0$, we have $u=v$ for $k_{x}>0$ and $u=-v$ for $k_{x}<0$. As a result, we have boundaries conditions as

[TABLE]

With these boundary conditions, we take the ansatz for $z<0$ into $H_{\mathrm{B}}\psi=E\psi$ and we have, for $k_{x}>0$,

[TABLE]

Equating the real parts and the imaginary parts separately, we have $E=k_{x}k_{y}$ and

[TABLE]

In order to have $\mathrm{Re}\lambda_{<}>0$, $k_{x}$ is limited by $k_{x}<\sqrt{k_{\mathrm{W}}^{2}+\alpha k_{y}^{2}}$. Similarly, for $k_{x}<0$, we have $E=-k_{x}k_{y}$ when $-\sqrt{k_{\mathrm{W}}^{2}+\alpha k_{y}^{2}}<k_{x}$. In conclusion, when putting back omitted constants, the surface states (the Fermi arcs) survive in $|k_{x}|<\sqrt{k_{\mathrm{W}}^{2}+\alpha k_{y}^{2}}$ and take energy

[TABLE]

The corresponding wave functions are

[TABLE]

where

[TABLE]

VI.2 Model A

With the same trick one can also solve the surface states for Model A. We skip the deviations and only show the result below: for $|k_{x}|<\sqrt{k_{\mathrm{W}}^{2}+\alpha k_{y}^{2}}$

[TABLE]

where

[TABLE]

with energy $E=-v_{\parallel}k_{y}$.

We point out main difference between the two models. Model A shows typical understanding of a Fermi arc connecting two Weyl nodes with linear dispersion $E\propto k_{y}$. However, in Model B the surface-state wave function is not continuous at $k_{x}=0$, which is proportional to $\left(1,i\right)^{T}$ on one side and $\left(1,-i\right)^{T}$ on the other. So it indicates that there are two Fermi arcs that do not connect these Weyl nodes (but to other pairs), presenting a hyperbola dispersion $E\propto|k_{x}|k_{y}$.

VII Annihilation of Weyl nodes

The pair of Weyl nodes move toward each other if we tune the parameter $k_{\mathrm{W}}^{2}$ smaller and eventually will collide with each other and annihilate. Before collision, the two Weyl nodes remain intact and the system is still gapless for both Model A and Model B. However, the outcomes are different for the two models after the collision, i.e. $k_{\mathrm{W}}^{2}<0$. For Model A, the system will be gapped out with minimum energy formed by a ring. On the contrary, model B still remain gapless, while the two nodes annihilate into a nodal ring with zero energies. To demonstrate this more clearly, we rescale the parameters to simplify the Hamiltonian for the two models as

[TABLE]

[TABLE]

Before collision where $k_{\mathrm{W}}>0$, the Weyl nodes are located at $(\pm k_{\mathrm{W}},0,0)$. By tuning $k_{\mathrm{W}}^{2}$ to become negative, Model A has gap $E_{g}=2|k_{\mathrm{W}}|$ at $k=(0,0,0)$ or gap of $E_{g}=\frac{\sqrt{4\alpha|k_{\mathrm{W}}^{2}|-1}}{\alpha}$ at ring of $\sqrt{k_{y}^{2}+k_{z}^{2}}=\sqrt{\frac{2\alpha k_{\mathrm{W}}^{2}-1}{2\alpha^{2}}}$ for $k_{x}=0$ and other positions are also gapped. Model A is fully gapped unless accidental cases like $|k_{\mathrm{W}}^{2}|=\frac{1}{4\alpha}$. However, Model B remains gapless in which the WNs annihilate into a nodal ring which has zero energies. The nodal ring has the radius $k=\sqrt{\frac{|k_{\mathrm{W}}^{2}|}{\alpha}}$ in the $k_{x}=0$ mirror plane.

VII.1 The introduction of mirror symmetry breaking effect

We can further see the effect of mirror symmetry breaking terms induced by some perturbations to the Weyl nodes before and after the collision. Such perturbation can be realized several way, such as applying magnetic field to the system. Depending on specific material system Hamiltonian, the spin operators can in different combination of Pauli matrices, determined by system symmetries. Since the mirror operators $\mathcal{M}_{x}$ for Model A and Model B are known, in conjunction with combined symmetry $\mathcal{C}_{2}\mathcal{T}$, the allowed forms of spin operators can be determined. Further restriction of allowed forms is possible if we have more symmetry constraints, but this definitely depends on details of the systems. The procedure to determine forms of spin operators that can be concordant with symmetry requirement is shown in the next section.

Here we demonstrate effects of some possible mirror symmetry breaking terms induced by applying magnetic field to the system. In Model A, the allowed spin operator $y$ component can be combination of $s_{y}=\{k_{x}\sigma_{0},k_{x}\sigma_{x},k_{x}\sigma_{y}\}$, which can break mirror symmetry if we apply magnetic field in $y$ direction. For simplicity, we restrict the discussion to the form of $k_{x}\sigma_{y}$. The mirror symmetry breaking perturbation is therefore $\Delta k_{x}\sigma_{y}$ where the $\Delta$ is small perturbation determined by the strength of magnetic field. As usual in the $k_{\mathrm{W}}^{2}>0$, the locations of WNs are determined by each component of Pauli matrices to be zeros. The Weyl nodes are still topologically protected to exist but shifted to position of $(\pm\sqrt{\frac{k_{\mathrm{W}}^{2}}{1-\alpha\Delta^{2}}},\mp\Delta\sqrt{\frac{k_{\mathrm{W}}^{2}}{1-\alpha\Delta^{2}}},0)$. As for the annihilation results by tuning $k_{\mathrm{W}}^{2}<0$ to make Weyl nodes to collide, the system is still fully gapped where the gap $E_{g}=2|k_{\mathrm{W}}|$ at the origin and $E_{g}=\frac{\sqrt{4\alpha|k_{\mathrm{W}}^{2}|-1}}{\alpha}$ at $k=(0,0,\pm\sqrt{\frac{2\alpha k_{\mathrm{W}}^{2}-1}{2\alpha^{2}}})$.

For Model B, the mirror symmetry breaking perturbation can be $\Delta\sigma_{y}$ induced by the magnetic field if the spin y component $s_{y}=\sigma_{y}$ fulfilling the symmetry requirement. The WNs are shifted to $k=(\pm\sqrt{\frac{k_{\mathrm{W}}^{2}+\sqrt{k_{\mathrm{W}}^{4}+4\alpha\Delta^{2}}}{2}},\mp\sqrt{\frac{-k_{\mathrm{W}}^{2}+\sqrt{k_{\mathrm{W}}^{4}+4\alpha\Delta^{2}}}{2\alpha}},0)$, where the relative sign of $k_{x}$ and $k_{y}$ are determined by the sign of $\Delta$. After WNs annihilation when tuning $k_{\mathrm{W}}^{2}<0$, the system remains gapless but no longer existing a nodal ring. The zero gapless positions are located at $k=(\pm\sqrt{\frac{-|k_{\mathrm{W}}^{2}|+\sqrt{k_{\mathrm{W}}^{4}+4\alpha\Delta^{2}}}{2}},\mp\sqrt{\frac{|k_{\mathrm{W}}^{2}|+\sqrt{k_{\mathrm{W}}^{4}+4\alpha\Delta^{2}}}{2\alpha}},0)$.

Here we demonstrate the differences between Model A and Model B in their behaviour of WNs annihilation when the parameter $k_{\mathrm{W}}^{2}$ are tuned from positive to negative. Although both of them describe the pair of WNs when $k_{\mathrm{W}}^{2}>0$, Model A generally will be gapped out when $k_{\mathrm{W}}^{2}<0$ while Model B remains gapless located at a nodal ring. Even if we apply the magnetic field to break mirror symmetry, the feature of gapful Model A and gapless Model B remains the same.

VIII Allowed forms of spin operators

First of all, we have to point out that the Pauli matrices in the Hamiltonians $H_{\mathrm{A}}$ and $H_{\mathrm{B}}$ stand for the pseudo-spin describing two bands’ degrees of freedom not real spin. With spin-orbit interaction, spin, orbital and momentum are strongly coupled in bands, so that there is no simple or universal relation between the pseudo-spin and spin. Here we will try to extract spin degrees of freedom from the pseudo-spin based on symmetry point of view and give readers an idea how spin is included in the pseudo-spin for our models. The relation will no doubt depend much on details of systems.

We start with Model B, in which the mirror operator chosen is $\mathcal{M}_{x}=\sigma_{x}$. Under the mirror reflection, momentum and spin change as $\mathbf{k}=(k_{x},~{}k_{y},~{}k_{z})\mapsto(-k_{x},~{}k_{y},~{}k_{z})$ and $\mathbf{s}=(s_{x},~{}s_{y},~{}s_{z})\mapsto(s_{x},~{}-s_{y},~{}-s_{z})$. At the same time, the pseudo-spin changes as $(\sigma_{x},~{}\sigma_{y},~{}\sigma_{z})\mapsto(\sigma_{x},~{}-\sigma_{y},~{}-\sigma_{z})$. We find that the pseudo-spin and spin have the same transformation and might conclude that they are identical. However, we cannot have such conclusion because we have not compare their transformations under all possible symmetry operations. As a results, with only the mirror symmetry, we can claim that the spin components might contain contributions as follows:

[TABLE]

where linear combinations of elements in the curly brackets are possible with proper normalization.

To reduce the complexity, we consider that there also exists the combined symmetry of twofold rotation about $z$ and time-reversal symmetry, denoted by $\mathcal{M}_{2}T$. Suppose that the $k_{z}$ is relative to $k_{z}=0$ or $\pi$, $\mathcal{C}_{2}\mathcal{T}$ makes $(k_{x},~{}k_{y},~{}k_{z})\mapsto(k_{x},~{}k_{y},~{}-k_{z})$, $(s_{x},~{}s_{y},~{}s_{z})\mapsto(s_{x},~{}s_{y},~{}-s_{z})$. Many antiunitary operators could be used for $\mathcal{C}_{2}\mathcal{T}$ with the restriction that

[TABLE]

either for spin-0 or spin-1/2 systems. However, when we refer to the transformation

[TABLE]

we find that it has to be $\mathcal{C}_{2}\mathcal{T}=\sigma_{x}K$, where $K$ is the complex conjugation operation. With $\mathcal{C}_{2}\mathcal{T}$, the spin will reduce its compositions as follows:

[TABLE]

With the spirit, we can obtain spin from the pseudo-spin for Model A too. Here $\mathcal{M}_{x}=\sigma_{0}$ and $\mathcal{C}_{2}\mathcal{T}=\sigma_{x}K$, we show them as

[TABLE]