Chirality flip of Weyl nodes and its manifestation in strained MoTe$_2$

Due to their topological charge, or chirality, the Weyl cones present in topological semimetals are considered robust against arbitrary perturbations. One well-understood exception to this robustness is the pairwise creation or annihilation of Weyl cones, which involves the overlap of two oppositely charged nodes in energy and momentum. Here we show that their topological charge can in fact change sign, in a process that involves the merging of not two, but three Weyl nodes. This is facilitated by the presence of rotation and time-reversal symmetries, which constrain the relative positions of Weyl cones in momentum space. We analyze the chirality flip process, showing that transport properties distinguish it from the conventional, double Weyl merging. Moreover, we predict that the chirality flip occurs in MoTe$_2$, where experimentally accessible strain leads to the merging of three Weyl cones close to the Fermi level. Our work sets the stage to further investigate and observe such chirality flipping processes in different topological materials.


I. INTRODUCTION
Weyl semimetals are arguably the most robust form of gapless topological matter [1][2][3][4][5][6][7][8][9]. They host point-like degeneracies of their energy bands, called Weyl nodes (or cones), whose low energy, linear dispersion relation is similar to that of elementary particles called Weyl fermions [10]. These band touching points are not accidental, but are a manifestation of the topologically nontrivial character of the semimetal phase. Each node has a chirality associated to it, a topological charge given by its Chern number, which means that an isolated Weyl cone cannot be gapped out by any perturbation. Instead, it can only be moved in energy and momentum space, or tilted [11]. This property sets Weyl semimetals apart from the many other types of gapless topological matter, in which degenerate points [12][13][14][15][16][17] (or lines [5,18]) require additional symmetries in order to remain protected.
Provided that Weyl cones are not isolated, it is possible to change their number: they can be created or annihilated pairwise when they overlap in energy and momentum. Similar to electrodynamics, these two-node processes require 'topological charge conservation,' meaning that the total Chern number of all band touching points in the Brillouin zone (BZ) must vanish [1]. The pairwise creation and annihilation of Weyl cones has been studied extensively, especially as a method to engineer Weyl semimetal phases by applying perturbations such as strain [19][20][21][22], magnetic field or changes in magnetization direction [23][24][25][26], disorder [27][28][29][30], phonons [31], or high-frequency illumination [32,33] to real materials.
In theory, there is no constraint limiting the number of Weyl nodes that can merge at a given point in the BZ. In practice, however, processes involving the simultaneous overlap of three or more cones of different chirality are highly improbable, especially when each node is allowed to occupy any point in energy and momentum space. To date, multi-Weyl merging has been studied mainly in the context of rather exotic types of topological semimetals [14,34], hosting 'unconventional fermions.' There, the symmetry-protected band degeneracies may split into multiple Weyl nodes if very particular symmetries are broken [14,35,36].
Here we study simple three-node processes occurring purely inside a Weyl semimetal phase, without unconventional fermions: a single Weyl cone overlaps with a pair of oppositely charged ones, causing the latter to disappear from the band structure. By topological charge conservation, the result is a single node with a flipped chirality. Interestingly, while lattice symmetries are not responsible for the existence of the degeneracy points, they play a fundamental role in enabling the three-node process to occur. This is because they restrict the relative positions and charges of the nodes in the BZ (for instance to planes or lines). Throughout the following, we will focus on systems obeying a combination of twofold rotation and time-reversal symmetry.
One of the advantages of studying multi-node processes in Weyl semimetals as opposed to more exotic types of gapless topological matter is that Weyl semimetal phases are much more abundant in real materials [37]. In fact, we find that one of the earliest predicted Weyl materials, MoTe 2 [38], hosts very close to the Fermi level a threenode process upon applying small uniaxial strain. As a result of this process, a Weyl node in a high-symmetry plane flips its chirality. This shows that the three-node process is within reach of photoemission experiments. Further, it shows that strain in real materials can act as a 'chirality switch', allowing to tune the transport properties associated with gapless topological systems.

