Sixfold fermion near the Fermi level in cubic PtBi2

We show that the cubic compound PtBi2, is a topological semimetal hosting a sixfold band touching point in close proximity to the Fermi level. Using angle-resolved photoemission spectroscopy, we map the bandstructure of the system, which is in good agreement with results from density functional theory. Further, by employing a low energy effective Hamiltonian valid close to the crossing point, we study the effect of a magnetic field on the sixfold fermion. The latter splits into a total of twenty Weyl cones for a Zeeman field oriented in the diagonal, [111] direction. Our results mark cubic PtBi2, as an ideal candidate to study the transport properties of gapless topological systems beyond Dirac and Weyl semimetals.


Introduction
Topological semimetals are materials whose bandstructure contains topologically protected intersections of energy bands. Well-known examples include Dirac and Weyl semimetals, that host point-like, cone-shaped band crossings [1][2][3][4]. They are governed by low energy effective Hamiltonians which take the same form as those describing their high energy physics namesakes, the elementary Dirac and Weyl fermions. Unlike in quantum field theory, however, Poincaré symmetry does not constrain the set of emergent quasiparticles allowed in condensed matter systems. As such, the low energy excitations of a crystal may come in flavors that go beyond those of the standard model of particle physics [5][6][7][8][9][10][11][12][13][14][15][16]. In terms of the bandstructure, these unconventional fermions arise at the intersections of multiple energy levels, such as 3, 4, 6, or even 8-fold points.
Beyond their appeal as a playground that harbors aspects of the fundamental theories of high energy physics, topological semimetals might also have potential technological applications. For the latter, one aims to exploit the unique electrical properties induced by the crossing points, such as the chiral anomaly [17], or the axial gravitational anomaly [18]. Since these effects involve mainly states close to the Fermi level, their observation imposes strict constraints on the energy at which band intersections should occur, thus reducing the range of viable material candidates. Indeed, despite the relatively large number of confirmed Weyl and Dirac systems , few show band crossings in close proximity to the Fermi level. For semimetals hosting unconventional fermions the list is smaller still, especially given that to date, few such materials have actually been synthesized and checked experimentally for the presence of nodal fermions [42][43][44][45][46][47][48][49][50].
In this work, we identify cubic PtBi 2 as an ideal candidate for transport studies on topological semimetals. We use a combination of angle-resolved photoemission spectroscopy (ARPES) and density functional theory (DFT) to map the bandstructure of the system. We show that cubic PtBi 2 hosts a sixfold band touching point -a triple Dirac point -in close proximity (10-20 meV) to the Fermi level. To study the essential topological features of the crossing, we use a low-energy effective Hamiltonian which reproduces the ab-initio data. Upon including a weak Zeeman field to the model, we find that the sixfold point can split into a total of twenty type-II Weyl cones.

