Introduction

Searching for elementary particles is originally a fundamental subject in high-energy physics, but recently attracts great interest in condensed matter physics due to the discovery of Majorana1,2,3, Dirac4,5,6, and Weyl fermions7,8,9,10,11 as the low-energy excitations in solid-state systems. Distinguished with Dirac fermions, Weyl fermions, derived from the solution of massless Dirac equation (so-called Weyl equation), are massless but with a definite chirality. Until now, as an intended elementary particle, the long-sought-after Weyl fermions have not yet been observed in particle physics. However, for the last decade, it is shown that the linear low-energy excitations near the twofold-degenerate points (Weyl nodes) in electronic band structures of solids can be described by the Weyl equation. As a consequence, the quasiparticles near the Weyl nodes are referred to as Weyl fermions. Besides the high mobility arising from their linear dispersion, Weyl fermions present some unusual properties, such as nontrivial surface Fermi arc states, chiral anomaly induced negative magnetoresistance effect, and anomalous Hall effect with broken time-reversal symmetry \({\mathcal{T}}\), thus arousing extensive effort.

Topologically protected Weyl nodes in solids behave as the sink or source of Berry flux and can be regarded as the magnetic monopoles of Berry curvature field. According to the “no-go theorem”12,13, Weyl nodes come in pairs with opposite chiralities and usually carry the topological chiral charge of ±1 given by Chern number \({\mathcal{C}}\). We hereafter denote this conventional Weyl node as “single-Weyl” node for convenience. However, unlike the elementary particle obeying the Poincare symmetry in nature, the quasiparticles in condensed matter systems are constrained by the space group symmetries of solid. In view of the diversity of space group symmetries, it is therefore likely to realize unconventional quasiparticles without high-energy counterparts in realistic materials. For example, two linearly dispersive Weyl nodes with the same chiral charge can become bound together at one point as constrained by the SU(2) symmetry14,15,16,17, resulting in linear fourfold-degenerate band-crossings with a double chiral charge of single-Weyl excitations, i.e., ±2. These band-crossings, possessing the ±2 chiral charge, are usually referred to as double-Weyl nodes. Another type of double-Weyl is the double-degenerate band-crossings with quadratic dispersions18,19, associating with the topological chiral charge of ±2. This unique quadratic double-Weyl can emerge at the high-symmetry points and lines in momentum space under the constraint of certain crystal symmetry. In other words, this double-Weyl is protected by the corresponding crystal symmetries and can split into two single-Weyl excitations by breaking these crystal symmetries. For instance, the inclusion of spin-orbit coupling (SOC) can lead to the splitting of the double-Weyl node18,19.

A typical hallmark of Weyl excitations is the topologically protected surface arc states7,8,9,10,11,20,21,22,23,24,25, which are terminated by the projections of Weyl nodes onto the surface Brillouin zone (BZ). When the quadratic double-Weyl and linear single-Weyl coexist in solids, one common case is that the surface states are terminated by the projections of Weyl nodes with opposite chiralities and equal values of topological charge, as shown in Fig. 1a. Another possible case is that the surface arc states are connected with opposite chiralities but unequal values of topological charge, such as the arc between Weyl node carrying +2 topological charge and Weyl node with −1 topological charge (see Fig. 1b). Therefore, in this case, one double-Weyl node and two single-Weyl nodes are grouped as a three-terminal Weyl complex with two surface arcs, guaranteeing that the no-go theorem still holds. Wang et al26. realize this exotic three-terminal Weyl complex very recently, and the corresponding double-Weyl nodes are protected by the combination of C3 and time-reversal symmetries or the C6 symmetry in the trigonal and hexagonal crystal systems. However, the physics of the three-terminal Weyl complex remains incomplete and many issues still have to be resolved. For example, is it possible for three-terminal Weyl complex to appear in other crystal systems? How does the three-terminal Weyl complex evolve when the key symmetries are broken? Motivated by these essential questions, we start our work.

