High-throughput screening strategy of topological phonons

As shown in several prototypical materials^23–42, identifying TPs requires several stages of manual selections and subjective human decisions. To enable the TPs discovery in an automatic manner, we present a HT screening and data-driven approach to discover and categorize TPs, as described in Fig. 1, including the following four steps.

Fig. 1

The schematic flowchart of high-throughput computational screening on topological phonons.

This workflow is capable of identifying the features of TPs, elucidating details of topology and constructing TPs database by computing and collecting phonons of a variety of solid materials in an automatic manner.

Phonon data collection. To obtain phonon spectra for a large volume of known materials, we first collected 10,000 materials’ force constant data from public phonon database^58,59. The data set was further augmented by our in-house computations for over 2000 materials belonging to 58 common structural prototypes. It is well known that the calculated force constants are numerically sensitive to the choices of several parameters, e.g., the supercell size, K-point mesh, and energy cutoff. To guarantee that the predictions are reliable, we filtered out the materials with notable imaginary frequencies (<−0.5 THz) in the whole phononic momentum space.

Nodal straight lines identification. We computed their band dispersions on the automatically generated high-symmetry band paths⁶⁰. If there exist degenerate phononic bands along high-symmetry paths,we would compute the Berry phase $$ {\gamma }_n={\oint }_c\mathrm{d}l\cdot {\mathcal{A}}_n$$, by an integral of Berry connection ($$ {\mathcal{A}}_n\left(q\right)=i<{\mu }_{n,q}\mid {\nabla }_q\mid {\mu }_{n,q}>$$) over a closed q path⁶¹, for 20 consecutive points on each of these bands. The bands possessing continuous points with Berry phase values, ±π, would be marked as the topologically nontrivial nodal straight lines.

Crossing points screening. For the rest band paths in the phonon spectra, we systematically scanned 50 points on each band path. In the entire frequency range, we considered the points possessing two adjacent eigenfrequencies < 0.5 THz. For each point, we performed a minimization based on the conjugate gradient algorithm to obtain the local minimum of the frequency difference (Δ_freq). The points with Δ_freq < 0.2 THz were checked if they possess with Berry phase values of ±π. After optimization, the identified crossing points may go anywhere in the entire reciprocal space. Therefore, we also checked if the points are at or off the high-symmetry paths.

Crossing points assignment. The identified phononic crossing points were then divided into two groups based on the presence of both inversion symmetry (P) and TRS (T) for each material. In a three-dimensional (3D) system with PT symmetries, the Berry curvatures of nondegenerate phononic bands are forced to be zero and the Weyl TPs would not occur in such system. Once the phononic bands at a degenerate point have opposite nonzero Berry phases, such topological nontrivial degenerate points have to occur continuously by forming nodal-line (ring) TPs, due to the continuity of phonon wave function in the 3D momentum space. As a result, when the PT symmetry is present, we just need to seek nodal-lines (rings) off high-symmetry line. In noncentrosymmetic materials, the phonon dispersions possibly form single Weyl or high degenerate Weyl TPs, in addition to nodal-lines (rings) TPs. In order to clarify these three types of TPs, we introduced another formula of Chern number, which can be derived by integrating Berry curvatures of a closed surface⁶¹ according to $$ n=\frac{1}{2\pi }{\int }_{\mathit{S}}\mathrm{d}\mathit{S}\cdot \mathrm{\Omega }\left(\mathit{q}\right).$$ Here, S is a closed surface which wraps the target crossing point, and the Ω(q) is the Berry curvature at the phonon momentum q on the selected closed surface. For the isolated crossing points at the high-symmetry band paths, we marked the points with nonzero integer Chern numbers (e.g., ±1, ±2) as single or high degenerate Weyl points. Otherwise, they would be labeled as nodal-ring points, given the fact that the nodal-line points were already extracted in step (2). Of course, it needs to be emphasized that many materials may yield multiple crossing bands along off high-symmetry paths, which can be separate Weyl TPs or nodal-line (ring) TPs (Supplementary Table 1). In principles, this approach can be easily extended to investigate these crossing points.

