Cobalt-Based Magnetic Weyl Semimetals with High-Thermodynamic Stabilities

Experiments identified Co3Sn2S2 as the first magnetic Weyl semimetal (MWSM). Using first-principles calculation with a global optimization approach, we explore the structural stabilities and topological electronic properties of cobalt (Co-based shandite and alloys, Co3MM-X2 (M/M-=Ge, Sn, Pb, X=S, Se, Te), and identify new stable structures with new Weyl phases. Using a tight-binding model, for the first time, we reveal that the physical origin of the nodal lines of a Co-based shandite structure is the interlayer coupling between Co atoms in different Kagome layers, while the number of Weyl points and their types are mainly governed by the interaction between Co and the metal atoms, Sn, Ge, and Pb. The Co3SnPbS2 alloy exhibits two distinguished topological phases, depending on the relative positions of the Sn and Pb atoms: a three-dimensional quantum anomalous Hall metal, and a MWSM phase with anomalous Hall conductivity (~1290) that is larger than that of Co2Sn2S2. Our work reveals the physical mechanism of the origination of Weyl fermions in Co-based shandite structures and proposes new topological quantum states with high thermal stability.


Introduction
Recent years have seen tremendous development in topological quantum materials (TQMs), including topological insulators (TIs) and semimetals with nontrivial band topology. However, there remain a few challenges to practical applications of TQMs. For example, the challenge in the development of TIs, 1-5 one of the most widely studied classes of topological materials, is to protect their key features against disorders or defects that may destroy their quantized conductance. 6,7 On the other hand, the discovery of new Chern insulators [8][9][10][11] with broken timereversal symmetry is a challenge because of the incompatibility between ferromagnetism and electronic insulation. Topological semimetals (TSMs), defined using a local topological invariant, have become a hot spot in the family of topological materials because of their rich physical properties. One of the most important TSMs is the Weyl semimetal (WSM), whose band dispersion near the Weyl point can be described by the "Weyl equation." 12 Although nonmagnetic WSMs have been thoroughly studied in recent years, [13][14][15][16][17][18] magnetic WSMs (MWSMs) are still rare, despite their great potential as building blocks for next-generation ultra-fast topological spintronics.
MWSMs exhibit many exotic physical properties, such as quantized anomalous Hall effects, 19 large anomalous Hall conductivity (AHC), [20][21][22] and domain wall physics. 23 Since the theoretical discovery of pyrochlore Y2Ir2O7, 24 the first MWSM, several other potential MWSMs have been reported. They include HgCr2Se4, 19 Heusler compounds, [25][26][27][28][29] Fe3Sn2, 30 Sn2Nb2O7, 31 spinel compounds 32 and rare earth materials. 33 Recent experiments confirmed the predicted shandite structure of Co3Sn2S2 as an MWSM by directly observing its Fermi-arc and linear dispersion bulk bands. 34,35 Simultaneously, the magnetic Weyl fermion line was confirmed by directly observing the drumhead surface states in the room-temperature magnet Co2MnGa. 36 However, our understanding of their structural stabilities, microscopic physical mechanism of Weyl points, topological phase transition, and novel properties of cobalt (Co)-based shandites and their alloys is still lacking.
In this letter, we investigate the topological properties of Co3MM'X2 (M/M'=Ge, Sn, Pb and X=S, Se, Te, denoted as Co-MM'-X for simplicity) compounds with high thermodynamic stabilities. Using a global structure search approach based on a particle swarm optimization (PSO) algorithm, we discovered that their structural stabilities as a shandite phase are well described by a structural tolerance factor-defined by the ratio between the atomic radii of the metal and the chalcogen atoms constituting a compound-below which the shandite structure becomes unstable. Shandite becomes a ground state for Co-Sn-S, Co-SnPb-S, Co-Pb-S, and Co-Pb-Se, metastable for Co-Sn-Se and Co-Ge-S, and unstable for Co-Ge-Se and Te-based systems. We determined that Co-SnPb-S alloys can be three-dimensional (3D) quantum anomalous Hall metals (QAHMs) or MWSMs, depending on the relative positions of the Sn and Pb atoms. The MWSM phase had AHC of 1290 Ω -1 cm -1 near the Fermi level, even higher than that of Co3Sn2S2. We established a tight-binding (TB) model that reproduced the number and types of Weyl points for a Co-M-X system, analyzed by first-principles calculations. Our TB model explains that the formation of nodal lines originates from the interlayer coupling between Co atoms in different Kagome layers, and the number and types of Weyl points are mainly governed by the interaction between interlayer M' and Co atoms.  calculations (see Methods section). The crystal structures are color coded and arranged by the energy hierarchy for a given composition; i.e., the one listed first is the most stable configuration. Here, the frequency indicates the number of times the given crystal was identified; thus, it relates the probability of discovering each structure. The energy differences (∆E) between the shandites and the ground state structures for the given compounds are listed in Fig. S1. The compounds Co-Sn-S, Co-Pb-S, Co-Pb-Se, and Co-SnPb-S have the shandite structure (the space group166, red bar in Fig. 2a) as their most stable structure, whereas the shandite structure becomes unstable for Co-Ge-Se and Co-M-Te. For Co-Sn-Se, Co-Ge-S, Co-GeSn-S, and Co-GePb-S, other structures are the most stable ones, with energy differences of less than ~60 meV per formula compared with that of the shandite, showing that the shandite is metastable for those compounds. Note that moderately metastable materials, with energies above the ground state commensurate with the found metastability range (~70 meV/atom, varying by chemistry and composition), are also reasonable candidates for synthesis. 39 We discovered that the ratio between the sizes of the metal atoms and the sizes of the chalcogens is related to the structural stability of a shandite; and we introduced a structure tolerance factor (s), s=r(M)/r(X) as a structure descriptor for shandite, where r(M) and r(X) indicate the radius of a metal and a chalcogen, respectively, and the average atomic radius of two metal atoms is taken as r(M) for metal alloys. As demonstrated in Fig. 2a, only the compounds with large s values ( ≥ 1.5) adopt shandite as a ground state structure. The relationship between the tolerance factor s and the stability of the shandite structure is presented in Fig. 2b. We further analyzed the magnetic stabilities categorized in the ferromagnetic (FM), A-type anti-ferromagnetic (AFM), C-type AFM, and nonmagnetic states (see Fig. S2 of SI), and concluded that all the shandite structures maintain an FM ground state with the spin oriented along the z direction (the energy differences between different magnetic states are summarized in Table. S1). Given the structural and magnetic stabilities, we propose Co-Pb-S and Co-SnPb-S as new magnetic shandite compounds that may be synthesized by experiments.

Crystal structures and magnetic ground states for Co-MM'-X
Interestingly, recently, the alloy system Co3Sn2-xInxS2 has been synthesized 40, 41 by the Bridgeman technique 42 and was identified to maintain a large AHC (1500±300 Ω -1 cm -1 at x=0. 15). Our calculations also confirmed the stability of Co-Sn-S, as it has already been synthesized and identified as an FM Weyl semimetal. 43 In the following discussion, we mainly focus on the Co-MM'-S systems.

Weyl points in Co-M-S (M=Ge, Sn, Pb) systems
The shandite Co-Sn-S has nodal lines located at the kb=kc (kb and ka range from −  where the Weyl points-crossing points between valence bands and conduction bands-are highlighted by the red dots. By plotting the inverses of energy gaps between the conduction and valence bands, we can clearly identify the existence of the Weyl points (Fig. 3c). On the kb=kc mirror plane, Co-Ge-S and Co-Sn-S each have a pair of Weyl points, while Co-Pb-S has three pairs. Note that the Weyl points in each pair are symmetric with respect to the inversion center (L point) in the BZ, thus they maintain opposite chirality. Figure 3d presents 1D band structures, tangential to the nodal lines at Weyl points, that identify type-II Weyl points 14 for Co-Ge-S and Co-Pb-S, and type-I Weyl points for Co-Sn-S. Given the C3z symmetry of the system, there are three pairs of type-II Weyl points in the BZ (another two pairs of Weyl points are located at the ka=kb and ka=kc mirror planes) for Co-Ge-S, nine type-II pairs for Co-Pb-S, and three type-I pairs for Co-Sn-S. Our results for Co-Sn-S are consistent with theory 20,22,44 and experiments 34,35 in the literature. The compound Co-Pb-S had two pairs of Weyl points so close that they almost formed short nodal lines; resolving those Weyl points might be an experimental challenge. We further explored the whole BZ in search of Weyl points located at general points and found additional Weyl points for Co-Pb-S. They were two nonequivalent type-I Weyl points, one located close to the Fermi level (~5 meV higher) and the other one at  Table S2.

3D quantum anomalous Hall phase and large AHC in Co-based alloy systems.
The shandite Co-MM'-S (M/M'=Sn and Pb) has two distinctive configurations, depending on the occupation of Sn and Pb atoms at either at in-plane (M) or interlayer (M') sites of the Kagome plane (see Fig. 1a). Both structures have shandite FM ground states (see Table S1 of SI), and the formation enthalpy of the Co-SnPb-S is lower than that of Co-PbSn-S by ~45 meV per formula unit. Interestingly, these two configurations display very distinctive topological properties: Co-SnPb-S belongs to the 3D quantum anomalous Hall phase, 45,46 whereas Co-PbSn-S is an MWSM with AHC (~1290 Ω #$ cm #$ ) higher than that of Co-Sn-S, as explained below in detail. In a stark contrast to the Co-SnPb-S with no Weyl points, Co-PbSn-S contained type-II Weyl points with SOC. There was a pair of type-II Weyl points (Fig. 4b, for their types and energy levels) in the kb=kc mirror plane and another pair (W " # and W " % ) located at the ka=kb=kc axis (i.e., the z axis in real space) in the BZ with opposite chirality (the specific coordinate of this pair of Weyl points can be seen in Table S2). Note that near this pair of Weyl points, we found that a strict flat band pierces two bands (these conduction and valance to the Fermi level (~15 meV above the Fermi level) (Fig. 4c), which makes it a promising candidate for topological spin-transport devices.
On the other hand, Co-SnPb-S has another distinctive topological property-a 3D QAHM; its band structure can adiabatically evolve into a 3D quantum anomalous Hall insulator 45 (QAHI) without any band crossing. 47,48 To confirm its nontrivial band topology, we followed the criteria proposed by Jin et al. 45 : (1) Table S3), which also satisfied the second condition. Thus, we conclude that Co-SnPb-S is a QAHM with an AHC of ~530 Ω #$ cm #$ , much lower than that of Co-PbSn-S, perhaps because of the lack of Weyl points (Fig. 4b).

Tight binding model for the Co-based shandite structure.
We analyzed the contribution of each atomic species to the emergence of nodal lines through band inversion between the conduction and valence bands. For all the compounds-Co-Sn-S, Co-SnPb-S, and Co-PbSn-S-the dxy and dz2 orbitals of Co atoms were the main contributors to the formation of the nodal lines, with some contributions from the pz orbitals of metals located at the interlayer sites (See Fig. S7), negligible contributions from any metals located at the intralayer sites and the S atoms for the band inversion. We tried to use dxy and dz2 orbitals for each Co atom to construct the TB model. However, we find it can not reproduce the nodal line in the kb=kc mirror plane. We noted that the crystal field which surrounds the Co atom maintains a low symmetry. Thus, the five d orbitals of a Co atom will split and each of them belongs to 1D irreducible representation, with orbital moments quenched (L eff =0) 49 . When the onsite SOC is "turned on," the total angular momentum equals to 1/2 (J eff =1/2). Hence, we try to use a pseudospin Jeff=1/2 (J1/2) orbital 50,51 to construct our TB model.
Our TB Hamiltonian included interactions between J1/2 orbitals on each Co atom, pz orbitals of interlayer metal atoms, and pz orbitals of S atom. For simplicity, we assumed an infinitely large exchange field and limited our discussion in the spin-majority subspace, which reduced to six orbitals in total in the unit cell. Our TB Hamiltonian reads = )* + +*, + -.+/0. + )*#1 + )*#2 + + *3+456 . (1) The first term, )* , describes the hopping interactions between the J1/2 orbitals of Co atoms: )* = ∑ where 47 8 and 97 represent creation and annihilation of an electron on the (J1/2) orbital at i site and j site, respectively; $ and " are the nearest neighbor (NN) and next nearest neighbor (NNN) hopping parameters between the orbitals in the Kagome plane (see Fig. 1), respectively; > is the NN hopping parameter between the interlayer orbitals (see Fig. 1); and < > and ≪ ≫ represent the NN and NNN sites for each, respectively. The effective SOC interaction between the Co atoms is incorporated in the second term, +*, , as the interaction between J1/2 orbitals in the Kagome plane: where +* is the effective SOC hopping strength between the NN J1/2 orbitals and the sign whereis an effective Rashba hopping strength between the NN J1/2 orbitals of Co atoms; where a and represent the J1/2 orbital of Co and the pz orbital of the M atom at the interlayer site, respectively, and C represents the hopping parameter between the nearest orbitals. The term, )*#2 , describes the interaction between Co and S atoms: where is the pz orbital of the S atom and E represents the hopping parameter between the nearest orbitals. The term *3+456 represents the onsite energies of B orbitals of M and S atoms: , where 1 and 2 are the onsite energies of the B of M and B of S atoms, respectively. Here, we set the onsite energies of J1/2 of Co to be zero.
Next, we investigated the role of each term of the Hamiltonian in the emergence of nodal lines, as well as the formation of Weyl points. The first two terms of )* in Eq. (2) describe the interactions between Co atoms in the Kagome lattice plane. We failed to capture nodal lines with those terms solely (see Fig. S8a), but with the addition of the interlayer interaction (the third term in Eq. [2]), we were able to monitor the emergence of the nodal lines centered at the L point of BZ (Fig. S8), even without considering any additional terms from metal atoms. Upon application of the NN SOC interaction between the Co atoms in the Kagome plane (Eq. [3]), the nodal line split into two pairs of Weyl points in the kb=kc mirror plane. Those Weyl points were robust against further application of the Rashba SOC (Eq. [4]). We found it interesting that the interaction between metals and Co atoms (Eq. [5]) resulted in the correct number of Weyl points, as well as their types, as identified by DFT ( Fig. 5a and b): one pair for Co-Ge-S and Co-Sn-S, three pairs for Co-Pb-S, and their types (type-I and type-II). Finally, the last term in the Hamiltonian (Eq. [6]), the interaction between S and Co, slightly modified the shape of the nodal lines and the types of Weyl points in the mirror plane. In summary, the nearest interlayer coupling between J1/2 effective orbitals of Co atoms resulted in the emergence of the nodal lines, whereas the number of Weyl points and their types were mainly determined by the interaction between pz orbitals of M and J1/2 of Co. Additionally, the number of Weyl points and their types were sensitive to the on-site energy of pz orbitals ( 1 ) of the M atoms. Figure   5c presents the phase diagram of a Co-based shandite structure in the parameter space of C and 1 , where the number and type of Weyl points are those located on the ka=kb mirror plane.
At the boundary between W-I and W-II regions, there is a phase transition from type-I to type-II, in which the band crossing is close to a critical tilt. Our DFT results revealed that Co-Sn-Se holds this kind of critical tilt band dispersion at this boundary (see Fig. S5c). These results indicate that our TB captured the essential physical mechanisms that governs the emergence of nodal lines as well as Weyl points in CoMS (M=Ge, Sn, Pb) systems. It is noted that although our TB model does consider the role of intralayer metal site, we can describe the QAHM state which the alloy system Co-SnPb-S maintains. However, for Co-PbSn-S alloy, based on our TB model, with changing the parameters (t4 and 1 ), we did not find the pair of Weyl points which locate at the z axis ( " # and " % , see Fig. 4(b)). We guess this may be attributed to two reasons.
One is that Weyl points emerge occasionally with fine-tuned parameters. Hence, it is challengeable to find them if the proper parameters are not given. The phase diagram (Fig. 5(c does not consider the interaction between Pb (intralayer) and other atoms. This simplification may also limit the finding of the pair of Weyl points ( " # and " % in Fig. 4(b)) for Co-PbSn-S system.
In summary, we systematically investigated the stability and electronic properties of shandite Co-based materials. For structural stability, we found a good "descriptor" (i.e., r[M]/r]X]) that can identify not only the stability of the shandite crystal structure but also the stability of the magnetic states. For electronic properties, we found that the behavior of Weyl points was strongly dependent on chemical composition. In particular, the Weyl semimetal Co-PbSn-S exhibited a large AHC very close to the Fermi level and is a good candidate as a spintronic transport device. Finally, with an effective TB model, we revealed that the band inversion nodal line in the Co-based Kagome lattice was derived from the interlayer coupling between Co atoms and that the behavior of Weyl points was mainly determined by the interaction between Co atoms and the interlayer M atoms. Our work gives a clear physical picture of the origination of the Weyl fermion in the Co-based Kagome lattice and provides guidance for synthesizing new stable TQMs.