Theory
The DFT calculations where performed using the full 4-component relativistic mode of the Full Potential Local Orbital (FPLO) code [51] in version 15.02 with a 10 × 10 × 10 k-mesh for the Brillouin zone (BZ) integration and the local density approximation (LDA) [52]. The space group is 205 (P a3) with lattice parameter a 0 = 6.6954 and Pt at Wyckoff position ( 1 2 , 0, 0) and Bi at (0.12913, 0.12913, 0.12913). Following the self consistent calculation a maximally projected Wannier function model [53,54] containing all Bi 6p-and Pt 5d-orbitals was constructed. The corresponding core energy window reaches from −5 eV to −3 eV with a lower and upper Gaussian fall-off of halfwidth 1 eV and 3 eV attached, respectively, ensuring a very good fit to the DFT band structure in an energy window from −6 eV to +5 eV with a maximum error around and below the Fermi level of 30meV.
This Wannier model was used to calculate the spectral density for an infinite bulk system by k z -integrating the bulk band structure including an imaginary part of the band energies of 20meV ((001)-bulk projected bands [BPB]). Furthermore, a semi-infinite slab (SIS) with a triple layer BiPtBi as termination was setup to calculate spectral densities of an idealized (001)-surface (simple Hamiltonian cutoff at the surface). The spectral densities were obtained via a Green's function recursion method [55,56].
For the subsequent comparison with the photoemission intensity from the (001) surface of PtBi 2 we have calculated the spectral function corresponding to the bulk (BPB) and surface states (SIS) in a repeated zone scheme. Since the k z -resolution is a strongly materialdependent quantity and thus the degree of surface contribution is not a priori known, we present both momentum- (Fig. 1d, 1f) and energy-momentum (Fig. 1e, 1g) intensity distributions along high symmetry directions by either including (SIS) or excluding (BPB) the surface states. Full integration of the bulk states along the (001) direction (BPB) will allow us a rough estimate of the k z -resolution when comparing to experiment. Comparison of panels 1d) and 1e) with 1f) and 1g) respectively reveals the contribution of the surface states, as predicted by theory. Fig. 1b suggests that there are three types of the Fermi surfaces in PtBi 2 and all of them are strongly three-dimensional. The structure near the Γ-point is due to the projection of large electron-like Fermi surface (FS) (Fig. 1d and 1f). The ones in the corners of the BZ are the projections of eight electron-like Fermi surfaces united to the shape similar to a stellated octahedron centered in the R-points. Other twelve "beans" are the hole-like pockets, made by the crossings along ΓM and ΓA. They are responsible for the four larger filled ellipses close to the sides of the surface Brillouin zone and for the four smaller ellipses overlapping with the projections of the FS in the corners. Note the C 2 symmetry of the intensity pattern in Fig. 1d, which stems from the symmetry reduction of the non-symmorphic cubic space group due to the surface termination/k z -integration. For instance, filled ellipses are closer to each other then to Γ along ΓY while the opposite is true for the ΓX direction. The projection with the center at the M-point is clearly asymmetric too. The surface adds features which make the C 2 -symmetry even more pronounced (Fig. 1f). These are the "brackets" connecting the Γ-centered FS-projection with ellipses along ΓX as well as the vertical arc-like connections of the "petals" near M-point. The surface also adds bright and well localized intensity spots at Γ-and M-points.

