Comprehensive scan for nonmagnetic Weyl semimetals with nonlinear optical response

As the development of topological band theory, comprehensive databases about time reversal and crystalline symmetries protected nonmagnetic topological materials were developed via first-principles calculations recently. However, owing to the low symmetry requirement of Weyl points, the symmetry-based topological indicator cannot be applied to Weyl semimetals (WSMs). Hitherto, the WSMs with Weyl points in arbitrary positions are still absent in the well-known databases. In this work, we develop an efficient algorithm to search for Weyl points automatically and establish a database of nonmagnetic WSMs with Weyl points near Fermi level based on the total experimental noncentrosymmetric crystal structures in the Inorganic Crystal Structure Database (ICSD). Totally 46 Weyl semimetals were discovered to have nearly clean Fermi surface and Weyl points near Fermi level within 300 meV, and 9 of them are chiral structures which may host the quantized circular photogalvanic effect. In addition, the nonlinear optical response is studied and giant shift current is explored in the end. Besides nonmagnetic WSMs, our powerful tools can also be used in the discovery of magnetic topological materials.

Generally speaking, topological semimetals are the systems whose conduction and valence bands cross with each other in the Brillouin zone and cannot be removed by perturbations of Hamiltonian without breaking specific symmetries. Topological WSM is such a material with double degeneracy Weyl points. Compared with other topological states, the WSM is extraordinary since the topology is not protected by any crystalline symmetry. It can only appear in the system whose spin degeneracy is removed by breaking time reversal or spacial inversion symmetry. In this case, the bands that create the Weyl points take the 2 × 2 Weyl equation form, H = ±v F k · σ, which can be viewed as half of the massless Dirac equation with defined chirality in three-dimensional space [33]. Because a two-bands model Hamiltonian can be described by a linear combination of three Pauli matrices generally, the linear band crossing can be obtained by tuning the parameters without considering any additional conditions. Therefore, besides time reversal or spacial inversion symmetry, the lattice periodic boundary condition is the unique requirement for Weyl points. Accordingly, the identification of topology based on the characteristics of high-symmetry points does not work anymore for WSMs, and that is why the WSMs are still absent in the comprehensive databases [30][31][32].
In this work, we only focus on the experimental structures to complete the databases by including nonmagnetic WSMs. The total information about their topological properties is listed in the Supplementary Information. At the same time, two effective tools used for generating Wannier function automatically and searching Weyl points are developed in Fig.   1. The accuracy of them is essential for finding materials with observable properties of Weyl points, Some of them have already been theoretically predicted and experimentally confirmed, such as the family of TaAs [23][24][25][26]. Besides type-I and type-II WSMs, 9 chiral WSMs with inversion and mirror symmetries breaking are also listed in our database. For chiral WSMs, the truly quantized circular photogalvanic effect (CPGE) [34][35][36][37] could be induced by the energy difference of Weyl points with opposite chirality. Moreover, we explored the bulk photovoltaic effect (BPVE) [38] of all these noncentrosymmetric WSMs.
A giant BPVE which reach up to 700µA/V 2 is predicted, and some of them can even extend to the regime of visible light.