II. CHIRALITY FLIP PROCESS IN SIMPLE MODELS
We start by discussing how the chirality of a Weyl node can flip using toy models. Given the requirement of topological charge conservation and considering only simple Weyl cones, chirality flips are only possible if at least two other, oppositely charged Weyl nodes are involved. For concreteness, we consider processes that start from three Weyl nodes with Chern numbers C = +1, −1, −1, and end with a single Weyl node having C = −1 [39]. We find that there are two scenarios through which this can happen: one that involves a three-node process and one that only involves two-node processes.
We begin by discussing the first scenario, corresponding to a three-node process, meaning that all three Weyl cones simultaneously overlap at the same point in the BZ. While this merging could in principle occur at generic momenta, its likelihood can be greatly increased when symmetries constrain the relative positions and charges of the band touching points. As mentioned above, we will consider a continuum model invariant under the combination of twofold rotation and time-reversal symmetries, with Hamiltonian Here, k = (k x , k y , k z ) is the quasimomentum and σ i are Pauli matrices encoding the degree of freedom associated with the two bands. For simplicity, we will set units such that the Hamiltonian and the quasimomentum are dimensionless. The Hamiltonian Eq. (1) obeys a twofold rotation symmetry along the k x = k y = 0 axis, with operator C 2 = −iσ y , as well as time-reversal symmetry, T = iσ y K (where K is complex conjugation), such that Their combination, C 2 T = K implies that and that the Hamiltonian is real for k z = 0. The constraint Eq. (4) means that if a Weyl cone is present in the C 2 T invariant plane, k z = 0, then it cannot exit the plane, due to topological charge conservation. Further, if there is a Weyl node at k z > 0, another one must be positioned symmetrically at k z < 0. Note that pairs of out-of-plane cones must have the same Chern number, since both time-reversal and twofold rotation are charge preserving.
For α < 0, H 1 hosts three Weyl cones, a positive chirality node at the Γ point, k = 0, and two negative chirality nodes at k = (0, 0, ± |α|). In contrast, for α > 0 there is a single Weyl node with negative chirality at Γ. As shown in Fig. 1(a, b), changing the parameter α continuously from negative to positive values causes the three nodes to merge at the same point, such that the central node flips its chirality. The three-node process occurring in the band structure of H 1 is protected by symmetry: the Weyl merging must involve all three nodes simultaneously. This is because the Hamiltonian obeys C 2 and T separately, and because all Weyl nodes are located on the twofold rotation axis. The band touching points cannot move away from k x,y = 0 without breaking the rotation symmetry, and they must be positioned symmetrically with respect to k z due to time-reversal symmetry. The only allowed overlap is thus at Γ.
The presence of a chirality flip in H 1 is also related to topological quantities. Due to the presence of C 2 T , the Hamiltonian is real at k z = 0. Therefore, a closed momentum-space loop in this plane will be characterized by a π-quantized Berry phase, provided that the loop encircles the central node. This Z 2 topological invariant indicates that the parity of nodes inside the loop must be conserved [40,41]. Therefore, the central node is not allowed to leave the plane, but it is allowed to reverse its chirality.
If the system instead hosts Weyl cones away from the rotation axis, or if it obeys only the combined C 2 T symmetry but not its individual components, the three-node process is not symmetry protected anymore. This does not mean that it is forbidden. Even if it is not forced to occur, it still can, and as we show later actually does happen in real materials.
We now move on to discuss the second scenario for the chirality flip, one that does not involve a three-node process. We illustrate this scenario in Fig. 1(c), which shows the Weyl points of the Hamiltonian Note that H 2 obeys the same constraint as H 1 in Eq. (4), but has neither rotation nor time-reversal symmetries. For α < 0 there are three Weyl nodes in the band structure of H 2 : one in the C 2 T invariant plane k c = (k 0 , 0, 0), where k 0 is the single real solution of k 3 + k 2 = −|α|, and two away from the C 2 T plane at k s = (0, 0, ± |α|). Due to C 2 T , the outside nodes can only enter the plane at the same point, as before. However, now they can meet without the third node, thus forming a double Weyl point. In H 2 , the double Weyl point occurs at k z = 0 when α = 0. For α 0 the overlapping nodes separate again, and all three cones are in the C 2 T plane, located at (k i , 0, 0), where k i are the real solutions of k 3 +k 2 = |α|. The three nodes remain until α = 4/27, at which point two of the opposite chirality Weyl points annihilate. At α > 4/27, we are left with a single Weyl node. The trajectory of all three Weyl cones is shown schematically in Fig. 1(c). Similar to the behavior of nodes in the spectrum of H 1 , going from negative to positive α converts a positive chirality Weyl point into a negative chirality one.
As before, considering a momentum-space loop in the C 2 T plane that encloses the projections of all nodes in Fig. 1(c), we find a π-quantized Berry phase for all α. This indicates that the two-node and three-node process cannot be distinguished topologically by using the Berry phase. This is expected, as one can imagine continuously moving the meeting point of the out-of-plane nodes until this point overlaps with the k z = 0 node. However, when considering an infinitesimally small loop around the initial red node, a blue node will have to cross this loop for some α, thus ensuring a discontinuity in the Berry phase.
We note that C 2 T -symmetry also relates to the Euler class [42][43][44][45][46][47][48], which quantifies the obstruction for pairs of Weyl nodes to leave the C 2 T plane. For a generic, multiband system obeying C 2 T symmetry, this obstruction is expressed by an additional set of topological invariants, independent of the Chern numbers, which characterize Weyl nodes inside the k z = 0 plane. These additional invariants are called "non-Abelian frame charges," and we denote them by q. A pair of nodes with opposite non-Abelian charges have an Euler class of zero, and a pair of nodes with the same non-Abelian charge carry a nonzero Euler class. Importantly, if the Euler class is nonzero then the two nodes are not allowed to exit the plane, neither by merging and moving away from k z = 0 nor by annihilating with each other.
In our two-band models, these additional invariants q are given simply by the winding number of the vector multiplying the σ-matrices [45], and take values q = ±1. For the two-node process shown in Fig. 1(c), the Euler class has several consequences. First, the out-of-plane nodes (blue) must acquire opposite non-Abelian charges (q = +1 and q = −1) once they enter the plane, reflecting the fact that this pair is allowed to leave the C 2 T plane again, through the reversed process. Further, the annihilation of the initial in-plane node (red) with one of two nodes that entered the plane tells us that they must have opposite non-Abelian charges. Therefore, the non-Abelian charge q of the remaining node must be the same as that of the initial in-plane node at the beginning of the process. In other words, even though the Weyl cone flips its chirality, the non-Abelian charge associated to it remains conserved. There could be no annihilation if the nodes of equal chirality (blue) were swapped, since the two merging nodes would then have the same value of q.