Experiment
PtBi 2 single crystals were grown as described in Ref. [57]. Angle-resolved photoemission spectroscopy (ARPES) measurements were carried out at "1 3 -ARPES" facility at BESSY at a temperature of 1 K within the wide range of photon energies (60 -100 eV) from the cleaved surface of high-quality single crystals. The overall energy and momentum resolutions were set to ∼8 meV and ∼0.013Å −1 , correspondingly.
In Fig. 2 we show the ARPES Fermi surface maps recorded using photons of different energies and compare them with the calculated momentum distribution including surface states (SIS). All maps shown in Fig. 2 and in particular the pattern within the first Brillouin zone (white dashed box) are in qualitative agreement between each other, with theory ( Fig. 2c), as well as with a previous study focusing on the surface states of PtBi 2 [58]. All the FS projections discussed earlier are seen in the experimental maps having very similar shapes and intensity patterns, including the C 2 -symmetry. Moreover, the surface related features, such as "brackets" and arc-like connections near the corners of the BZ are clearly visible as well. The similarity of all experimental maps taken using different photon energies clearly implies a very low k z -resolution. The differences are mostly due to unequal intensity distribution rather than different peak positions of the spectral function. Therefore, one can consider any of the maps as representing the total spectral weight -fully integrated over k z bulk states plus surface states. Intensity variations are due to matrix elements highlighting a particular portion of the k-space.
Keeping this in mind, a closer inspection of the maps reveals certain discrepancies with theory. Ellipses are closer to each other in the experiment near Y-points and even overlap or at least touch there. They also seem to be larger than their near-X counterparts. Another difference is the absence of the sharp intensity maxima due to surface states at the Γ-and M-points.
To get a better insight into the origin of the observed similarities and discrepancies, we show in Fig. 3 the comparison of the experimental momentum-energy cuts with the corresponding calculations of the total intensity along the ΓY-and ΓX-directions. In Figs. 3 and 4 we use binding energy for the vertical axis, which has a sign opposite to one of energy in Fig. 1. As in the case with the FS maps, the overall agreement is very good. All the well-pronounced features in the calculations can be tracked in the experimental data. Even the intensity distribution is captured by the calculations correctly in most cases, but due to a big number of features their relative intensities are not always reproduced exactly. The most pronounced difference is the absence of the crossing of two dispersing features at the Fermi level at the Γ-point. These features can be clearly identified as surface states, as follows from the comparison of Fig. 1e and Fig. 1g. We noticed the corresponding absence of the bright spot in the Fermi surface maps earlier.
The data shown in Fig. 3b provides an explanation for this discrepancy. Indeed, there are two features close to the Γ point with the dispersion opposite to the electron-like parabola supporting the large Γ-centered Fermi surface, but they cross the Fermi level before they reach the Γ point. This is best seen when considering the left-most Γ point in Fig. 3b. To highlight these dispersions we have schematically extended them above the Fermi level using dashed lines. The corresponding set of features is seen also in all experimental maps in Fig. 2 inside the large Γ-centered FS. Also, one has to keep in mind that the idealized surface termination in the SIS calculations may not accurately reproduce the experimentally realized energy dispersion or the ones corresponding to a self-consistent slab calculation including all surface relaxations. Nonetheless, we note, that the other set of surface states, the one responsible for "brackets," is reproduced remarkably well: the linear dispersions starting approximately at 200 meV below the Fermi level at Γ are clearly seen in Fig. 3d for the left-most and right-most BZs.
In Fig. 4 we show the comparison of theory with experiment for the ΓM-high-symmetry direction. We have selected two representative cuts with different intensity distributions to highlight different parts of the spectra. Both measurements again show an overall agreement with the calculations. Also this time the most noticeable variations can be attributed to the surface states. The difference between the plots in Fig. 1e and 1g suggests that purely surfaceoriginated features (except the already discussed ones at Γ) are the parabolic dispersion with the top located roughly in between Γ and M at +150 meV and features dispersing between −100 and +500 meV and crossing near the Fermi level, thus highlighting the M-point in the calculated map (dashed lines in Fig. 1f). The experiment clearly shows that the mentioned parabolic dispersion becomes fully occupied. Indeed, the two "X"-like features are clearly seen between Γ and M with the top of the above mentioned parabola residing at 90 meV below the Fermi level.
Taking into account that also the energetics of the Γ-centered surface states is not fully accounted for by the calculations, it is not surprising that the surface states appear to be missing at the M-point -the bright spot is not there in the maps and we do not observe any corresponding crossing of features in Fig. 4b and 4d. However, a weak feature at the Fermi level can be seen at the M-point (e.g. left M-point in Fig. 4d). Taking into account that all surface states discussed above are the strongest features in terms of intensity, we attribute this one to the bulk states. Its shape corresponds to the bottom of the parabolic projection seen slightly above the Fermi level in Fig. 1e and thus to the triple Dirac point.
In spite of deviations in energy positions of some surface states the still observed overall agreement implies a remarkable correspondence between theory and experiment for the bulk states. Thus, the 3D Dirac points revealed by the calculations turn out to be experimentally realized. Using the positions of the Dirac points in the calculations one can trace their location in the experimental data. Because of the low k z -resolution, we are not able to see the dispersions crossing in a single point directly, but the pattern of the features nearby allows us to identify the location of their projections (see crosses in Fig. 4b). As mentioned above, the triple Dirac point is located at approximately 10-20 meV below the Fermi level, as indicated by one of the two white crosses in Fig. 4b.