RESULTS
Our automatic workflow for searching nonmagnetic WSMs is shown in Fig. 1. It includes two crucial algorithms, the Wannier function generation and Weyl points search.
Constrained by the accuracy of the density functional theory (DFT), we don't consider the systems with strongly correlated f-elections, namely Lanthanide and Actinide series except element La, while beginning our work with the experimental noncentrosymmetric crystal structures in the Inorganic Crystal Structure Database (ICSD) [39]. Here we remain some systems with 3d-elections, which may also be strongly correlated systems, to make sure that we don't omit some potential candidates. After removing alloys and repeated structures, totally 8896 compounds left. Among them, all the systems with the total magnetic moments exceeding 0.05 µ B per unit cell and magnetic elements V, Mn, Fe, Co and Ni are regarded as magnetic compounds. Furthermore, the gap and the charge carrier density (N ) of the nonmagnetic phases are checked. Because we are only interested in WSMs with observable topological Weyl-points-related properties, the insulators and good metals (with N > 2 × 10 22 /cm 3 , which is obtained by the comparison with Co 2 VGa [40]) will not be considered in the further calculations.
After all the screening processes above, we obtain the list of the nonmagnetic noncentrosymmetric experimental crystal structures who may possess Weyl points near Fermi level.
For these compounds, the atomic-orbital-like Wannier functions are constructed by an effective algorithm automatically as shown in Fig. 1b. The crucial problem to sole for this step is the setting of the fitting windows with appropriate projection orbital basis sets. In the nonmagnetic first-principles calculations, we select s orbitals of alkaline-earth metals, d orbitals of transition metals and s + p orbitals of p-block elements as the minimum basis sets. Thus, we can label the maximum weight value among the selected band as r. The weight represent the ratio of selected basis projection. If the weight r is smaller than 0.8, the minimum basis sets will be insufficient to produce the DFT band structure; hence, we switch the bases to the maximum basis sets which include s + p + d orbitals of all the elements.
With the selected basis sets, the inner projection energy window can be easily fixed as 90% of the maximum weight value r. Subsequently, we loop over the outer projection energy window and check the mean absolute error between the DFT and Wannier band structures around the Fermi level until the error is smaller than 0.02 eV. Meanwhile, it is noteworthy that a good projection can be always guaranteed because of the large band gap between valence and core states.
Using the high-symmetry tight-binding (TB) model Hamiltonians obtained in the construction of Wannier functions, the Weyl points can be scanned. Since the existence of Weyl points in noncentrosymmetric nonmagnetic system only requires the translational symmetry, the fact that the Weyl points would locate at generic k-points without any special symmetry renders the search difficult. The typical way to search for the Weyl points use the band gap as a criterion. However, this method is much easy to reach the local minimum value, thereby losing some Weyl points. In this work, we use the Berry phase of the parallelepiped corresponding to the k-mesh of the reciprocal space as the criterion. Initially, based on the TB model Hamiltonian, we generate a relatively sparse 51 × 51 × 51 k-mesh along the three reciprocal vectors to calculate the Berry phase of each parallelepiped. If the value is larger than 0.01 × 2π, Weyl points might be present inside the parallelepiped. For this type of parallelepiped, we increase the k-mesh finer and repeat calculating Berry phase. Eventually, the Berry phase will increase in some fine parallelepiped and finally reach a value n × 2π, where n is an integer, corresponding to a single Weyl point. In our calculations, the parallelepiped is refined to around 1/1600 3 of the reciprocal lattice to obtain the Berry phase with n 0.99.
After searching the Weyl points, we have a raw list of the WSMs with all the information of Weyl points. It may also have some magnetic materials such as antiferromagnets. Thus, we process a double check here to modify the raw list, and the list of nonmagnetic WSMs are obtained.
Despite having as many as 8896 noncentrosymmetric candidates, we only find 46 nonmagnetic WSMs with Weyl points near the Fermi level within 300 meV, including 9 chiral crystal structures. This result express that the good nonmagnetic WSMs are difficult to find. The details of every WSMs that we find by our algorithm are given in Supplementary   Table 1. Until now, the Weyl points have been classified into type-I and type-II according to the shape of the linear band crossings. The Fermi surface would shrink to a singularity point for the type-I Weyl point, while the type-II Weyl point presents a Fermi surface with a linear crossing which is composed of electron and hole [41]. On the other hand, another main difference is the precondition of chiral anomaly. It can exist along any direction in the type-I WSM as long as electric and magnetic fields are not perpendicular to each other, but only exist when the electric and magnetic fields point to some selected direction in the type-II WSM.  Fig. 2e. The Fermi arc length can be ∼ 15% of the reciprocal lattice vector, which makes it easy to detect by surface detection methods, such as angle-resolved photoemission spectroscopy (ARPES) and scanning tunnelling microscopy (STM).
The other example Ag 2 S is a hybrid WSM [42] containing both two kinds of Weyl points.
It also has an orthorhombic crystal structure with space group P 2 1 2 1 2 1 (No. 19), as shown in One point we need to emphasized is that Ag 2 S does not present any mirror symmetry.
Although the mirror symmetry is not necessary for the existence of Weyl points, the chirality of Weyl points must be odd with mirror operations. Therefore, the mirror symmetry forces one pair of Weyl points with opposite chirality locate at the same energy. In another words, WSM without mirror symmetry exhibits a net chirality if the Fermi energy is located between one pair of Weyl points with opposite chirality. One of the most important properties for mirror-symmetry-broken WSMs is the quantized CPGE, which is the only possible quantized signal induced by the monopole Weyl point so far. This can also be applied to Ag 2 S with W1 ∼ 23 meV above the charge neutral point and W2 ∼ 28 meV below the charge neutral point, as illustrated in Fig. 3d. Because the charge carrier density is only around 5 × 10 20 /cm 3 and the contribution of the trivial Fermi surfaces can be neglected, Ag 2 S may be an ideal candidate to observe the quantized CPGE.
The chiral WSM Ag 2 S is further confirmed by the chiral surface Fermi arcs. Fixing the energy at W1 and W2, we can obtain the extremely long chiral Fermi arcs in (001) surface in Fig. 3f and 3g, respectively. Owing to the large spin splitting, the spin texture should be easily observed by both ARPES and STM. Since two W2 Weyl points overlap with each other when they are projected into (001) surface, two Fermi arcs start from same projected Weyl point. These extremely clear topological properties offer a good platform for both surface detection and bulk nonlinear optical transports.