III. TRANSPORT PROPERTIES
The initial and final stages of both chirality flipping processes are the same. This means that the observables for large |α| have the same qualitative behavior for the two processes. The difference in the two scenarios is prominent close to α = 0 where the two Weyl nodes outside of the C 2 T invariant plane reach the plane.
We study the differences between the two processes by computing the conductivity using semiclassical Boltzmann transport theory in the constant relaxation time approximation. The details of the calculations are shown in the supplemental material (SM) [49]. At large |α| and low energies, both systems show a similar behavior: the conductivity of isolated Weyl nodes is σ jj ∝ µ 2 [50], with j = x, y, z and µ the chemical potential measured relative to the Weyl node energy. At α = 0, however, the two systems are significantly different. In the three-node process we have a C = −1 node, with linear dispersion in k x and k y and cubic dispersion in the k z direction. In the twonode process we have a double Weyl point with C = −2 and with quadratic dispersion in k x and k z , but linear dispersion in k y . This difference in the dispersion at low energies results in different conductivities as a function of chemical potential. At α = 0 the three-node process is characterized by whereas the two-node process shows Note that in the latter case we only focused on the double Weyl node forming at k z = 0. The third Weyl node will contribute additively to the conductivity [49]. These differences are prominent when looking at the conductivities in the plane of the three nodes (σ xx and σ zz ). In Fig. 2 we plot the ratio σ xx /σ zz as a function of α for different chemical potential values. Panel (a) shows the three-node process, in which the ratio is enhanced close to µ = 0 and α = 0, which is due to the different dispersions in the two directions. On the contrary, panel (b) shows that σ xx /σ zz = 1 at α = 0 for the two-node process, because of the identical dispersion relations in the two directions.