In order to experimentally observe phononic surface states, a TP material is expected to possess distinct Weyl or nodal-line (ring) TPs and possible clean nontrivial surface TPs. For this purpose, we mathematically define the clean TPs for the nontrivial crossing points satisfying two conditions in bulk phonon spectra, (i) the crossing points have to be located at local minima with negligible or zero phononic density of states (DOS; <0.01 states/atom/THz) and (ii) the dispersion at the local minima is sufficiently large (∂E/∂q > 3.0 THz ⋅ Å). On basis of these two criteria, we have filtered 322 clean TP materials (Supplementary Table 5). For instance, single type-I Weyl TPs in LiCaAs are clean, because it does not overlap with the other bulk phonon branches and has a zero phonon density at the frequency of 4.590 THz (Fig. 2i, j).

Fig. 2

Phonon dispersion and single Weyl TPs of LiCaAs.

a The unit cell and primitive cell of LiCaAs (space group $$ F\overline{4}3m$$ 216). b The phonon dispersion along the high-symmetry lines. The blue circles denote the phononic crossing points (the single Weyl TPs). c The bulk BZ and the (111) surface BZ. d, e The 3D phonon dispersions centered at the Weyl point on the q_xy and q_yz planes, respectively. f, g The Wannier center evolution and Berry curvatures distributions around this Weyl phonon. h The surface phononic states along the high-symmetry lines. The red and yellow circles represent the projected Weyl TPs with positive and negative topological charges of 1 and −1, respectively. i, j the Weyl phonon induced nontrivial surface phononic arc states at 4.50 and 4.43 THz, respectively.

Topological phonons of materials

In total, our approach revealed that 5014 materials exhibit TPs states (Supplementary Table 1). Among them, we have identified two main categories of nodal-line (ring) and Weyl TP materials. Among nodal-line (ring) TPs (Supplementary Table 2), there are several possible subclasses, including nodal-link, nodal-chain, and nodal-net TPs according to different symmetries. Among Weyl TPs, there are two main subclasses of single Weyl TPs (Supplementary Table 3) and high degenerate Weyl TPs (Supplementary Table 4) according to the degree of phononic band degeneracy. In particular, we note that, different from electronic system, it is impossible to have the intrinsic TP insulators without any tunable external field. This is mainly because the phonon spectrum always preserves the time-reversal symmetry (TRS). The TRS for phonon spectrum was suggested to be broken only in artificial lattices, such as in ionic lattices using Lorentz force³⁵ and magnetic lattices using spin–lattice interactions⁶². Therefore, we did not attempt to identify the intrinsic topological phonon insulators in our current work.

To understand the occurrence of TPs in various materials and their correlations, we chose to investigate four representative cases, including (1) half-Heusler LiCaAs alloy for single Weyl TPs, (2) superconducting BeAu for high degenerate Weyl TPs, (3) ScZn with the PT symmetries for nodal-line (ring) TPs, and (4) TeO₃ with both PT and nonsymmorphic symmetries for HNN TPs. To conveniently analyze their TP states, we use a general k ⋅ p Hamiltonian as follows,

Applying different constraints of crystal and PT symmetries to this Hamiltonian, we can describe various geometries of the phonon crossings.

Single Weyl TPs in noncentrosymmetric half-Heusler alloys

Half-Heusler compounds A^IB^IIX^V (where A ∈ {Li, Na, K, Rb, Cs}; B ∈ {Mg, Ca, Zn}, and X ∈ {P, As}) with noncentrosymmetric space group F43m (no. 216) have been widely studied ^63,64. Since they share similar phonon dispersions, here we focus on the case of LiCaAs. As shown in Fig. 2a, the As atoms are in the 4a Wyckoff sites, while both Li and Ca are in the 8c sites. At 4.590 THz, the highest longitudinal acoustic (LA) and the lowest transverse optical (TO) branches have a band cross at the $$ \left(0.435,0.0,1.0\right)\frac{2\pi }{a}$$ point along the W–X high-symmetry path Fig. 2b. Interestingly, the phononic spectra around this crossing point exhibit a Weyl cone-like shape in both q_xy and q_yz planes, as shown in Fig. 2d, e. To further confirm its topological nature, the Wannier center evolution has been derived for the third phononic band (Fig. 2f). It gives a topological charge of −1 and this crossing point acts like the sink of Berry curvatures in Fig. 2g. These results indicate that the crossing point along the X–W path is an ideal type-I Weyl TP with the topological charge of ±1.

