Predicted Realization of Cubic Dirac Fermion in Quasi-One-Dimensional Transition-Metal MonoChalcogenides

We show that the previously predicted Fermion particle that has no analogue in the standard model of particle theory - the cubically dispersed Dirac semimetal (CDSM) - is realized in a specific, stable solid state system that has been made years ago, but was not appreciated to host such a unique Fermion, composed of six Weyl Fermions, 3 with left-handed and 3 with right-handed chirality. We identified the crystal symmetry constraints and found the space group P63/m as one of the two that can support a CDSM, of which the characteristic band crossing has linear dispersion along the principle axis but cubic dispersion in the plane perpendicular to it. We then conducted a material search using density functional theory identifying a group of quasi-one-dimensional molybdenum mono-chalcogenide compounds A(MoX)3 (A = Na, K, Rb, In, Tl, X = S, Se, Te) as ideal CDSM candidates. Studying the stability of the A(MoX)3 family reveals a few candidates such as Rb(MoTe)3 and Tl(MoTe)3 that are predicted to be resilient to Peierls distortion, thus retaining the metallic character. The combination of one-dimensionality and metallic nature in this family provides a platform for unusual optical signature - polarization dependent metallic vs insulating response.


The topological significance of band crossings in solids:
The crossing on energy bands in complex materials showing dense manifold of states is a ubiquitous effect routinely reported in the past ~ 50 years in countless publications, a visual effect often referred to as "band spaghetti". Such crossings have been known for a long time to result from specific space group symmetries [1,2], involving various band dispersion E(k) with powers k n and to be characterized by a band degeneracy g(k) at the crossing wavevector k.
With the renewed interest in materials with strong spin-orbit coupling (SOC) it became recently clear that such band crossing points could carry interesting information on the topological behavior of the system, leading to specific behaviors of surface/edge states in reduced dimensions [3][4][5]. In SOC system respecting time-reversal symmetry, the degeneracy of multi-band crossings in momentum space could be g(k) = 2, 3, 4, 6, and 8, and the dispersion power k n with n = 1, 2, 3 at the crossing point could be linear, quadratic or cubic, respectively. Crossings between conduction and valence band with linear dispersion in at least one k-space direction and g = 2-fold degeneracy corresponds to Weyl semimetal (WSM) [3,[6][7][8][9][10] and g = 4-fold degeneracy corresponds to Dirac semimetal (DSM) [4,[11][12][13][14][15][16]. The remaining degeneracies at band crossing points g = 3, 6, and 8 correspond to quasi-particles without analogous states in the standard model of particle physics [5,[17][18][19]. The latter respects Poincare symmetry and has but three is linear, whereas the dispersion along the a-b plane could be either linear (n = 1), quadratic (n = 2), or cubic (n = 3). A Weyl point with in-plane dispersion power n carries a Chern number n or -n, corresponding to a degeneracy of n conventional Weyl Fermions (left-or right-handed) all with the same chirality [7]. In contrast, a Dirac point with specific n, has zero net Chern number, and is 2n-fold degenerate Weyl Fermion with half left-handed and half right-handed chirality [15,20,21]. Such Weyl and Dirac Fermions with high-order dispersions, caused by crystalline symmetry in solids, also do not have counterparts in high energy physics. While the surface states of DSM are not topologically protected (unlike those of TI and WSM), such quadratic (n = 2) and cubic Dirac Fermions (n = 3), especially the latter, are interesting relative to the conventional n = 1 Dirac Fermions: Their distinguishing features (as discussed in the concluding remarks) include creation of special WSM with multiple Fermi-arcs; characteristic quantum transport signatures; quantum criticality and phase transitions.

Search for materials that realize cubic Dirac Fermions:
In addition to the efforts to define and classify such specific "new Fermion" band crossings induced by crystalline symmetries and topology in condensed matter systems [5], an important challenge is to systematically identify material realization of such unusual Fermions. We summarize all the different types of Fermions in solid state physics, classified by the degree of band degeneracy (g) and the highest power of band dispersion (n), with example materials in Appendix A. Whereas quadratic dispersion has been predicted to exist in SrSi2 [22] (g = 2), band-inverted -Sn (g = 4), and PdSb2 (g = 6) [5], the cubic dispersion, expected to exist both in Weyl [7] and Dirac semimetals [15], has not been realized as yet in any material candidates. In part, the difficulty to find such cubically-dispersed Dirac semimetals (CDSM) is related to the multitude of nontrivial requirements one needs to impose on such a material search, including appropriate crystal symmetry, angular momentum and electron filling. Finding such materials by accidental discovery or simple trial-and-error would thus be unlikely.
Here we establish understanding-based design principles for CDSM and use these to deliberately screen the candidates that satisfy such conditions by exhaustively looking through all 230 space groups in 3D. We find that only materials in two space groups P63/m (No. 176) and P6/mcc (No. 192) have the potential to host cubic Dirac Fermions.
This narrowing down of the possibilities is then followed up by a (yet non-exhaustive) material search of compounds belonging to these two CDSM-hosting space groups using density functional theory (DFT, see Appendix B for computational details). We identify a group of molybdenum transition-metal chalcogenide compounds A I (MoX VI )3 (A I = Na, K, Rb, In, Tl; X VI = S, Se, Te) with space group P63/m as ideal candidates. The structure of this type of compound ( Fig. 1) is basically quasi-1D chains (Mo3X3) 1-

Crystal structure and chemical bonding of the A(MoX)3 family:
The quasi-1D structures A(MoX)3 system is derived from the general family, known as "Chevrel clusters" [37], that has been synthesized with extended MoX clusters Mo6X8. The Mo6X8 unit can be viewed as a Mo6 octahedral surrounded by eight X chalcogen atoms, or two In addition to the odd number of electron filling, we have identified this family of compound as a CDSM (Table I)  (a more detailed derivation could be found in Ref. [7,15]). Therefore, the leading order of k-dispersion at the DP can thus be written as p-q mod n.
In A(MoX)3 compounds, the little group of the A point is C6h with a C6 rotation axis.
where ± = ± with the origin at the DP, and , , are independent real coefficients. Such a block-diagonal matrix can be decomposed into two Weyl Hamiltonians, with each being characterized by a topological invariant, i.e., a Chern number, which is defined as the number of monopoles of Berry curvature of a closed 2D surface enclosing the Weyl nodes. In A(MoX)3, the DP at A is composed of two opposite cubic Weyl Fermions [7] carrying Chern numbers +3 and -3 joining without annihilation.
In other words, such as a DP can be viewed as a being composed of 6 conventional Weyl Fermions with 3 having left-handed and 3 having right-handed chirality. It is known in high energy physics that the Weyl Fermions with the same chirality cannot be degenerate, so such a 6-fold DP is indeed a "new Fermion", caused by symmetries present specifically in crystals.
Next we consider the L point, which has the little group C2h with a C2 rotation axis.
Thus, the band eigenvalue of rotation operator is ± , and {p, q} = {0, 1}. Therefore, the DP at L has linear dispersion along kx and ky direction, as confirmed by DFT calculation shown in Fig. 2(b). The type of linear DP at L in A(MoX)3 is identical to that of hypothetical β-cristobalite BiO2 [11] and distorted spinel BaZnSiO4 [40], but is more experimental accessible because of its stability (will discuss later) and the success history of synthetization [23]. In addition, the DP at L in Tl(MoTe)3 is only 7 meV below the Fermi level, which hopefully could be observed by ARPES.

Dirac Fermions within the conduction bands:
The crystal symmetry of A (MoX)3 compounds can also host another type of Dirac Fermions, which is induced by the band inversion between conduction and valence bands with an accidental band crossing inside the BZ. Examples include Na3Bi, Cd3As2 that were verified by ARPES measurement [13,14], and ternary honeycomb materials such as BaYBi (Y = Au, Ag and Cu) [41] and metastable allotropes of Ge and Sn [42] predicted by first-principles calculation.
However, because the conduction and valence band of A(MoX)3 family only meet at A and L points, such DPs are actually band crossing of two conduction bands or two valence bands at points along Γ-A direction. In Tl(MoTe)3 there are two inequivalent DPs having the energy ~ 300 meV above the Fermi level, as marked by the purple circles in Fig. 2a. Interestingly, although they located close to each other in terms of both momentum and energy, the dispersion properties of the two DPs are quite different. In contrast to the band crossings at A and L that are protected by non-symmophic symmetry, the DPs at originate from the inversion of bands with different parities.
Thus, the inversion operator takes the form P = ± , and the low-energy Hamiltonian of the DP deviates from Eq. (1) accordingly [38]. originates from band crossing due to band inversion and happens in pairs around a TRI kpoint (see Fig. 1f); while DP at L is the touching point between conduction and valence band at the boundary of BZ due to non-symmophic symmetry; second, as well as A and 1, DP at 2 show isotropic dispersion between kx and ky direction because the DPs along z axis feels SO(2) symmetry at small in-plane k, while DP at L has anisotropic linear dispersions along kx and ky direction, as shown in Fig. 2b. We summarize the physical characters of the different DPs in Tl(MoTe)3 in Table I.

Peierls distortion and stable structures in quasi-1D A(MoX)3 compounds: As noted
above, the group of materials A(MoX)3 with reported P63/m structure has been made and is laboratory stable. However, the structure determination done so far [23,43] did not shed light on the possibility of possible Peirles distortion. It is natural to expect that the ideal 1D metallic structures are unstable against Peierls distortion, and could have in the undistorted structure either soft phonon modes or higher energy than the distorted phase.
Such distortion could destroy the crystal symmetries responsible for the DP and thus open a gap. To predict realistic quasi-1D DSM candidates we thus studied systematically by DFT the phonon spectra (shown in Fig. 6) and thermodynamic stability (shown in Table III  and Se, including the four compounds having soft phonon modes, are highly unstable in undistorted structure. As a result, the ground states of these materials experience Peierls transition to lower the total energy by 5.5-6.9 meV/f.u., and thus become semiconducting. Mo triangles tend to become pairs by moving towards to each other, forming alternative short-long-short-long bonding with each other. Interestingly, the 3 Se atoms within the same plane of each Mo triangle tend to move oppositely and thus forms buckled in-plane structure. We applied such distortion mode to the undistorted structure and after relaxation we find that such Peierls distortion indeed eliminates the negative phonon modes (see Fig. 4c) as well as lowers the total energy by 6.8 meV per f.u.. The distorted structure (see Fig. 1d) has a reduced symmetry with a space group of P3 ̅ (No. 147), in which the screw axis symmetry is no longer preserved. As a result, a band gap is opened at the A and L points. Fig. 4d shows that in distorted K(MoSe)3 there is a 280 meV band gap throughout kz = π plane, indicating that the relative small change in band length between Mo triangles (0.03 Å) induces remarkable effect in electronic structure.
Indeed, the existence of Peierls distortion is basically the competition between band eigenvalues and elastic energy. If the gain of occupied band eigenvalues induced by creating a gap is less than the cost of elastic energy by modulating the atomic positions, Peierls distortion will not happen. The subtlety of whether a quasi-1D system would experience Peierls transition is closely related how "1D" the system is. To demonstrate this, we investigate the relationship between the stability of undistorted A(MoX)3 compounds and their character of one-dimensionality. Since the linear Dirac bands along extend along the chain direction through the inter-triangle Mo-Mo bonding, leading to large dispersion. On the other hand, the energy bands within a-b plane is rather narrow, implying weak in-plane hopping. Especially, the flatness of the bands within kz = π plane reflects the instability due to Fermi nesting. Therefore, we define a parameter = The case of Peierls semiconductor, e.g., K(MoSe)3, is slightly different in that the zdirection behaves like a small-gap semiconductor, evidenced by nonzero energy onset of the absorption peak shown in Fig. 5b. Future measurement of the polarization dependence will serve to disentangle the basically different optical properties of 1Dmetal and the semiconductor.

Discussion
low temperature, while in high temperature the high-symmetry metallic state is usually more stable because the gain of band eigenvalues is reduced by the thermal excitation of electrons across the band gap. Given that the structural determination of single crystal assigns A(MoX)3 to the high-symmetry space group P63/m at high temperature, the message from Fig. 3 could provide some implications on the metal-insulator transition of this system that can be compared with the transport measurements [27,28] and theoretical explanations using electron-phonon coupling [25]. We conclude that the compounds in the "Peierls semiconductor" region should undergo a phase transition from metal to semiconductor when the temperature goes below a critical value, while at low temperature the compounds in the "Coexistence" region could appear both phases based on the growth condition. The phase diagram with three regions shown in Fig. 3  at least two types of sample with one metallic down to low temperature and another having phase transition a critical temperature [27]. Our results on thermodynamically stability are also consistent with the recent calculation using the electron-phonon coupling, concluding that the main mechanism of the metal-insulator transition is the dynamic charge density wave that corresponds to Peierls-type displacement [25].

Experimental accessibility: The A(MoX)3 compounds have been synthesized over 30
years ago [23] at 1000-1200 ºC in sealed molybdenum crucibles under low argon pressure. The single crystals were needle-like growing along the c axis, indicating the 1D character. Most of this family of compounds were found to be stable in air. It is noticeable that Na(MoSe)3, In(MoSe)3 and Tl(MoSe)3 were synthesized and reported to be a quasi-1D superconductor with relatively high critical fields [27,44,45]. Some compounds from this family were shown to have metal-insulator transition below a critical temperature [25,27,28]. According to our DFT calculations, the compounds that do not have Peierls distortion and thus tend to stay metallic are Tl(MoTe)3 and Rb(MoTe)3. These tellurides are also predicted to have stronger dispersion than selenides and sulfides (see Fig. 7), providing a better chance to resolve the bands by ARPES.
Another challenge for detecting the cubic Dirac fermion is the possibility that the A cation in A(MoX)3 will show some off-stoichiometry, e.g., deficiency which is a natural consequence of the high temperatures required in crystal growth. Indeed, single crystal diffraction showed that while the structures of the quasi-1D MoX chain is nearly perfect, the A-site shows up to 15% under stoichiometric [45]. One might be able to re-plenish the A-site deficiency on the surface by deposition of the cation ultra-high vacuum evaporation so as to make the DP measureable by ARPES.  [46]. Then the surface behavior of a DSM will follow the direction where a small perturbation will lead the system to. For DSMs with a pair of DPs located away from the TRI k-points and P = ± (parity inversion), such as Na3Bi and Cd3As2, a small perturbation can open the gap while preserve the band inversion, leading to a topological insulator (TI). Thus, the surface states are robust as a closed Fermi contour.

Why are cubic Dirac Fermions interesting relative to conventional Dirac Fermions?
DSMs with DP located at the TRI k-points, such as BiO2 (in a hypothetical SiO2 structure platform to realize other topological phases by the symmetry tuning [11,15,47]. The main points of interest are: (i) By breaking time-reversal symmetry in DSM, e.g., via introduction of magnetic ions, such systems are known to transform to WSMs. More specifically, for CDSM this was predicted theoretically [15] to lead to WSM with the unique occurrence of three Fermi arcs connecting the surface projection of the Weyl node and its antinode. This is the largest number of pairs of Weyl nodes that can theoretically be accommodated, leading to an enhanced conductivity step for the quantum anomalous Hall effect. (ii) The dispersion power n (= numbers of monopole charges) present in symmetry-broken DSM was predicted to produce n-dependent quantum interference effects [48], leading to dispersion-dependent quantum transport phenomena. For conventional DSM/WSM (n = 1), a destructive quantum interference results in a weak antilocalization correction proportional to −√ in the weak field limit. Such negative longitudinal magnetoresistance, also known as chiral anomaly [49,50], has been confirmed by various transport measurements [51][52][53][54][55]. In contrast, for DSM/WSM with n = 2, it has been predicted that a weak localization correction proportional to +√ applies on the magnetoconductivity [48], the Coulomb interactions along the in-plane directions are screened with a faster decaying than that along the rotation axis (r -1 ). Recently, it was predicted that WSM with n = 3, in the presence of short-range interactions, can easily undergo a continuous quantum phase transition into either a translational symmetry breaking axion insulator or a rotational symmetry breaking nematic state [59]. Furthermore, the nonlinear dispersion and the 1D nature of a condensed matter system would cause a breakdown of the interacting Fermi liquid theory for electron behavior, leading to Luttinger liquid instead.      Table II shows the examples of materials that host different types of Fermions classified by the degree of degeneracy (g) and the highest power of band dispersion (n).

A. Classification of different types of Fermions in solid state physics
Here we consider three-dimensional crystals with spin-orbit coupling, respecting timereversal symmetry. We consider single-point degeneracy in the k-space, so the materials with line node is not included. Some cases, e.g., g = 8 and n = 3, are forbidden because of the restriction of crystal symmetries. Some cases are predicted to exist but no material realization yet, marked by "?" in the

B. First-principles calculations
All The calculations including total energy, electronic structure and phonon dispersion were performed by density functional theory (DFT) where the geometrical and total energies are calculated by the projector-augmented wave (PAW) pseudopotential [60] and the exchange correlation is described by the generalized gradient approximation of g n Perdew, Burke and Ernzerhof (PBE) [61] as implemented in the Vienna ab initio package (VASP) [62]. The plane wave energy cutoff is set to 500 eV, and the electronic energy minimization was performed with a tolerance of 10 −5 eV. Spin-orbit coupling is taken into account self-consistently throughout the electronic structure calculations. The atomic projection on band structure is calculated by projecting the wave functions with plan wave expansion on the orbital basis (spherical harmonics) of each atomic site. The phonon spectra were identified by PHONOPY package [63], in which the force constants are calculated in the framework of density-functional perturbation theory (DFPT) [64].
Phonon calculations were performed within a 1×1×4 supercell (56 atoms). To evaluate the anisotropic optical properties of quasi-1D A(MoX)3 compounds, we calculate the frequency ( ) dependent dielectric function ( ) based on DFT. Then the optical absorption coefficient is evaluated from the dielectric function: where is the speed of light and [ ( )] is the real part of ( ).
In Table III,

C. Space groups that host cubic Dirac fermions
Here we use symmetry analysis to show that how non-symmophic symmetry ensures In a spin-orbit system, the anti-unitary operator T behaves as T 2 = -1, leading to Kramers degeneracy. Together with inversion symmetry, it turns out that all the energy bands are two-fold degenerate with the two components related each other by PT, i.e.,