IV. WEYL NODES IN STRAINED MoTe2
Many materials reported to host Weyl cones have C 2 T symmetry, and in several cases Weyl nodes are placed inside the C 2 T -invariant plane, including {W,Mo}Te 2 [38,51], Ta{Ir,Rh}Te 4 [52] and ZrTe [53]. Under a suitable external perturbation, these compounds are natural candidates for the three-node process described above.
As an example, we show based on density-functional calculations that uniaxial strain in MoTe 2 controls a process of the type sketched in Fig. 1(b). We perform fully relativistic calculations as implemented in the FPLO code [54]. We use as equilibrium lattice parameters those reported in Ref. [38] at 100 K (a = 3.468 A, b = 6.31 A and c = 13.861 A). We focus on uniaxial strain such that the lattice parameter a is enlarged. For a fixed deformation δ a (measured in percent of the original value of a), we determine the deformations of b and c according to the Poisson ratios γ ab = 0.19 and γ ac = 0.96 reported in Ref. [55]. Previous works have shown that treating the Mo-4d shell with the GGA+U method improves the description of photoemission experiments provided by GGA [31,51,56,57]. Thus we here use the former method with J = 0 and U = 2.6 eV. Further details, including the robustness of our results to all these choices are presented in the SM [49]. Fig. 3(a) shows the band structure of MoTe 2 without strain (δ a =0) along a path that includes Weyl nodes connecting the higher valence and lower conducting bands (W v,c 1 and W v,c 2 ) and nodes between the two upper valence bands (W v,v 1 and W v,v 2 ). The results of interest in this work originate in the latter two which are found 25 meV (W v,v 1 ) and 34 meV (W v,v 2 ) above the Fermi level. Note that both are type-II, or over-tilted nodes, unlike the band touching points of the toy models presented above. W v,v 1 has seven other partners in the BZ, which together form a set of eight nodes at generic momentum coordinates, which are related to each other by the reflection symmetries M x and M y , the twofold rotation around the c axis, C 2 , and time-reversal T . On the other hand, there are only four nodes of the type W v,v 2 , all of which are positioned inside the C 2 T -invariant plane [ Fig. 3(b)]. As δ a increases, the energy difference between W v,v 1 and W v,v 2 decreases and vanishes at the crit- There are a total of four such processes throughout the entire BZ, one at each quadrant of the (k x , k y ) plane, due to the presence of reflection symmetries. In a given quadrant, before the Weyl merging, one finds a pair of nodes at finite k z and Chern number C = −1, and a single node having k z = 0 and C = 1. After the merging, only a single Weyl node having C = −1 and k z = 0 exists [ Fig. 3(d)]. Thus, MoTe 2 realizes a three-node process of the type discussed previously.
Note that none of the nodes are on the twofold rotation axis, so the triple-Weyl merging is not forced to occur, i.e., it is not protected by symmetry. Instead, our results suggest that the three-node process is favored energetically, since it involves a smaller change in the local gap between the bands forming the Weyl cones (see SM [49]).