Shift Current Analysis
Very recently, a giant shift current was observed experimentally in WSM TaAs [43], which suggests a topological effect in BPVE and the practical applications of WSMs in optical detectors and solar energy conversion. Inspired by this work, we screened the nonlinear optical response for all the nonmagnetic WSMs in our database. Half of them can host a large shift current above 100 µA/V 2 with photon energy larger than 0.1 eV (Supplementary Table 2 and Supplementary Figure 1 Through perturbative analysis, the expression of shift current [44] is written as where i, j, k = x, y, z. The second order optical conductivity tensor can be computed as where n and m are band indexes, f nm = f n − f m is the Fermi-Dirac distribution difference, In the experimental measurement of WSM TaAs (ref. [43]), an apparent measured value σ xxz = 34 ± 3.7µAV −2 with photon energy ω = 0.117eV is obtained. Our calculation gives a well-fitting theoretical value σ xxz = 27µAV −2 with the same photon energy. Furthermore, ref. [43] added the reflectance of TaAs at 10.6µm to the apparent measured value manually and the giant shift current σ xxz = 154 ± 17µAV −2 is obtained in the end. At the same time, σ zzz has the biggest value for TaAs in both our calculation and ref. [43]. Beside the largest shift current candidate BiTeI, we also find another two kind of candidates. One is multifunctional candidates combining topology, strong shift current, and superconductivity. As indicated in Fig. 4b and 4c, the shift current tensor can reach up to σ zzz = 564 µA/V 2 and σ zzz = 250 µA/V 2 for the new found WSMs LaIrP and LaRhP, respectively. In addition, since both LaIrP and LaRhP are superconductors [45] with transition temperatures T c = 5.3 K and 2.5 K, these two materials will provide good platform for the study of the interplay among topology, superconductivity and nonlinear optical responses. The other one is the candidates with large photon energy. For all these WSMs with large shift current in Supplementary Table 2 and Supplementary Figure 1-3, we find that the peak value of σ ijk locates at an energy below 1 eV, except SrHgSn, CaBiAg and CaCdPb. Fig. 4d shows the shift current tensor of SrHgSn with space group P 6 3 mc (No. 186) as an example. One can find a peak value σ zzz = 165 µA/V 2 at photon energy ∼ 1.8eV , which is already located at the range of red light.
The shift current can be considered as a joint effect of the frequency-dependent density of states (DOS) and the wave function. Taking SrHgSn as an example, we can try to understand the microscopic mechanism in these two aspects. From the DOS of SrHgSn shown in Fig.   5b, one can find two peaks locating bellow and above the charge neutral point, respectively.
The energy difference is around 1.8 eV, which is very close to the two peaks of σ zzz in Fig.   4d. σ zzz ( k; 0, ω, −ω) is also calculated for each k-point with the photon energy 1.8 eV and 1.95 eV near the two peaks. Thus we obtained that the positive contribution is mainly from the special plane k z ∼ 0.31 2π/c at photon energy around 1.8 eV, while k z ∼ 0.13 2π/c plane at photon energy around 1.95 eV has a negative contribution. The distribution of σ zzz ( k; 0, ω, −ω) at photon energy 1.8 eV is dominated by twelve positive blue lines caused by the excitation from band n − 1 to bands n + 3 and n + 4, as shown in Fig. 5c and 5f. The distribution of σ zzz ( k; 0, ω, −ω) changes sharply when photon energy increasing from 1.8 eV to 1.95 eV. At photon energy 1.95 eV, σ zzz ( k; 0, ω, −ω) is dominated by two negative red rings as shown in Fig. 5d and 5e. These two red rings are derived by the light excitation from two occupied bands n and n − 1 to the non-occupied bands n + 7, n + 8, n + 9 and n + 10, as presented in Fig. 5g and 5h. Therefore, the two peaks of σ zzz (0, ω, −ω) in Fig. 4d originate from large DOS located around E F − 1.2eV and E F + 0.6eV , and the sign changes because of the different band excitation at different k-points. One point need to emphasize is that the strong shift current at high photon energy is not caused by the Weyl points here.
In summary, through high-throughput calculations, we have complemented the nonmag- We take these crystal structures as the input of the full-potential local-orbital (FPLO) minimum-basis code [46]. Here we choose the generalized gradient approximation of the Perdew-Burke-Ernzerhof functional [47] as the exchange-correlation potential. A dense 12 × 12 × 12 k-mesh is used in the static nonmagnetic calculation. SOC is also included self-consistently during the workflow. In order to get the surface states of WSMs, Green's function method [48,49] is applied with a half-infinite boundary condition. During the DFT calculations, we assume that the ground state is nonmagnetic after excluding the materials with magnetic elements and large total magnetic moment. The final database which may contains antiferromagnets is also double-checked to make sure they are indeed nonmagnetic.

DATA AVAILABILITY
The data generated during this study are available from the corresponding author on reasonable request.

CODE AVAILABILITY
The codes used in this study are available from the corresponding author on reasonable request.  f-h Light excitation from band n − 1 to bands n + 3 and n + 4 in k z = 0.31 2π/c plane with ω = 1.8eV , from band n − 1 to bands n + 7, n + 8, n + 9 and n + 10 in k z = 0.13 2π/c plane with ω = 1.95eV , and from band n to bands n + 7, n + 8, n + 9 and n + 10 in k z = 0.13 2π/c plane with ω = 1.95eV , respectively. n is the number of electrons in the primitive cell. The blue and red arrows in f-h represent positive and negative contributions to σ zzz .