Fig. 1: Schematic of the surface arc states.
figure 1

Schematic of the surface arc states terminated by the projection of Weyl nodes with equal (a) and unequal (b) topological charge. The large and small spheres denote the Weyl nodes with chiral charge ±2 and ±1, respectively.

To realize these unique surface arc states, the key is to find out quadratic double-Weyl in solids. As discussed above, the quadratic double-Weyl nodes in electronic systems can readily split into two single-Weyl nodes due to the SOC effect18,19. Unlike the case in electronic systems, the phonon systems with negligible SOC effect, not constrained by the Pauli exclusion principle, provide a feasible platform to realize the unique surface arc states in a wide frequency range. In this work, we propose that the quadratic double-Weyl nodes with the higher chiral charge +2 emerge in the phonon of SrSi2. This quadratic double-Weyl is demonstrated to be protected by the screw rotational symmetry \({\widetilde{C}}_{4}\). As a result, the three-terminal Weyl complex with two surface arcs connecting the Weyl nodes with unequal topological charge is realized.

Results and discussion

SrSi2 crystallizes in a cubic Bravais lattice and belongs to the nonsymmorphic space group of P4332 (no. 212). Each unit cell contains four Sr atoms and eight Si atoms, in which the Si atoms form three-dimensional three-connected nets, as shown in Fig. 2a. The optimized lattice constant is a = 6.564 Å, which is in good agreement with the experimental value (a = 6.535 Å)27. The calculated Wyckoff Positions of Sr and Si atoms are (1/8, 1/8, 1/8) and (0.4223, 0.4223, 0.4223), showing excellent consistency with the experimental results27. The space group of SrSi2 has twofold screw rotational symmetries \({\tilde{C}}_{2}\) along the [011], [01\(\bar{1}\)], [101], [10\(\bar{1}\)], [110], and [1\(\bar{1}\)0] directions and fourfold screw rotational symmetries \({\tilde{C}}_{4}\) along the [100], [010], and [001] directions, which are respectively responsible for the linear single-Weyl and quadratic double-Weyl phonons, consequently leading to the unique surface arc states. The bulk BZ and the projected (001) surface BZ of SrSi2, as well as the corresponding high-symmetry points, are schematized in Fig. 2b.

Fig. 2: Crystal structure, BZs, and phonon spectrum of SrSi2.
figure 2

a Crystal structure of SrSi2. b Bulk BZ and the (001) surface BZ of SrSi2, where the high-symmetry points are marked with black and teal dots. The quadratic double-Weyl and linear single-Weyl are marked with teal (large) and orange (small) spheres. c Phonon spectrum of SrSi2, in which the two crossing branches forming the double-Weyl W1 and single-Weyl W2 are denoted by the corresponding irreducible representations of C4 and C2, respectively.

As shown in Fig. 2c, we calculated the phonon spectrum of SrSi2 in terms of the finite displacement method. In the low-frequency range of <8 THz, there are many trivial or nontrivial phonon band-crossings of the acoustic and optical phonons irrelevant to our topic. Here, we mainly pay attention to the branches of optical phonons in the high-frequency range. We can observe two twofold-degenerate band-crossings (Weyl nodes) at the frequency of 9.432 THz. One is (0.3674, 0.0, 0.0) Å−1 along Γ-X, and another is (0.2431, 0.2431, 0.0) Å−1 along Γ-M, denoted as W1 and W2, respectively. Along the Γ-X direction, i.e., the qx axis, all momentum points are invariant under the symmetry operator of fourfold rotational symmetry C4[100] along [100], and the little group therefore is C4. As a result, the two crossing branches belong to the irreducible representations Γ3 and Γ4 of C4[100]. Along the Γ-M direction, twofold rotational symmetry C2[110] along [110] is present. Thus we can denote these two crossing branches as the irreducible representations Γ1 and Γ2 of C2[110] (see Fig. 2c). Considering the C2 and C4 rotational symmetries in SrSi2, there should be 6 Weyl nodes corresponding to W1 along the C4 axes (qx, qy and qz) and 12 Weyl nodes corresponding to W2 along the C2 axes, as depicted in Fig. 2b.