To elucidate the underlaying physics for the type-I Weyl TPs in LiCaAs, we constructed the Hamiltonian from Eq. (1) according to the fact that LiCaAs possesses both the TRS T and the twofold rotational symmetry C₂. As shown in Fig. 2d, e, the type-I WPs do not have the tilt term (expressed as d₀σ₀ in Eq. (1)). Combining T and $$ C_2^z$$ symmetries, the operator can be represented by σ_zκ (where κ is the complex conjugate operator). The Hamiltonian

This condition implies that the crossing points on the square plane ($$ q_z=\frac{2\pi }{a}$$) are constrained to the high-symmetry lines along either q_x or q_y direction. As a result, the type-I Weyl TPs on the W–X line of the BZ boundary are protected by both C₂ and T symmetries, which can produce well-separated WPs and give rise to the extremely long open arcs on the surface. The topologically protected surface states of (111) surface along high-symmetry line have been derived in Fig. 2h and the Weyl points marked as the red and yellow spheres are projected on the $$ \overline{L}$$–$$ \overline{X}$$ line. The open surface arcs have also been clearly observed at 4.43 and 4.50 THz (Fig. 2i, j), respectively. Each pair of projected Weyl points with opposite charges is clearly connected by the open surface arc, which is long enough and can provide a robust one-way phonon propagation channel on the surface without backscattering from defects. To our best knowledge, the Weyl TPs in LiCaAs exhibit two unique features different from the known Weyl TPs. First, this type-I Weyl point is the only phononic band crossing between the LA and TO branches. Second, it is the only known clean Weyl TP, which doesn’t overlap with any other phonon branches. The calculations further reveal that the family of half-Heusler A^IB^IIX^V compounds host similar Weyl TPs and their nontrivial long arc surface states. These materials are suitable for detecting both single Weyl TPs and nontrivial long arc surface phonons by experiments.

High degenerate Weyl TPs in noncentrosymmetric superconductor BeAu

High degenerate Weyl points refer to the nontrivial crossings, which have a degeneracy higher than two. Our screening has found the existence of threefold and fourfold nontrivial high degenerate Weyl points in 447 TP materials (see Supplementary Table 4). For threefold degenerate Weyl TPs, the Hamiltonian can be written as H₃(k) ∝ q ⋅ S, where q is a wavevector and S_i are the rotation generators for spin-1 bosons. Those three bands at the crossing point will have the Chern numbers of +2, 0, and −2 and they are also called “spin-1 Weyl point”²⁴. For fourfold degenerate Weyl TPs, their Hamiltonian can be written as H₄(k)~I₂⨂(q ⋅ σ), where I₂ is the second-order identity matrix and σ_i(i = x, y, z) are the three Pauli matrices. Those high degenerate Weyl TPs can be regarded as the sum of identical spin-$$ \frac{1}{2}$$ Weyl points and they can be called “charge-2 Dirac point”²⁴ due to their Chern numbers.