V. CONCLUSION
We have shown that the topological charge associated to Weyl cones can change sign. This chirality flip involves three nodes, and may occur generically in Weyl semimetal phases. We discussed the two generic mechanisms by which this can happen: a three-node process where a triplet of Weyl cones merge at the same point in the BZ, and a two-node process, which occurs by means of successive, pairwise mergings of nodes. Additional symmetries, in our case time-reversal and twofold rotation, increase the likelihood of chirality flips by constraining the positions and relative charges of Weyl nodes. In some cases, symmetries may even enforce the occurrence of the three-node process.
Our results indicate that the chirality flip occurs in one of the most well-studied Weyl materials, MoTe 2 , where moderate uniaxial strain leads to a simultaneous merging of three nodes close to the Fermi level. Depending on the doping level of different MoTe 2 samples, this indicates that chirality flips are within reach of photoemission and transport experiments. It could be interesting to extend our study of electronic transport to models relevant to MoTe 2 , in particular to magnetotransport properties sensitive to the chiral anomaly.
More generally, our work sets the stage for further investigating such chirality converting processes in a wide range of materials and scenarios. As we have mentioned, there is an abundance of Weyl materials obeying C 2 T symmetry, many of which host Weyl cones pinned to the C 2 T -invariant plane. Furthermore, similar processes could be also observed in topological metamaterials that host Weyl nodes such as interacting spin systems [58,59] or multiterminal Josephson junctions [60,61].
An interesting direction for future work could be to examine the behavior of chiral topological metals, which may obey time-reversal and rotation symmetries, but in which inversion symmetry is strongly broken. In these systems, Kramers' theorem guarantees the presence of Weyl points at the time-reversal invariant momenta (TRIM) of the BZ. When TRIM points lie on a rotation axis, we have shown that chirality flips must occur via a three-node process.

ACKNOWLEDGMENTS
We thank Ulrike Nitzsche for technical assistance and János K. Asbóth for useful discussions. In this supplemental material we discuss the details of the results presented in the main paper. We give a detailed explanation on the calculation of the density of states and conductivity. We present the details of the ab-initio calculation. We provide relevant input files and data in the data repository [1].