To further explore the topological nature of these Weyl nodes, we plot the three-dimensional (3D) representations of phonon spectra around the W1 and W2 nodes, as shown in Fig. 3a, b. From Fig. 2c, it seems that W1 and W2 Weyl nodes cross linearly. However, the 3D phonon spectrum in the (100) plane across W1 shows that the W1 band-crossing is quadratic along qy and qz axes despite the linear dispersions along the qx axis, suggesting that W1 is a quadratic double-Weyl phonon node. Meanwhile, Fig. 3b indicates that W2 is indeed a conventional single-Weyl phonon node with linear dispersions. To determine the chiralities of these Weyl nodes, we employ the Wilson loop approach28,29 to calculate the Chern number \({\mathcal{C}}\) on a sphere enclosing the Weyl node by tracking the evolution of the hybrid Wannier charge centers (HWCCs). We find that the double-Weyl W1 contributes a Chern number \({\mathcal{C}}=+2\) and the Chern number of single-Weyl W2 is \({\mathcal{C}}=-1\) as expected. In other words, the double-Weyl W1 behaves as a source of the Berry curvature field, whereas the single-Weyl W2 acts as a sink of the Berry curvature field.

Fig. 3: 3D phonon spectra and the evolution of HWCCs.
figure 3

3D representation of phonon spectra around the W1 (a) and W2 (b) nodes onto the (100) plane. The sum of HWCCs as a function of azimuthal angle θ for W1 (c) and W2 (d).

Next, we show a symmetry analysis of these single-Weyl and double-Weyl nodes based on the two-band k p model. The two crossing branches around the Weyl node K can be described by the effective Hamiltonian

$${\mathcal{H}}({\bf{K}}+{\bf{q}})=m({\bf{q}}){\sigma }_{+}+{m}^{* }({\bf{q}}){\sigma }_{-}+n({\bf{q}}){\sigma }_{3},$$
(1)

where σ± = σ1 ± iσ2, and σj(j = 1, 2, 3) are the Pauli matrices. m(q) is a complex function and n(q) denotes a real function, in which q = (qx, qy, qz) are three components of the momentum q. The kinetic term of the k p Hamiltonian has been ignored because it is irrelevant with the band-crossing.

Along the Γ-X direction, SrSi2 possesses the screw rotational symmetry \({\tilde{C}}_{4[100]}=\{{C}_{4[100]}| {\bf{t}}\}\), where C4[100] is the fourfold rotational symmetry and t = (1/4, 3/4, 3/4) is the fractional translational vector in units of the lattice constant a. The eigenvalues of \({\tilde{C}}_{4[100]}\) can be given by \({\tilde{g}}_{{\alpha }_{k}}={g}_{{\alpha }_{k}}{e}^{-i{q}_{x}a/4}\), where \({g}_{{\alpha }_{k}}={e}^{i{\alpha }_{k}\pi /2}({\alpha }_{k}=0,1,2,3)\) are the eigenvalues of C4[100]. Under the \({\tilde{C}}_{4[100]}\) symmetry, the k p Hamiltonian satisfies

$${\tilde{C}}_{4[100]}{\mathcal{H}}({\bf{q}}){\tilde{C}}_{4[100]}^{-1}={\mathcal{H}}({R}_{4[100]}{\bf{q}}),$$
(2)

and R4[100] represents the rotational matrix of C4[100]. When two phonon branches with eigenvalues of \({g}_{{\alpha }_{1}}\) and \({g}_{{\alpha }_{2}}\) cross on this \({\tilde{C}}_{4[100]}\)-invariant path, the matrix representation of \({\tilde{C}}_{4[100]}\) can be expressed as

