Topologically distinct Weyl fermion pairs

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 tunneling does not always expel LLs from zero energy because the number of zero-energy modes is protected by 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 more investigation in the future.


Results
Models. We consider a pair of WNs with opposite chiralities sitting on two sides of the mirror plane M x .
Their separation 2k 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 present, in other regions of the BZ can be ignored as they are much farther away. To have we find that there are two possible choice of M x and thus two possible models. We dub them Model A and Model B: Model A is for M x = σ 0 that the two bands have equal mirror parity in the mirror plane, while Model B is for M x = σ x that two bands have opposite mirror parities. Specifically, the Hamiltonians are written as where k 2 � = k 2 y + k 2 z . In general, coefficients for k y and k z can be different, but the physics are the same. A k -linear term is allowed to appear in the σ x term, but it is omitted for an elegance reason. We have confirmed that its presence once being small does not change qualitative conclusions. The Planck constant ℏ is set as unity throughout the paper. We note that the origin is meaninglessly specified and might not be at a time-reversalinvariant momentum.
Both models contain two WNs located at (±k W , 0, 0) . Expanding around the WNs, they approximate, to linear order, as here χ labels the chirality and also position of the WN, and v x = k W /m and v y = v z = v � . The values of k W , m, v y and v z are all assumed to be positive without loss of generality. α ( > 0 ) in H B is required in order to have a saddle point at energy E VH = k 2 W /2m , above which close energy contours are assured. Without special regard to the symmetry or phase, Model A was usually adopted to study the effect of a pair of nodes 38,45 . Applying the magnetic field along the perpendicular z direction to this system, we solve the LL spectrum by substituting k in Eq. (2) by � = � k + e � A . We choose the Landau gauge � A = Bxŷ , so we make k x → x = k x , k z → z = k z , and k y → y = l −2 Bx , where the magnetic length l B = 1 eB and x is the coordinate relative to the guiding center x 0 = −l 2 B k y . x is conjugate to k x by quantizing x → i ∂ ∂k x . The WN separation is used as a measure to define the dimensionless momentum scale as q = k x /k 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 ω c = 2v x v y l −1 B and decrease with the energy barrier E VH , defined as g = . g is proportional to the magnetic field B, and the appreciable coupling g ≈ 1 is achieved when the magnetic length l B is comparable to the scale defined by k −1 W . [Take Weyl semimetal TaAs for example. The nodes separation of W 1 is k W = 0.0072 2π a , and the Fermi velocity in the conduction band is (v x , v y , v z ) = (2.5, 1.2, 0.2) × 10 5 m/s 46 . The coupling g = 1 under the field strength B ≈ 11.88 T having magnetic length l B = 7.44 nm, which is close to the length 7.60 nm defined by k −1 W .] By defining the remaining variables into dimensionless quantities, q z = v � k z α , we can study the Hamiltonian under magnetic field as a function of q and q z , which is written as We numerically solve this system with raising and lowering operators a † and a (where q, ∂ ∂q ∼ a ± a † ). Special treatment is developed to solve it more efficiently and the details are shown in the "Methods".
The LL spectrum at k z = 0 , i.e. q z = 0 , with respect to g is shown in Fig. 1a. In the limit of g → 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 ≈ 0.3 (decrease with levels) in Fig. 1a. We have ascribed the band splits to the χ v x q x σ x + v y q y σ y + v z q z σ z . www.nature.com/scientificreports/ magnetic tunneling in Ref. 38 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 {H, σ z } = 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. Then, 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 and k x k z → (� x � z + � z � x )/2 = k x k z . In terms of dimesionless parameters defined above, the Hamiltonian under magnetic field becomes 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. 1b. For small g, the effect of α is small, so we provide the results for α = 0 as shown by red dots in Fig. 1b. The analytical solutions are g − n with n ∈ Z ≥ 0 and each energy has twofold degeneracy. The LL energies emerge from the typical LL spectrum as E n = √ nω c in the limit of g → 0 and deform with finite g. However, the analytical solutions for α = 0 are only applicable in a low-energy region. At high energies, with large g or n, the LL Figure 1. The LL spectra at k z = 0 with respect to the coupling measure g with α = 0.05 for both Model A (a), and for Model B (b). In (b), Labels "e" and "o" denote even and odd parity of states, and the red dots are the analytical solutions, except zero-energy ones, when α = 0 . www.nature.com/scientificreports/ quantization makes no sense when E > E 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 ≤ 1 2 1 g − 1 2 . 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. 1b. 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 α = 0 . The band splitting is reasonable as seen in a symmetric double-well with finite tunneling 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 α . They can be directly verified by solving the zero-eigenvalue problem. The eigenfunctions are found to be 0, e −κq 2 �(q) T , where �(q) can be either 1 F 1 − 1 4 ; 1 2 ; ξ 2 q 2 or H /2 (ξ q) , the former being the Kummer confluent hypergeometric function and the latter the Hermite polynomial, indicating a double degeneracy. Here κ = 1+ The derivation of the analytical form can be found in the "Methods".
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. 1a. 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. For a 2D space with a defect (a Weyl node), a 1D topological number is required. Because of the presence of chiral symmetry, the systems belong to the AIII symmetry class and are classified by the winding number, a Z-type topological invariant 47 . In the paper, we dub it vorticity. In the chiral basis, the phase φ in the off-diagonal entry of the Hamiltonian characterizes the vorticity defined to be ν = 1 2π S 1 dk � · ∇φ , 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 Fig. 1c,d. Therefore, the net vorticity in Model A is ν A = 0 and is ν B = 2 in Model B. According to the index theorem [48][49][50][51] 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 in the "Methods".
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. It should be answered whether two models show distinct consequences. To investigate this phenomenon, we consider two pairs of WNs separated in k z with each pair as studied before. We modify our models by k z → 1 2k V k 2 z − k 2 V , and dub the modified models as Model Ã and B , respectively; in this way the separation is 2k 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 Ã , as the zeroth LLs at both k z = ±k V are gapped in large magnetic fields, a 3D insulating phase is present, as shown in Fig. 2a. In Model B , by contrast, the protected zero-energy states extend into the chiral Landau bands and result in a 3D metallic  www.nature.com/scientificreports/ phase, Fig. 2b. Moreover, the fact that the two chiral Landau bands crossing at either k z = k V or −k 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 π-Berry phase from the WN to suppress backscattering 52 . The antilocalization 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, intervalley 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 | ≡ 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 ′ for a scattering from q to q ′ (momentum relative to WNs). Under Born approximation, the average scattering rate is given by where ξ q = ℏv F q . We realize that when the impurity is anisotropic, for example p-wave, the differences in the two models is evident. Take a p y -wave impurity for instance, where the scattering potential U +− k,k ′ ∼ (q y − q ′ y ) changes sign in y. As the Fermi velocities v y have opposite signs at the two valleys in Model B, which indicates a π-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 summarize 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 � k y , and for Model B the energy is E = v � k 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 the "Methods").
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 2 W 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. These are simply consequences of symmetry-guaranteed anti-band crossing and band crossing. Table 1. Average inter-valley scattering rates 1 τ I for anisotropic impurities potentials. Here p-wave ( p x , p y , and p z ) impurities potentials are considered. The p x -impurity potential (second row) does not differentiate the two models, but the form of ∼ p � = p y or p z can tell the difference (last row). N F is the density of states at the Fermi energy and u I characterizes the strength of the impurity potential. Figure 3. www.nature.com/scientificreports/ However, we point out that these are also consistent with the topological conditions. Since the annihilation process does not break chiral symmetry, vorticity is conserved. For k 2 W < 0 , the gapped phase in Model A assures ν A = 0 , and the nodal ring in Model B piercing through the plane perpendicular to the mirror plane accounts for ν 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 44 . 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.