DFT calculations
Electronic structures were determined by DFT calculations as implemented in the Vienna ab initio simulation package. 53 The projector augmented wave 54, 55 was used to treat the ionelectron interactions, and a generalized gradient approximation of the Perdew-Burke-Ernzerhof 56 functional was applied to consider the exchange-correlation potential. The cutoff energy and k-mesh were 500 eV and 8×8×8, respectively. All the structures were optimized using a conjugate gradient approach with the atomic forces smaller than 0.01 eV/Å.

PSO simulations
Crystal structure predictions for each chemical composition were carried out by employing the PSO algorithm implemented in the software Crystal Structure Analysis by Particle Swarm Optimization (CALYPSO). 38 For a single prediction, each generation consisted of 30 structures, and structural evolution proceeded to at least the tenth generation. For each generation, 60 percent of the structures were constructed through the PSO algorithm, and the remaining 40 percent were randomly produced to prevent premature convergence of structural prediction.
During the structure search process, structural relaxations and total energy calculations were performed using first-principles calculations.

Anomalous Hall conductivity
First, for T 3 MM'X 2 , we use 3d orbitals of T, 5p orbitals of M(M'), and 3p orbitals of X to construct the TB Hamiltonian as implemented in the Wannier90 package. 57 Then, based on this TB Hamiltonian, we used the WannierTools 58 to calculate the AHC for the T 3 MM'X 2 system.
The AHC can be obtained as the sum of Berry curvatures for the occupied bands 59 :

IJ
Here, is the Berry curvatures, n is the index of the occupied bands, and f is the Fermi-Dirac distribution function. We used a 150×150×150 k-points grid to calculate the value of AHC.

Effective TB model
The numerical TB model was implemented in the Pythtb 60 package, which uses the orthogonal TB approach. The Bloch-like basis function can be constructed as follows: R and i are the site and orbital index, respectively; ti is the position of orbital i in the home unit cell; is the TB basis orbital in cell R and should satisfy the following condition: which is the so-called orthogonal TB approach. The eigenstates can be expanded as Then, the Hamiltonian at k point can be written as Here, the phase factor 4 •S % $ # " T is included, which is called convention I. 60 By solving the matrix equation &• 3 = 3 3 , we can get the eigenvalue and eigenstate at each k point and thus get the whole bands in the BZ.

Data availability
The data that support the findings of this study are available from the corresponding authors upon reasonable request.