$${\tilde{C}}_{4[100]}=\exp \left(i\pi \frac{{\alpha }_{1}+{\alpha }_{2}}{4}\right)\exp \left(i\pi \frac{{\alpha }_{1}-{\alpha }_{2}}{4}{\sigma }_{3}\right){e}^{-i{q}_{x}a/4}.$$
(3)

Note that here α1 ≠ α2, otherwise these two branches will open a gap. Substituting \({\tilde{C}}_{4[100]}\) with Eq. (3), we can simplify the left of Eq. (2) into

$$\begin{array}{lll}{\tilde{C}}_{4[100]}{\mathcal{H}}({\bf{q}}){\tilde{C}}_{4[100]}^{-1}&=&\displaystyle m({\bf{q}})\exp \left(i\pi \frac{{\alpha }_{1}-{\alpha }_{2}}{2}\right){\sigma }_{+}\\ &&+\;\displaystyle {m}^{* }({\bf{q}})\exp \left(-i\pi \frac{{\alpha }_{1}-{\alpha }_{2}}{2}\right){\sigma }_{-}\\ &&+\;n({\bf{q}}){\sigma }_{3}.\end{array}$$
(4)

For convenience, we can choose the basis of (q+, q), where q± = qy ± iqz. Consequently, R4[100]q gives

$${R}_{4[100]}({q}_{+},{q}_{-})=({e}^{i\pi /2}{q}_{+},{e}^{-i\pi /2}{q}_{-}).$$
(5)

Based on Eqs. (4) and (5), we can obtain the symmetry-constrained expressions of m(q) and n(q) as follows:

$$\begin{array}{lll}{e}^{i\pi ({\alpha }_{1}-{\alpha }_{2})/2}m({q}_{+},{q}_{-})=m({q}_{+}{e}^{i\pi /2},{q}_{-}{e}^{-i\pi /2})\\ n({q}_{+},{q}_{-})=n({q}_{+}{e}^{i\pi /2},{q}_{-}{e}^{-i\pi /2}).\end{array}$$
(6)

Clearly, the detailed expressions of m(q) and n(q) depend on \({e}^{i\pi ({\alpha }_{1}-{\alpha }_{2})/2}\). It is easy to find that the possible values of \({e}^{i\pi ({\alpha }_{1}-{\alpha }_{2})/2}\) are ±i and −1. To determine the detailed expressions of m(q) and n(q), we can expand them with the basis of (q+, q) as follows:

$$\begin{array}{lll}m({q}_{+},{q}_{-})=\mathop{\sum }\limits_{{l}_{1}{l}_{2}}{M}_{{l}_{1}{l}_{2}}{q}_{+}^{{l}_{1}}{q}_{-}^{{l}_{2}}\\ n({q}_{+},{q}_{-})=\mathop{\sum }\limits_{{l}_{1}{l}_{2}}{N}_{{l}_{1}{l}_{2}}{q}_{+}^{{l}_{1}}{q}_{-}^{{l}_{2}}.\end{array}$$
(7)

In the case of \({e}^{i\pi ({\alpha }_{1}-{\alpha }_{2})/2}=\pm i\), the lowest order terms of m(q) constrained with \({\tilde{C}}_{4[100]}\) symmetry from Eqs. (6) and (7) are M10 or M01. We therefore write out

$$m({\bf{q}})={a}_{\pm }{q}_{\pm },\qquad n({\bf{q}})={a}_{x}{q}_{x},$$
(8)

which obviously supports the linear dispersions of the two crossing phonon branches. In the case of \({e}^{i\pi ({\alpha }_{1}-{\alpha }_{2})/2}=-1\), from Eqs. (6) and (7), the lowest order terms of m(q) that can survive are M20 and M02. Thus, m(q) and n(q) can be described as

$$m({\bf{q}})={b}_{+}{q}_{+}^{2}+{b}_{-}{q}_{-}^{2},\qquad n({\bf{q}})={a}_{x}{q}_{x},$$
(9)