Among 447 TP materials with high degenerate Weyl TPs, BeAu is a typical noncentrosymmetric superconductor with a critical temperature of 3.2 K (ref. ⁶⁵). BeAu crystallizes in a cubic B20 structure with the space group of P2₁3 (no. 198). Both Be and Au atoms are located at 4a Wyckoff sites (Fig. 3a). The phonon spectrum along the high-symmetry lines in Fig. 3b shows six threefold degenerate Weyl TPs (blue circles at Γ) and six fourfold degenerate Weyl TPs (green circels at R). These high degenerate Weyl TPs are protected by both the lattice symmetries (twofold screw rotations and threefold rotations) and TRS. To determine their topological natures, we calculate the Chern numbers for high degenerate Weyl TPs at both Γ and R in the BZ (Fig. 3c, d). The Chern numbers for threefold degenerate Weyl TPs are −2, 0, and 2, while those are ±2 for fourfold degenerate Weyl TPs. For the Weyl TPs at Γ, we can derive its Wannier center evolutions because of the well-separated bands. The threefold degenerate Weyl TPs in Fig. 3c are contributed from three phonon branche nos. 10, 11, and 12. As shown in Fig. 3f, the Wannier center evolutions for those three branches indicate that both nos. 10 and 12 bands are topologically nontrivial, whereas no. 11 band is trivial. This fact reveals that the threefold degenerate point is similar to a single Weyl point. However, it has a high topological charge of +2 and the corresponding Berry curvatures at q_xy plane give the source behaviors at Γ (ω = 3.669 THz), which is protected by the twofold screw rotation axis at Γ in the BZ. The fourfold degenerate Weyl TPs in Fig. 3d are contributed from two doubly degenerate bands, with the Chern numbers, c = ±2. These fourfold degenerate Weyl TPs at R have the topological charge −2, which is also protected by the twofold screw rotation axis at R. Both the threefold and fourfold degenerate Weyl TPs have the opposite topological charges, indicating the conservation of the topological charges for Weyl TPs in the first BZ. Meanwhile, as shown in Fig. 3h, we have derived the topological nontrivial surface TPs along the high-symmetry lines on the (001) surface. Clearly, the topological nontrivial surface states connect both threefold and fourfold degenerate Weyl TPs at $$ \overline{\mathrm{\Gamma }}$$ and $$ \overline{M}$$. Consequently, it can be clearly seen that two long arcs connect the projected double Weyl points at center-$$ \overline{\mathrm{\Gamma }}$$ and corner-$$ \overline{M}$$, and those surface arcs are constrained by the TRS. The nontrivial open arc phononic surface states can be used to design the ideal negative refraction materials and this concept was recently realized in macroscopic metamaterials of Weyl phononic crystal⁶⁶.

Fig. 3

Phonon band structures and surface states for topological high degenerate Weyl TPs in BeAu.

a The crystal structure of BeAu (space group P2₁3 198). b The phonon dispersion along the high-symmetry lines. The blue circles are the threefold degenerate Weyl TPs at Γ and the green circles are the fourfold degenerate Weyl TPs at R. c Threefold degenerate Weyl point of 3.669 THz at Γ. d Fourfold degenerate Weyl TPs of 3.956 THz at R. e Bulk and surface BZ of BeAu. f The Wannier center evolution for three branches 10, 11, and 12 centered at the Γ. g The Berry curvature distributions of three branches 10, 11, and 12 at the centered Γ. h The surface local density of states for (001) surface along the high-symmetry directions. i, j The corresponding surface arcs at 3.65 and 3.55 THz, respectively. Even though BeAu exists in reality, we still found that around Γ point one acoustic branch of BeAu shows the extremely small imaginary frequency, which cannot be removed in our current calculations, possibly due to misconsideration of long range interatomic interaction in the force constant construction or anharmonic effects.

Nodal-line (ring) TPs in centrosymmetric ScZn

Nodal-line (ring) TPs, formed by continuous phononic band crossings, have also been predicted in phonon spectra of materials^30–34. Our results suggest that many materials exhibit such nodal-line (ring) TPs. Here, we introduce ScZn, with a B2 lattice structure (see Fig. 4a), that hosts both nodal-ring and straight line TPs as shown in Fig. 4c.

Fig. 4

Phonon band structures and topological properties of phonon of ScZn.