DENSITY OF STATES AND CONDUCTIVITY OF THE SIMPLE MODELS
Here we give details on the calculation of the density of states (DOS) and conductivity of the simple models introduced in the main text. The two models are The eigenvalues of these Hamiltonians are The dispersion relation of these models together with the Berry curvature are shown in Figs S1 and S2. For simplicity we will focus on the part of the processes where the two Weyl nodes at k z = 0 enter the C 2 T -invariant plane. In the H 2 case this means that we ignore the third Weyl node in the plane and use β = 0.
The density of states is calculated using where δ is the Dirac delta function. Since the spectrum is symmetric with respect to E = 0, we only discuss D(E ≥ 0). We start with the H 1 system. The numerically calculated DOS is shown in Fig. S3.
In the α > 0 case at low energies the linear term dominates and we get an anisotropic Weyl node with velocity α in the z direction. In the special case of β = 0 we have linear dispersion up to arbitrary energies and the DOS can be calculated analytically: This result matches the curves in Fig. S3b when α > β and the energy is sufficiently low (the darker curves). At α = 0 we are exactly at the chirality-flip point. Here in the z direction the linear term vanishes and the cubic term becomes dominant. In the case of α = 0 the DOS can again be calculated analytically and we get This exact solution matches perfectly the yellow curve in Fig. S3b. At α < 0 there are three Weyl nodes. At low enough energies we see the quadratic increase of the DOS as in Eq. (S6). Once the energy reaches the saddle points in between Weyl nodes we get a kink in the DOS as seen in Fig. S3a. The larger the separation between the Weyl nodes the higher the energy where the kink happens. Thus, for large negative α the qualitative behavior of the DOS is D ∝ E 2 in a large energy range. Close to α = 0, however, we get a qualitatively different behavior, with D ∝ E 4/3 .
We continue with the H 2 system at β = 0. The numerically calculated DOS is shown in Fig. S4. Since the system is symmetric with respect to changing α → −α and k x ↔ k z the DOS will be symmetric in α. In the case of α = 0 at low energies we always have two Weyl nodes and in the DOS we see the quadratic dependency on energy as in Eq. (S6). Similarly to the H 1 system we see kinks in the DOS at the energies of the saddle points in between Weyl nodes. The major difference between the DOS of H 1 and H 2 is at α = 0 where the two Weyl nodes outside the C 2 T invariant plane reach the plane. It can be shown that the DOS for H 2 is where γ is an energy-independent integral The numerical value of this integral is γ ≈ 0.05463. As we can see the DOS is proportional to energy, in contrast to the D ∝ E 4/3 dependence we found for the three-node process. We continue with the calculation of the conductivity of the above Hamiltonians. To get a qualitative picture on how the conductivity behaves as a function of the chemical potential we use the semiclassical Boltzmann transport theory with relaxation time approximation [2]. We also assume zero temperature and use the relaxation time at the Fermi energy as a parameter τ ≡ τ (µ). The DC conductivity is calculated as The numerically calculated conductivity for the H 1 system is shown in Fig. S5.
In the α = 1, β = 0 case we have a single isotropic Weyl node and for the conductivity we get This is the usual result for the conductivity of a single Weyl node [3]. As it was discussed in Ref. [3], the Boltzmann approach strictly holds for µ > Γ where Γ = /2τ is the scattering rate. Close to the Weyl point energy a more detailed description would be necessary, based on the specific type of scattering effects [4,5]. Since we are interested only in the qualitative behavior we limit ourselves to the simple method of constant relaxation time approximation. At a general α with β = 0 we have an anisotropic Weyl node and the conductivity can be calculated as We can see this behavior in Fig. S5 at α > 0 and α β. At the chirality flip point α = 0 the conductivity is σ zz (µ) = 9 22π 2 e 2 τ β 1/3 µ 8/3 . (S15) At α < 0 we have three Weyl nodes and at small enough chemical potentials the total conductivity is the sum of conductivities of the three nodes. The conductivity will be anisotropic but the chemical potential dependence will be quadratic in all directions as in Eq. (S11). This can be seen in Fig. S5c where the ratio of σ xx to σ zz is shown. Close to α = ±1 the ratio is independent of the chemical potential. Close to the chirality flip the anisotropy becomes chemical potential dependent, and the ratio becomes enhanced at low chemical potentials. We continue with the calculation of the conductivity of the H 2 system at β = 0. The numerically calculated conductivity for the H 2 system is shown in Fig. S6. At β = 0 because of the symmetries of the system σ xx (α) = σ zz (−α). For large |α| we have two Weyl nodes and the conductivity has a quadratic chemical potential dependence as in Eq. (S11). At α = 0, σ xx = σ zz . It can be shown similarly to the DOS that σ xx/zz (µ) = λ 1 e 2 τ µ 2 , (S16) σ yy (µ) = λ 2 e 2 τ |µ|, where λ 1 and λ 2 are given by where η 0 was defined in Eq. (S9). The numerical values of these integrals are λ 1 ≈ 0.04222 and λ 2 ≈ 0.02731. The chemical potential dependence of the conductivity for large |α| is qualitatively the same for the H 1 and H 2 Hamiltonians, as expected. Close to α = 0, however, the chemical potential difference is qualitatively different, with different power laws in the two cases.