which definitely implies a quadratic double-Weyl node that possesses the linear dispersions along the qx direction and the quadratic dispersions along the qy − qz plane, suggesting that these quadratic double-Weyl nodes are protected by the \({\widetilde{C}}_{4}\) symmetry.

As discussed above, the two crossing bands associated with the quadratic double-Weyl W1 are identified as the irreducible representations Γ3 and Γ4 of C4 with eigenvalues of ±i, which corresponds to \({e}^{i\pi ({\alpha }_{1}-{\alpha }_{2})/2}=-1\). Therefore, the Weyl nodes along the qx direction are definitely quadratic double-Weyl nodes, which should be protected by the screw rotational symmetry \({\widetilde{C}}_{4}\).

Along the Γ-M direction, the screw rotational symmetry \({\tilde{C}}_{2[110]}=\{{C}_{2[110]}| {\bf{t}}\}\) (C2[110] is the twofold rotational symmetry along the [110] direction) is present in SrSi2. The eigenvalues of \({\tilde{C}}_{2[110]}\) are expressed as \({\tilde{g}}_{{\beta }_{k}}={g}_{{\beta }_{k}}{e}^{-i({q}_{x}+{q}_{y})a/2}\), in which \({g}_{{\beta }_{k}}=\pm\! 1\) are the eigenvalues of C2[110]. Under the \({\tilde{C}}_{2[110]}\) symmetry, the two-band k p are constrained by

$${\tilde{C}}_{2[110]}{\mathcal{H}}({\bf{q}}){\tilde{C}}_{2[110]}^{-1}={\mathcal{H}}({R}_{2[110]}{\bf{q}}),$$
(10)

which can give

$$\begin{array}{lll}-m({q}_{+},{q}_{-})=m(-{q}_{+},-{q}_{-})\\ n({q}_{+},{q}_{-})=n(-{q}_{+},-{q}_{-})\end{array}$$
(11)

when [110] is assigned as qx. Therefore, only considering the lowest terms of m(q) and n(q), we can obtain

$$m({\bf{q}})={c}_{+}{q}_{+}+{c}_{-}{q}_{-},\qquad n({\bf{q}})={a}_{x}{q}_{x}.$$
(12)

This expression shows that the Weyl node along Γ-M is a single-Weyl with linear dispersions, and these single-Weyl nodes are protected by the \({\widetilde{C}}_{2}\) symmetry.

To further examine the \({\widetilde{C}}_{4}\)-symmetry protected double-Weyl in SrSi2, we applied a 0.5% uniaxial tensile strain along z-axis, as shown in Fig. 4a. The results (Fig. 4b) show that the double-Weyl nodes along the qz direction are robust but with slightly moved positions. The reason is that the C4 rotational symmetry along the qz direction is still present under this uniaxial strain, which guarantees these two double-Weyl nodes. However, it is found that each double-Weyl node along qx and qy splits into two single-Weyl nodes, symmetrically located with respect to the qx = 0 or qy = 0 plane (see Fig. 4b). For example, the double-Weyl W1 splits into two single-Weyl (0.3739, 0.0731, 0.0) Å−1 and (0.3739, −0.0731, 0.0) Å−1 at 8.685 THz carrying the topological charge \({\mathcal{C}}=+1\), denoted as W1\({}^{\prime}\) and W1, respectively. This can be attributed to the fact that the C4 rotational symmetries along qx and qy are broken under this uniaxial strain. Furthermore, under the uniaxial strain along the z-axis, all the original single-Weyl points still exist. However, because the C2 rotational symmetries along [011], [01\(\bar{1}\)], [101] and [\(\bar{1}\)01] are broken, the positions of the single-Weyl nodes along these directions move away from these high-symmetry lines slightly. Meanwhile, the single-Weyl nodes along [110] and [1\(\bar{1}\)0] are robust because the C2 rotational symmetries are preserved along these directions.