a The crystal structure of ScZn (space group $$ Pm\overline{3}m$$ 221). b The BZ of ScZn and the illustration of the nodal-ring TPs (red curve). c The phonon spectrum of ScZn. d The illustrations of the straight nodal-line TPs along the Γ–R and Γ–X lines in ScZn. e, f are the 3D phonon bands around the nodal-line TPs surrounding the M point. g, h The derived phononic surface states at the frequencies of 4.12 and 4.05 THz, respectively, of the (001) surface BZ (as defined in b). i The surface phononic spectrum of the (001) surface. j, k The phononic surface states at the frequencies of 3.60 and 3.55 THz, respectively.

Firstly, we have noted that the phononic band crossing point between LA and TO on the X–M line of the BZ (Fig. 4c). Unlike the single Weyl TP in LiCaAs, this crossing is not an isolated point, but belongs to a closed ring formed by the continuous linear band crossings, as shown in Fig. 4e, f. The occurrence of the nodal-ring TPs in ScZn can be attributed to the presence of the PT symmetry. Considering the mirror symmetry, ScZn totally hosts six nodal-ring TPs, which are located at the boundary planes of the BZ surrounding the M point in Fig. 4b. Secondly, the nodal-line TPs have been observed along the Γ–X and the Γ–R directions, as shown in Fig. 4c. They can be viewed as countless linear band crossings along the high-symmetry lines and extend through the whole BZ, as illustrated in Fig. 4d, similar to MgB₂ (ref. ³¹). In particular, our further studies reveal that the nodal-ring TPs centered at the M point and the nodal-line TPs along the high-symmetry directions in ScZn are both protected by the PT and mirror symmetries. For inversion symmetry P, we can introduce the inversion operator as $$ \widehat{P}={\sigma }_z$$ to Eq. (1) to satisfy

We can simplify this equation into the following relations

At this stage, the eigenvalues of Eq. (2) can be simplified to $$ \omega \left(\mathbf{q}\right)=d_0\pm \sqrt{d_2^2+d_3^2}$$ and the band crossing points require d₂(q) = d₃(q) = 0. Here, we take the nodal-ring TPs centered at M as an example to elaborate the role of the mirror symmetry. First, d₂(q) = α_xq_x + α_yq_y + α_zq_z = 0 can determine a plane passing the central point of the circle and the d₃(q) = b + $$ {\beta }_x{q}_x^2$$ + $$ {\beta }_y{q}_y^2$$ + $$ {\beta }_z{q}_z^2$$ = 0 is an equation of an ellipsoidal surface centered at M. The crossings between the plane and the ellipsoidal surface form a closed loop, which is the nodal ring centered at M. However, this nodal ring can tilt to arbitrary direction. Here, we can choose the mirror symmetry as $$ \widehat{M_z}={\sigma }_z$$ to give an additional constraint to the Hamiltonian,

Using σ_z to replace $$ \widehat{M_z}$$, we can obtain the following expressions constrained by mirror symmetry

Therefore, the M-centered nodal-ring TPs can be confined to the plane of q_z = π. For nodal-line TPs, the mirrors can also constrain them along one high-symmetry line. Since the nodal-line TPs possess nonzero Berry phases, there exist phononic nontrivial drumhead-like surface states. To elucidate this feature, we have calculated the phonon spectrum on the (001) surface with Green’s function method in Fig. 4i. It shows the nontrivial surface states along the $$ \overline{\mathrm{\Gamma }}$$–$$ \overline{X}$$ line and the $$ \overline{X}$$–$$ \overline{M}$$ line, respectively. When the bulk phonons are projected to the (001) surface, the nodal-ring TPs, perpendicular to the q_z direction, can hold co-dimensions, which form a closed circle around the $$ \overline{X}$$ point, as illustrated in Fig. 4b. Within the projected nodal-ring TPs, the nontrivial drumhead-like surface states form closed loops around the $$ \overline{X}$$ point, which evolve with the phonon frequency (Fig. 4g with ω = 4.12 THz and Fig. 4h with ω = 4.05 THz). When the straight nodal-lines are projected onto the (001) surface, they are still along the high-symmetry paths to form the triangle region on the surface 2D BZ in Fig. 4d. Within this region, the nontrivial drumhead-like surface phononic states would occur. The nontrivial drumhead-like surface states induced by straight nodal lines are between 3.55 and 4.00 THz, and they connect two points on the neighboring projected straight nodal lines, as shown in Fig. 4j, k.