Low energy model
For the construction of the low energy model around the triple Dirac fermion point we used a fine resolution of the Wannier model (see Appendix).
The minimal model for a sixfold fermion of PtBi 2 (space group 205) can be found in Ref. [7]. The Hamiltonian is quadratic in momentum, k = (k x , k y , k z ), always measured relative to the band crossing point: where The symmetries which are relevant for the sixfold fermion are: a 3-fold screw in the [111] direction, denoted S 3 , a 2-fold screw in the x direction [S 2 ], inversion [I], and time-reversal symmetry (TRS, T ). The representation of these operators is where K denotes complex conjugation, σ are Pauli matrices in spin space, the remaining 3 × 3 matrices act in orbital space, and 1 3×3 is the identity matrix.
Using ab-initio data on a momentum grid around the 6-fold fermion, we extracted the parameters E 0,1,2,3 , a, φ a , and b of the Hamiltonian Eq. (1) (see Appendix). Note that φ b can be gauged away by the basis transformation U = diag(1, 1, 1, e iφ b , e iφ b , e iφ b ), and as such it will not affect the spectrum of the low-energy model. For this reason we set φ b = 0 throughout the following.
The magnetic field (here considered to be a Zeeman field) has two main effects. First, it breaks TRS, such that different spin states have in general different energies. Second, it breaks some of the lattice symmetries, depending on its orientation. On the level of the low-energy Hamiltonian, we model a Zeeman field as which is separated into an orbital part that acts on the orbital degrees of freedom, H orb z , and a spin part which is taken to be the usual Zeeman term. For the [111] direction, B x = B y = B z ≡ B, and H orb z is obtained as the most general 3 × 3 matrix which respects the 3-fold screw symmetry [H z , S 3 ] = 0 (see Eq. (3)). Following Ref. [7] we obtain As before, we determine the g factors by comparing the model to the ab-initio data, with details shown in the Appendix.
When the Zeeman field is turned on, the 6-fold crossing splits into a total of 20 Weyl cones. There are 8 Weyl nodes on the [111] axis, and 4 on each of the three principal directions (totaling 12), which are related to each other by the 3-fold screw symmetry. Because the Zeeman term preserves inversion symmetry, the band crossings are k → −k symmetric, with Weyl cones at opposite k having opposite charge.
In the following, we fix B = 3 · 10 −4 and show the Weyl node positions and charges. We find that all 20 Weyl cones are type-II, or over-tilted. Figure 5 shows the spectrum as a function of k x for k x = k y = 0 (left panel), as well as along k x = k y = k z ≡ k d (right panel). In both cases, Weyl points are marked with arrows.
To show that there are no other band crossing points except for the ones mentioned above, we have performed a 3D scan in momentum space, as shown in Fig. 6. In this figure we also indicate the monopole charge of each node, obtained by numerically integrating the Berry curvature over a closed surface surrounding each Weyl cone.

Conclusion
Crossings of energy levels such as Dirac points, Weyl points, and their multi-fold generalizations occur in a variety of bandstructures, but few materials host them in the vicinity of the Fermi level. We have shown that cubic PtBi 2 is one such material. We have performed angle-resolved photoemission spectroscopy to map its bandstructure, finding good agreement with density functional theory calculations. Our results show that PtBi 2 hosts a symmetry protected crossing of six energy bands, 50 meV above E F in DFT and 10-20 meV below E F according to ARPES.
Starting from the ab-initio results, we have used a low-energy effective model for the sixfold fermion, which allowed us to study the effect of a Zeeman field on the band crossing. Remarkably, for a magnetic field pointing in the [111] direction, the band crossing splits into a total of twenty Weyl cones, marking PtBi 2 as an ideal candidate for magnetotransport studies. We hope that our work will motivate further research in this direction, including measurements of the chiral anomaly and studies of non-linear transport effects generated by Berry curvature dipoles [59][60][61].

Acknowledgements
We thank Ulrike Nitzsche for technical assistance.