AB-INITIO RESULTS
We performed fully-relativistic calculations as implemented in the FPLO code version 18.57, using a mesh of 18 × 12 × 6 subdivisions for Brillouin zone integrations. Fig. S7a shows the band structures obtained with GGA and with GGA+U . In the latter case, we use J = 0 and U = 2.6 eV. In agreement with Ref. [6], the most prominent effect of including U is the elimination of the electron pockets at Y .
Tight-binding Hamiltonians were obtained by constructing Wannier functions associated with the Mo-4d and Te-5p orbitals, using the so-called projecting method. The eigenenergies computed with these models typically differ less than 5 meV from the DFT results. Starting from the Wannier model, Weyl nodes were obtained with a bisection-like algorithm that searches for Berry curvature sources, implemented in Ref. [9]. The algorithm provides a search of candidate Weyl nodes. Second, in order to establish or discard Weyl node candidates, we compute the corresponding Chern number by integrating the Berry curvature around a small sphere centered at the Weyl node candidate position. Table I lists relevant information of the Weyl nodes presented in the main text at zero strain. As an additional check, in certain cases we have numerically confirmed the existence of an apparent band-crossing directly in the DFT results by computing the eigenenergies in a sufficiently dense mesh of k-points near the Weyl node.
For the calculations of uniaxial strain presented in the main text, we first took steps of δ a of 0.03%. After that, we made a refinement close to the critical point in steps of 0.006%. Very close to the critical point, artifacts of the numerical resolution become evident, e.g., leading to numerically-calculated Chern numbers that deviate significantly from integer values and the corresponding data has been excluded. To explore the robustness of the predicted three-node process as a function of strain, we have repeated the search of Weyl nodes by varying the value of U and the values of the equilibrium lattice parameters. Fig. S7(b-e) show the results of different tests. Fig. S7b show data based on GGA+U with U = 2.6 eV and with equilibrium lattice parameters from Ref. [7] (a = 3.468 A, b = 6.31 A and c = 13.861 A ). Fig. S7c presents data based on the equilibrium lattice parameter reported in [8] (a = 3.477 A,b = 6.335 A and c = 13.883 A). Variations of this order do not lead to visible changes in the Weyl-node dynamics. Fig. S7d and e present results based on U = 2.5 eV and plain GGA (U = 0), respectively. The three-body chirality-flip process takes place also for these parameters. The precise value of U affects the critical value of strain at which the three-body process takes place.
In order to further confirm that the process we observe is indeed the three-node process, we focus on the Weyl node in the C 2 T -invariant plane. If the chirality flip is carried out through two-node processes the annihilation of the W v,v 2 node with one of the W v,v 1 nodes must happen in the C 2 T plane. In such a case the velocity difference between the two bands forming the Weyl node in the C 2 T -invariant plane must vanish at the point of annihilation. To exclude this possibility we calculated the minimal, in-plane velocity difference of the bands forming the W v,v 2 node. For this analysis, we start from two ab-initio constructed Wannier Hamiltonians, one at equilibrium δ a = 0 (H 0 ) and one corresponding to δ a = 0.6% (H 0.6 ), a value larger than the critical strain at which the three-body process occurs (∼ 0.5%). Then, we introduce a parameter λ that linearly interpolates between the two: H[λ] = (1 − λ)H 0 + λH 0.6 . The results as a function of λ are shown in Fig. S8. As we can see no annihilation occurs in the plane, which would be indicated by a dip of the velocity difference. Instead, we observe a monotonous increase. node, computed in the C2T plane, as a function of strain. λ is a parameter that linearly interpolates between the Wannier Hamiltonians at δa = 0 and δa = 0.6, values of strain below and above the critical strain (∼ 0.5%). The velocity difference is measured based on the the energy difference on a circle with radius k = 0.001 · 2π/a around the Weyl node as v = ∆E/k.
Finally, to gain more understanding on why MoTe 2 exhibits a three-node, as opposed to a two-node process, we plot in Fig. S9 the momentum dependent energy gap in the plane containing the three merging Weyl cones. We observe that the smallest energy difference between bands occurs on an arc connecting all three nodes. This suggests that the three-node process (solid arrows) is energetically favored compared to a hypothetical two-node process, for instance as indicated by the dashed arrows. This is due to the fact that having the out-of-plane nodes W v,v 1 merge at k z = 0 without the in-plane node, W v,v 2 , would require a larger change in the gap. nodes upon increasing strain. The dashed lines indicate a hypothetical path that would achieve the two-node process.