Our HT calculations further reveal that 114 CsCl-type materials isostructural to ScZn exhibit very similar phonon dispersions (see Supplementary Table 1). All of them host the nodal-line TPs within their acoustic branches along the Γ–X and Γ–M lines, but only some of them host nodal-ring TPs between the LA and TO branches. The existence of the nodal-ring TPs depends on the relative atomic masses of the constituents in compounds. If the atomic masses differ greatly, the nodal-ring TPs at zero acoustic-optical gap will disappear.

Hourglass nodal-net TPs in centrosymmetric TeO₃

Besides the above nodal-line (ring) TPs associated with the symmorphic symmetries, the nonsymmorphic space groups can produce the symmetry-enforced nodal-line (ring) TPs. Nonsymmorphic symmetries are formed by combining the operations of the point group g and translations t by fraction. TeO₃ is such a case, which crystallizes in the Pnna space group (no. 52)⁶⁷. In its unit cell, there are four Te atoms and each of them forms the TeO₆ octahedra (Fig. 5a). In this space group, six nonsymmorphic symmetries are present as follows: g₁ = {M_z∣a/2}, g₂ = {M_y∣a/2 + b/2 + c/2}, g₃ = {M_x∣b/2 + c/2}, g₄ = {C_2y∣a/2 + b/2 + c/2}, g₅ = {C_2z∣a/2}, and g₆ = {C_2x∣b/2 + c/2}. Our results reveal that TeO₃ possesses the HNN TPs in its phonon dispersion in Fig. 5b. To show the key feature of the HNN TPs, here we focus on the four phonon branches (bands 25, 26, 27, and 28) with a frequency range from 12 to 14 THz. These four bands form two nodal-line TPs (marked in blue in Fig. 5b) along the boundaries of the BZ (Fig. 5c). The occurrence of those nodal-line TPs are protected by the crystal symmetries. By combining the TRS T and g₁ = {M_z∣a/2}, the lattice momentum can be transformed to

Fig. 5

Hourglass nodal-net TPs of TeO₃.

a The unit cell of TeO₃ (space group Pnna 52). b The phonon dispersion of TeO₃. c The bulk BZ and the (001) and (100) surface BZs. Both the red and black arrows denote the hourglass nodal-line TPs, whereas the blue lines along the BZ boundary represent the twofold degenerate nodal-line TPs. All points on the q_y = π plane are twofold degenerate nodal points (called twofold degenerate nodal surface), as shown in the yellow plane. d The phonon spectrum along X–U–R–S–X at the q_yz plane. The eigenvalues of {M_x∣b/2 + c/2} are given at both U and R points. e The 3D phonon dispersions at the q_yz plane. The solid white line represents the hourglass nodal line on which any point represents a hourglass nodal point. f The phononic hourglass nodal dispersions protected by g₁ and g₄. The eigenvalues of {M_z∣a/2} are given at both Y and X points in the BZ. g Hourglass nodal chains (HNCs) at the q_xy plane (red solid lines) and at the q_yz plane (blue solid lines). Both these HNCs connect to each other only at the S point to form an hourglass nodal net (HNN) in the extended BZ. h The shape of the HNN in the three-dimensional view. Two black circles are used to calculate the Berry phase of the HNCs and the color represents the energy dispersion of the HNN in THz. i, j The phononic surface states along the high-symmetry lines for the (001) and (100) surfaces, respectively.