Fig. 4: The Weyl points of SrSi2 on the tensile strain and the surface states.
figure 4

a Schematic of the uniaxial tensile strain along z-axis of SrSi2. b The distribution of the nodes of SrSi2 under this uniaxial tensile strain, where the double-Weyl and single-Weyl are marked with teal (large) and orange (small) spheres, respectively. c The (001) surface states along the high-symmetry paths and surface arcs at different frequencies.

The surface states and constant-energy arc states provide a straightforward way to observe the topological properties of materials in the experiment. We here employ the iterative Green function method30 to calculate the (001) surface states along high-symmetry paths and iso-frequency arc states, as shown in Fig. 4c. The (001) surface states are visible, which facilitates the experimental observation by the high-resolution electron energy loss spectroscopy31 or other powerful techniques20. Moreover, the (001) surface states further show the linear dispersions of the single-Weyl node W2 and the quadratic dispersions of the double-Weyl node W1, showing consistency with the 3D phonon spectra. The iso-frequency arc states are shown at five different frequencies. Since SrSi2 possesses a cubic lattice, which makes an overlap of the arc states from different Weyl nodes. However, we can still identify the topologically protected surface arcs terminated with the projections of the double-Weyl nodes and two single-Weyl nodes. It is worth noting that the surface arcs shown here are distinctly different from the surface state of conventional Weyl quasiparticle, that is, the surface arcs here connect the Weyl quasiparticles with unequal charge, consequently leading to the three-terminal Weyl complex with two surface arcs. In addition, due to the presence of C4 symmetry, we can see that the arc states onto the (001) plane present a distinct tetragonal symmetry with the increase of frequency.

In summary, using first-principles calculations and symmetry analysis, we show that the unconventional quadratic double-Weyl phonons can be realized in SrSi2. Based on the two-band k p Hamiltonian, we demonstrate that these quadratic double-Weyl are protected by the fourfold screw rotational symmetry \({\widetilde{C}}_{4}\). The perturbation of 0.5% uniaxial tensile strain further confirms the \({\widetilde{C}}_{4}\)-protected double-Weyl. Moreover, the Weyl complex, consisting of one quadratic double-Weyl and two linear single-Weyl, presents unique surface arc states, which are different from conventional paired Weyl nodes. Distinguished with the three-terminal Weyl complex in the trigonal and hexagonal crystal systems in which the time-reversal symmetry may be needed to protect the quadratic double-Weyl, our results suggest that it is highly probable to realize this exotic Weyl complex protected by the fourfold screw rotational symmetry in ferromagnetic fermionic systems in the absence of time-reversal symmetry. In addition, SrSi2 is also a Weyl semimetal candidate with topological electronic surface states19, and therefore the topological phonon surface states predicted here can potentially couple with the electronic surface states and enhance the electron-phonon coupling on the surfaces of SrSi2, consequently leading to exotic physical phenomena.

Methods

All the first-principles calculations were carried out within the framework of density functional theory32,33 by using VASP34. The Perdew–Burke–Ernzerhof functional35 within the generalized gradient approximation was employed to describe the exchange-correlation interactions. We adopted the projector augmented-wave method36 to describe the electron and core interactions. This method treats the 4s24p65s2 and 3s23p2 electrons as the valence electrons for Sr and Si atoms in SrSi2, respectively. The cutoff energy for the plane wave was 500 eV, and the BZ for structural optimization was sampled with a grid of 8 × 8 × 8 by employing the Monkhorst–Pack scheme37. The phonon spectrum of SrSi2 was calculated by using the finite displacement method with the 2 × 2 × 2 supercell as implemented in Phonopy code38. To calculate the topological properties of the phonons for SrSi2, we constructed a phonon tight-binding Hamiltonian based on the phonon force constants. On the basis of the obtained tight-binding Hamiltonian and the iterative Green function method, we calculated the surface states and arcs using the WannierTools package 39.