Impurity potential Model
In 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 for 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.
Closing remarks. Our theory is focus on local Weyl fermions on two sides of a mirror plane and it is not subject to the type of Weyl semimetals-time-reversal or inversion symmetry breaking. The theory adopts k · p two-band models for describing low-energy physics, which is sufficient for discussion of magnetic field effects since the magnetic field in labs is not strong enough to couple far-distant Weyl fermions. The Pauli matrices in the models are acting to the pseudo-spin in the band space; the possible relationship between spin and the pseido-spin can be found in the "Methods".

Methods
Algorithm for finding Landau spectrum. To solve the Hamiltonian with variable q and its derivative ∂/∂q , we can use the language of raising and lowering operators. The replacement is which guarantees the [a, a † ] = 1 . By imposing a|0� = 0 , all the necessary relations are then found as a|n� = √ n|n − 1� , a † |n� = √ n + 1|n + 1� and a † a|n� = n� , where n = 0, . . . , L, . . . 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 † always break canonical commutation relations from the highest few levels. For operators of order k, such as a † a † a † of order 3, the levels which break the commutation relations happen at the highest k levels. Especially, when written in the basis of |0�, . . . , |L − 1� , the form of a † is of the L × L matrix as The highest L basis breaks the [a, a † ] = 1 and produces pseudo zero eigenvalues from the state of [0, 0, . . . , 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�, . . . , |L − k� of the Hilbert space. This can be effectively achieved by truncating the operator of L × L matrix into matrix of L × (L − k) . The full Hamiltonian of size of 2L × 2L then becomes matrix of size of 2L × 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 × n matrix M, there exists the SVD factorization to be of the form of M = U V † , where U is an m × m unitary matrix whose columns are called the left-singular vectors of M, V is an n × n unitary matrix with columns called right-singular vectors of M, and is a diagonal m × n matrix with non-negative real numbers on the diagonal. The right-singular vectors of M are a set of orthonormal eigenvectors of M † 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 www.nature.com/scientificreports/ 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 where the (0, . . . , 0) T are located at the last k Landau levels to be truncated at the up and down spin space separately. Both φ and χ are columns of size (L − k) × 1 and have support up to L k . Therefore, the weighting components of the {L k + 1, . . . , 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 × 2L matrix can be written as where the 2k columns filled with × denotes the columns to be truncated. For the eigenstate ψ , we can have Hψ = Eψ , in which we are interested in low energy E regime. Since components of the last {L − k + 1, . . . , L} levels for the eigenstates of interest are zero, The columns with × for the Hamiltonian actually have no effect. We can then drop them and collect the truncated Hamiltonian denoted as H to be which is 2L × 2(L − k) matrix. We then do the SVD factorization to have H = U V † , with each column vector of V denoted as which is a column of size 2(L − k) × 1 . Padding with zeros in the form of Eq. (9) for ψ to become ψ , we can have ψ as the eigenstate of H 2 with eigenvalue of E 2 . Doing some linear combination of eigenstates with the same eigenvalue = E 2 , we obtain the eigenstate ψ ′ satisfying the eigenequation Hψ ′ = Eψ ′ . 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. (9) 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 ψ constructed from ψ from the SVD can be eigenstate of H 2 is simple. The Hermitian conjugate of the Hamiltonian H is such that (9) ψ = φ , 0, . . . , 0,χ , 0, . . . , 0 T , 2g� (q) into the differential equations and obtain where 2 = 1 g 2 − ε 2 . Note that the exponential factor e − q 2 2g is to eliminate the quartic term q 4 . Assuming > 0 , we continue to take χ(q) = q f 1 (q) and � (q) = q f 2 (q) and have Then we rescale the length by q = q √ g , and obtain The final step is to change the variable from q to ρ =q 2 , which results in www.nature.com/scientificreports/ In the above, we already used Reformulating the last two equations, we have the so called associated (generalized) Laguerre equations where n = 1 2g − 2 = 0, 1, 2, . . . . The solutions f 1 and f 2 will be the associated Laguerre polynomials: f 1 (ρ) = L n−1 (ρ) and f 2 (ρ) = L n (ρ) with n = 0, 1, 2, . . . . For n = 0 , we have to take f 1 = 0 . As n = 1 2g − 1 2 1 g 2 − ε 2 ≥ 0 is a non-negative integer, the energy eigenvalues are here E n ≤ 1 g as 2 = 1 g 2 − ε 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, �) T with where L n is a polynomial function of q 2 to degree n. Therefore, the normalizability demands ≥ 1 2 , that is, 0 ≤ n ≤ 1 2 1 g − 1 2 . 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 α = 0 which omit the k y and k z dependence in the σ x term, since it produces an open equienergy contour for E ≥ E VH and hence the cyclotron orbit is not confined. The prime means the system in a magnetic field. We could try solutions either as � = (0, ψ) T or � = (ψ, 0) T , but it turns out that the second choice is not normalizable. Then we are going to solve the differential equation

Analytical form of zero energy solutions for Model
Its large-q limit can be conquered by setting ψ(q) = e −κq 2 �(q) and taking it into Eq. (22) to obtain κ . As κ > 0 for normalizability, we have κ = 1+ , 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 κ = 1+ √ 1−α gα and therefore we focus on the plus sign choice. It follows the differential equation for �(q) , which is and θ = ξ q , the equation is then transformed into the Hermite differential equation The solution for is 1 F 1 − 1 4 ; 1 2 ; θ 2 and H /2 (θ)) , where 1 F 1 (a; b; x) and H µ (x) are the Kummer confluent hypergeometric function and the Hermite polynomial, respectively.
It is known that H /2 (ξ q) is not purely even or odd with respect to q. As one can find that the even part of H /2 (ξ q) is actually 1 F 1 − 1 4 ; 1 2 ; ξ 2 q 2 , it is consistent that after reconstruction the zero-energy eigenfunctions are either of even or odd parity in q. (16)  www.nature.com/scientificreports/ 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.

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 A and H B to k z = k V and include the other pair of Weyl nodes at k z = −k V , we can have the modified models, denoted as HÃ and HB , in the form of The Weyl nodes of the first pair are then located at (±k W , 0, k V ) , while the other pair is located at (±k W , 0, −k V ) . Since we are looking at physics near the Weyl nodes and their low energy spectrum, the αk 2 term does not play much role. In most of the time they can be even dropped. The value of α under discussion is therefore small, and the αk 2 term affects mainly the dispersion in higher energy and does not influence the low energy of interest much.
Here Model Ã and Model 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 = σ 0 or σ 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 W , i.e. k V > k W . Due to the other pair of nodes, the term of q z σ z in Model A Hamiltonian under magnetic field H ′ A is replaced by . 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 Ã is then while the Hamiltonian for Model B is then 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 www.nature.com/scientificreports/ discussions of model H B and HB . 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 αk 2 have the same coefficients for k y and k z while linear terms ∼ v y k x k y σ y and ∼ v z k x k z σ 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 α can be small but cannot be zero. For single pair of Weyl nodes, the Hamiltonian is which is already defined in the main text. The coefficient α mainly influence dispersion in higher energies. Usually α 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 ′ y and k ′ z where the new ẑ direction is along the field. Suppose the field B = B(sin θŷ + cos θẑ) , then the new momentum coordinates are defined as k ′ y = cos θk y − sin θk z and k ′ z = sin θk y + cos θk z , and still k 2 In such definition, we take the advantage of k ′ z still being a good quantum number and k x k ′ y → l −2 B (k xx + i 2 ) similar to model H B with modified guiding center x 0 = l 2 B k ′ y . With definition of σ ′ y = cos θσ y − sin θσ z and σ ′ z = sin θσ z + cos θσ z , the Hamiltonian in new coordinates can be written as The dimensionless momentum in the field direction is defined as Therefore, the Hamiltonian under field in rotated coordinate is then where the prime on the left side refers to the Hamiltonian under magnetic field. The independent parameters are α , g, and rotated angle θ . Among them, α and g are related to materials properties, i.e. Weyl nodes quantities, and g and θ are related to field amplitude and direction respectively. As an example, we present the case of α = 0.05 and with fixed value of g = 0.8 , we rotate the angle θ of field with respect to ẑ axis from 0 • to 180 • and present the result in Fig. 5a. It can be found that the zero energies persist in all angles θ in the plane of k ′ 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 = ±k V , the Hamiltonian for Model A and Model B are Eq. (25), and we focus on discussing Model B.
Similarly, when written in the new coordinate (k x , k ′ y , k ′ z ) of the rotated frame, Model B Hamiltonian is then www.nature.com/scientificreports/ Different to H B , the dimensionless momentum along the field direction is defined as q z = k ′ z /k V . Therefore, the dimensionless parameters involving k ′ z would change the dependence from v x /v to k V /k W . Explicitly, the dimensionless parameters are defined as follows.
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 α , k W /k V , and v x v y . The Hamiltonian under magnetic field to solve is then With the reasonable choice of α = 0.05 , v y /v x = v z /v x = 0.5 , and k W /k V = 0.2 , we present the case of g = 0.8 and θ = 30 • in Fig. 5b. The corresponding values are α y = 0.2 , α 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.

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.
Model B. To model a vacuum-semimetal interface, we introduce a mass term in the Hamiltonian as where M(z) = 0 for z < 0 and M(z) = M → ∞ 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: where Re < and Re > are positive. For z > 0 , in the limit of M → ∞ , we have, by taking the ansatz into H B ψ = 0 and neglecting small numbers, leading to Since > > 0 , we have u = v for k x > 0 and u = −v for k x < 0 . As a result, we have boundaries conditions as With these boundary conditions, we take the ansatz for z < 0 into H B ψ = Eψ and we have, for k x > 0, Equating the real parts and the imaginary parts separately, we have E = k x k y and In order to have Re < > 0 , k x is limited by k x < k 2 W + αk 2 y . Similarly, for k x < 0 , we have E = −k x k y when − k 2 W + αk 2 y < k x . In conclusion, when putting back omitted constants, the surface states (the Fermi arcs) survive in |k x | < k 2 W + αk 2 y and take energy with energy E = −v � 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 ∝ k y . However, in Model B the surface-state wave function is not continuous at k x = 0 , which is proportional to (1, i) T on one side and (1, −i) 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 ∝ |k x |k y .
Annihilation of Weyl nodes. The pair of Weyl nodes move toward each other if we tune the parameter k 2 W 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 2 W < 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 Before collision where k W > 0 , the Weyl nodes are located at (±k W , 0, 0) . By tuning k 2 W to become negative, Model A has gap E g = 2|k W | at k = (0, 0, 0) or gap of E g = 4α|k 2 W |−1 α at ring of k 2 y + k 2 z = 2αk 2 W −1 2α 2 for k x = 0 and other positions are also gapped. Model A is fully gapped unless accidental cases like |k 2 W | = 1 4α . 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 = |k 2 W | α in the k x = 0 mirror plane.
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 M x for Model A and Model B are known, in conjunction with combined symmetry C 2 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 σ 0 , k x σ x , k x σ 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 σ y . The mirror symmetry breaking perturbation is therefore �k x σ y where the is small perturbation determined by the strength of magnetic field. As usual in the k 2 W > 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 ± k 2 W 1−α� 2 , ∓� k 2 W 1−α� 2 , 0 . As for the annihilation results (41) ψ(z < 0) = B e ik x x e ik y y e B z 1 √ 2 x − k 2 W − α(k 2 y + k 2 z )]σ x + k x k y σ y + k x k z σ z . www.nature.com/scientificreports/ by tuning k 2 W < 0 to make Weyl nodes to collide, the system is still fully gapped where the gap E g = 2|k W | at the origin and E g = √ , 0 , where the relative sign of k x and k y are determined by the sign of . After WNs annihilation when tuning k 2 W < 0 , the system remains gapless but no longer existing a nodal ring. The zero gapless positions are located at k = ± , 0 .
Here we demonstrate the differences between Model A and Model B in their behaviour of WNs annihilation when the parameter k 2 W are tuned from positive to negative. Although both of them describe the pair of WNs when k 2 W > 0 , Model A generally will be gapped out when k 2 W < 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.
Allowed forms of spin operators. First of all, we have to point out that the Pauli matrices in the Hamiltonians H A and H 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 pseudospin 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 M x = σ x . Under the mirror reflection, momentum and spin change as k = (k x , k y , k z ) � → (−k x , k y , k z ) and s = (s x , s y , s z ) � → (s x , − s y , − s z ) . At the same time, the pseudo-spin changes as (σ x , σ y , σ z ) → (σ x , − σ y , − σ 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: 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 M 2 T . Suppose that the k z is relative to k z = 0 or π , C 2 T makes (k x , k y , k z ) → (k x , k y , − k z ) , (s x , s y , s z ) → (s x , s y , − s z ) . Many antiunitary operators could be used for C 2 T with the restriction that either for spin-0 or spin-1/2 systems. However, when we refer to the transformation we find that it has to be C 2 T = σ x K , where K is the complex conjugation operation. With C 2 T , the spin will reduce its compositions as follows: With the spirit, we can obtain spin from the pseudo-spin for Model A too. Here M x = σ 0 and C 2 T = σ x K , we show them as (47) s x = σ 0 , k y σ 0 , k z σ 0 , σ x , k y σ x , k z σ x , k x σ y , k x σ z , s y = k x σ 0 , k x σ x , σ y , k y σ y , k z σ y , σ z , k y σ z , k z σ z , s z = k x σ 0 , k x σ x , σ y , k y σ y , k z σ y , σ z , k y σ z , k z σ z , (48) (C 2 T ) 2 = 1, s x = σ 0 , k y σ 0 , σ x , k y σ x , k x σ y , s y = k x σ 0 , k x σ x , σ y , k y σ y , k z σ z , s z = k z σ y , σ z , k y σ z .