It needs to be emphasized that the hourglass nodal points (Fig. 5b) are protected by crystal symmetries. Here, we focus on the hourglass nodal point along the U–R line in q_x = π plane (XURS) to elucidate the role of nonsymmorphic symmetries. For the operation g₃ = {M_x∣b/2 + c/2}, the eigenvalues are g_±(q_y, q_z) = $$ \pm e^{\left(i{q}_y/2+i{q}_z/2\right)}$$. Hence, at the point U(π, 0, π), the degenerate phonon points have opposite eigenvalues (±i). Similarly, for the point R(π, π, π), the eigenvalues would be +1 or −1. However, the g₄T operator protects the band degeneracy of the S–R line and it commutes with g₃ at the R point, which leads the same eigenvalues, (+1, +1) or (−1, −1) for the degenerate phonon bands at R. As a result, from U to R, phonon bands have to switch partners and inevitably cross each other, which causes the occurrence of the hourglass phonon bands in Fig. 5d. Interestingly, this hourglass phonon nodal point is not an isolated one and there exists an hourglass nodal ring on the XURS plane. As illustrated by the black curve starting from the S point in Fig. 5c, the hourglass nodal point along the U–R line is a point on this ring. To clearly visualize the shape of the hourglass phonon nodal ring, we plot the two crossing phononic bands on the XURS plane in Fig. 5e. Because the start and end points of this nodal ring are at the same point S in the corner of BZ, this hourglass nodal ring can form the hourglass nodal chain (HNC; see the blue curve in Fig. 5g), which goes through the whole BZ along the q_z direction. Similarly, the HNC can also be observed on the q_z = 0 plane (ΓXSY). For the glide reflection operation g₁, the eigenvalues of g₁ are $$ \pm e^{i{q}_x/2}$$. The eigenvalues at X(π, 0, 0) and Y(0, π, 0) can be identified to ±i and ±1, respectively. Since the twofold screw axis g₄ commutes with g₁, the eigenvalues at Y can be confirmed to be (+1, +1) or (−1, −1), while the X point holds the opposite eigenvalues (+i, −i), as shown in Fig. 5h. Hence, the hourglass nodal point appears along the X–Y line. On the ΓXSY plane, the similar hourglass nodal points appear to form a closed ring also starting from the S point (see the red curve in Fig. 5c). Thus, the HNC can also be observed in the extended BZ (see the red curve in Fig. 5g), which also goes through the whole BZ on the q_z = 0 plane. In particular, both HNCs originating from S are perpendicular to each other which form the HNN as shown in Fig. 5g, h.

Furthermore, the topological nontrivial surface states have been calculated on the (001) and (100) surfaces. Due to the co-dimension rule⁷¹, the surface states of the HNC on q_z = 0 plane can be obtained on the (001) surface, whereas the surface states of the HNC on the q_x = π plane emerge on the (100) surface. As the feature of nodal rings, the drumhead-like surface states can be observed inside or outside the projected region of HNCs, which are determined by the region with nonzero Berry phase³². As shown in Fig. 5i, j, the nontrivial surface states occurs outside the projected HNC on the (001) surface (Fig. 5i), while the drumhead-like surface states appear inside the HNC on the (100) surface (Fig. 5j). HNNs are protected by the nonsymmorphic symmetry and can be extremely stable. Due to the diversity of crystal symmetries, more interesting topological phonon nodal lines or rings can be expected.

Phonon calculation

All DFT calculations have been performed by Vienna ab initio simulation package, based on the projector augmented wave potentials and the generalized gradient approximation within the Perdew-Burke-Ernzerhof for the exchange correlation treatment. The force constants downloaded from the PHONONPY database were calculated by the finite displacement method. For the in-house phonon calculations, we used the density function perturbation theory. We performed the geometry optimization of the lattice constants by minimizing the forces within 0.001 eV/Å. The cutoff energy for the expansion of the wave function into the plane waves was set to 1.5 times of the ENMAX in the POTCAR. For the topological analysis, we used the conjugate gradient method in SciPy¹⁰⁷ to get the crossings, and calculated the Berry phase and Chern number to identify the nontrivial topological natures. To determine the topological charge, the Wilson-loop method^108,109 was chosen. A sphere centered at a WP was sliced into independent orbitals by a constant polar angle θ and the evolution of Wannier centers (phase factor ϕ) on orbitals can give topological charges of WPs. For the surface DOS, we used force constants as tight-binding parameters to construct surface and bulk Green’s functions, and the imaginary part of the Green’s function produces the DOS^110,111.