A Model parameters from ab-initio data
Here we describe the method used to extract the parameters of the low-energy model Eq. (1), given the bandstructure generated in DFT.
The simplest to estimate is E 0 , the energy of the band crossing point, obtained by setting k = 0. Setting k y = k z = 0 instead, the Hamiltonian is diagonal, allowing us to obtain E 1,2,3 , which are just the curvatures of the three parabolic bands (each band is doubly degenerate because of inversion and TRS). Note that the ab-initio data is on a discrete grid, so to get the curvatures we used the k x points closest to the 6-fold crossing (k x = 0.02 in units of 2π/a). This is necessarily an approximation, however, because Eq. (1) is only valid infinitesimally away from k = 0, so that for any finite distance away there will be a (small) contribution from higher-order momentum terms. We find E 0 0.027, E 1 12.192, E 2 13.855, and E 3 20.927, where we have chosen E 1 < E 2 < E 3 without loss of generality. All parabolas have positive curvature.
To extract the remaining parameters, we use the eigenvalues of the system along the diagonal momentum direction k d ≡ k x = k y = k z , which now only depend on b, a cos(φ a ), and a sin(φ a ). As before, each level ε is 2-fold degenerate Again, comparing to the levels closest to the crossing, at k d = 0.02, allows to extract the Hamiltonian parameters. Using the E 0,1,2,3 found before, this is a system of 3 equations with 3 unknowns, but it does not have an exact solution, again because of higher order momentum terms which are nonzero away from k = 0. For this reason, we use instead a numerical minimization routine to find the best a, b, and φ a . The routine minimizes j=1,2,3 |ε k·p j −ε DFT j |, where ε DFT j are the energies obtained by DFT, finding a −24.666, b 14.365, and φ a −0.846. The largest relative error introduced by this process is max j=1,2,3 |ε j − ε DFT j |/ε DFT j = 0.01529, so less than 1.6%. Finally, note that one can always change φ a → φ a + π and a → −a without changing the spectrum, simply because the two always appear as ae ±iφa in the Hamiltonian.
To estimate the g factors entering the Zeeman term (4), we use a similar numerical minimization algorithm. To this end, we added to the non-magnetic ab-initio Hamiltonian a term Σ · B where Σ is the matrix of the spin operator in the DFT basis, which allows to have a Zeeman field that acts correctly on the individual orbitals. We used a diagonal B = 0.01 and extracted the eigenvalues at k = 0: With the E 0 found before, we obtain g 0 0.498, g 1 0.237, and g 2 0.076. The relative errors of energy levels obtained through this numerical minimization (as compared to DFT) are 5.13%, 2.82%, 1.68%, 1.29%, 0.63%, and 0.18%. To get a better estimate, one can also include a linear in B total energy shift of the bands, by adding a term Bg shift 1 6×6 and fitting g shift together with the other g factors. Doing this results in g 0 0.391, g 1 0.348, g 2 0.194, and g shift 0.047. Notice that g shift is significantly smaller than the other terms. The relative errors in this case are also smaller: 1.09%, 0.97%, 0.55%, 0.017%, and the remaining two smaller than 10 −6 %. In the plots shown in Figs. 5 and 6, we have included this overall energy shift.      Gray lines indicate the three principal k x,y,z directions, as well as the diagonal k x = k y = k z direction. Ordering the eigenvalues in increasing energy from 1 to 6, colored patches indicate crossings between levels 2 and 3 (blue), 3 and 4 (red), 4 and 5 (green), as well as between levels 5 and 6 (black). Each colored patch is determined by the set of momenta for which the energy difference of neighboring bands is less than 10 −5 . Six inequivalent Weyl cones are labeled W j , j = 1, . . . , 6. All other crossing points can be inferred by applying inversion or 3-fold rotation. Their monopole charge is Q = +1 for W 1 , W 4 , and W 6 , whereas it takes the value Q = −1 for W 2 , W 3 , and W 5 .