Hear the Sound of Weyl Fermions

Quasiparticles and collective modes are two fundamental aspects that characterize a quantum matter in addition to its ground state features. For example, the low energy physics for Fermi liquid phase in He-III was featured not only by Fermionic quasiparticles near the chemical potential but also by fruitful collective modes in the long-wave limit, including several different sound waves that can propagate through it under different circumstances. On the other hand, it is very difficult for sound waves to be carried by the electron liquid in the ordinary metals, due to the fact that long-range Coulomb interaction among electrons will generate plasmon gap for ordinary electron density fluctuation and thus prohibits the propagation of sound waves through it. In the present paper, we propose a unique type of acoustic collective modes in Weyl semimetals under the magnetic field called chiral zero sound. The chiral zero sound can be stabilized under so-called"chiral limit", where the intra-valley scattering time is much shorter than the inter-valley one, and only propagates along an external magnetic field for Weyl semimetals with multiple-pairs of Weyl points. The sound velocity of the chiral zero sound is proportional to the field strength in the weak field limit, whereas it oscillates dramatically in the strong field limit, generating an entirely new mechanism for quantum oscillations through the dynamics of neutral bosonic excitation, which may manifest itself in the thermal conductivity measurements under magnetic field.

So far the Weyl semimetal is considered as a new topological state in condensed matter physics only because of its unique quasiparticle dynamics, which manifests itself in various of transport experiments [18][19][20][21][22]. On the other hand, the unique collective mode is another type of features to characterize a new state of matter, which is yet to be revealed for Weyl semimetal systems [27][28][29][30][31][32][33][34]. The most common collective mode in a liquid system is sound, which usually requires collisions to propagate. For a neutral Fermi liquid such as He-III [35][36][37][38], the ordinary sound can only exist when ωτ 1, where τ is the lifetime of the quasiparticles. For a clean system, the low energy quasiparticle lifetime approaches infinity with reducing temperature, which prohibits the existence of the normal sound modes at enough low temperature when ωτ 1. However, there is an entirely different type of sound that emerges at the above "collision-less regime" called zero sound, which is purely generated by the quantum mechanical many-body dynamics under the clean limit [35,[39][40][41][42]. In a typical Fermi liquid system, the zero sound can be simply viewed as the deformation of the Fermi surface that oscillates and propagates in the system with the "restoration force" provided by the residual interaction among the quasiparticles around the Fermi surface. Like other types of elementary excitation in condensed matters, the form functions of zero sound carry the irreducible representations ( irreps) of the symmetry groups of the particular system. For an electron liquid in a normal metal, the density fluctuation corresponding to the trivial representation is always governed by the long-range Coulomb interaction and becomes the well-known plasmon excitation with a finite gap in the longwave limit. Thus the zero sound modes can only exist in high multipolar channels ascribing to the non-trivial representation of the symmetry group, within which the residual interactions among the quasiparticles is positive definite. The above condition requires strong and anisotropic residual interaction in solids, which is difficult to be realized in normal metals.
One of the exotic phenomena of Weyl semimetal is the chiral magnetic effect (CME) [23,[43][44][45][46], where each valley will contribute a charge current under the external magnetic field. The amplitude of such "anomalous current" is simply proportional to the occupation number of that particular valley with its direction being parallel or anti-parallel to the field direction depending on the chirality of the valley. The above CME immediately cause an interesting consequence, the particle number imbal-ance among different valleys will induce particle transport and thus makes it possible to form coherent oscillations of the valley occupation numbers over space and time, which is an entirely new type of collective modes induced by the CME.
On the other hand, the most well known collective mode in a charged Fermi liquid system is plasmon, and for a Weyl semimetal under a magnetic field, they are such collective modes where the oscillation of the valley occupation numbers cannot cancel each other exactly and thus generates net charge fluctuation in real space [27][28][29][31][32][33][34]. (However, in for neutral Weyl Fermions where the interaction is short-ranged, the oscillating CME current generates a gapless collective mode [30].) Since these modes are coupled to the CME current, the plasmon frequency significantly depends on the magnetic field [34]. Following Ref. [34], in this paper we call them "chiral plasmons" (CP). In general, CPs form the trivial irrep of the symmetry group. As discussed in detail below, among all the CPs there are only two branches are fully gapped (with opposite frequencies), whereas the other branches are gapless. For a simplest Weyl semimetal with only a single pair of WPs, the little group at finite wave vector q contains only identity operator under parallel magnetic filed indicating that all the electronic collective modes propagating with wave vector q will generally cause net charge density fluctuation and thus belong to different branches of CP modes.
The situation is completely different for a Weyl semimetal with multiple pairs of WPs. Now we can have collective "breathing modes" for Fermi surfaces in different WP valleys so that they oscillate in an anti-phase way and cancel out the net charge oscillation exactly as illustrated in Fig. 1(d) for two pairs of WPs. Since these anti-phase modes don't cause any net charge current, the collective oscillations of the valley charge and valley current will be completely decoupled from the plasmon modes, and their dispersion relation remains linear in the longwave limit, which is called "chiral zero sound" (CZS) in this paper. As introduced in more detail below, the CZS modes carry all the non-trivial irreps of the corresponding little group, with which we can figure out how many CZS modes can exist with the magnetic field being applied in some particular crystal direction.
In order to clearly describe the physical process in WSM systems, we divide the charge current contributed by the i-th Weyl valley j i into two parts, the "anomalous current" j a i caused by the change of the valley occupation number through the CME and the "normal current" j n i caused by the deformation of the ith Fermi surface. For the general situation, both the two types of current will contribute to the collective dynamics. Therefore both the plasmon and zero sound modes are not purely contributed by the anomalous current and CME. However, we can consider a specific limit, under which only the anomalous current can survive, and both the CP and CZS are purely contributed by the CME. Such a limit is proposed previously by D. Son et al. [23], requiring the intra-valley relaxation time to be much shorter than the inter-valley one, which guarantees the intra-valley relaxation process is fast enough so that any deformation of the Fermi surface from its equilibrium shape can be neglected. In the following, we will call this limit as "chiral limit" and mainly discuss the physics of CZS under it.
Let us first introduce the Boltzmann's equation in the chiral limit. The Boltzmann's method is valid only in the semi-classical limit, where ω B τ 1, ω B µ, and µτ 1.
Here ω B = v F √ eB is the magnetic frequency, e = ∓|e| is the charge of electron-like (hole-like) quasiparticle, v F is the Fermi velocity, µ is the chemical potential measured from the energy of WP, and τ is the quasiparticle lifetime. (In this paper we set = 1.) In the semiclassical limit, the level smearing is much larger than the Landau-level splitting but is much smaller than chemical potential, such that the Landau-level quantization can be ignored and the Fermi surface is well defined. The collective dynamics of these systems can be described by the quasiparticle distribution function n ν (k, r, t), and the total energy reads where the second and third terms denote the long range Coulomb and residual short range interaction among the electronic quasiparticles around the WPs respectively [47]. Here δn ν (k, r, t) = n(k, r, t) − n F ( 0 ν (k) − µ) is the deviation from Fermi-Dirac distribution, ν is the valley index, δn(r, t) = ν´d 3 k (2π) 3 γ ν (k, B)n ν (k, r, t) is the charge density at position r, γ ν (k, B) = 1 − eB · Ω ν (k) is the momentum space measure under external magnetic field B [48], where Ω ν (k) is the Berry's curvature. In general, the short-range interaction matrix f ν,ν in above equation should be written as f ν,ν (k, k ) [47]. However, here we consider the case where the Fermi surfaces are small enough such that the k-dependence in f ν,ν (k, k ) can be omitted. Under the external magnetic field, the semi-classical equations of motion of the quasiparticle can be written as [49,50] where v ν (k) = ∂ k 0 ν (k) is the band velocity, and ν (k, r, t) = δE total /δn ν (k, r, t) is the renormalized quasiparticle dispersion given by the functional derivative of the total energy. Then the Boltzmann's equation is given by where the first term in the right hand side describes the evolution of distribution function in phase space under external field, and the second term describes the scattering process. An explicit derivation of this equation can be found in Appendix A. Now we decompose δn(k, r, t) into two parts: the part that preserves quasiparticle number of each valley, δn ν (k, r, t), and the part that changes valley particle number δn ν (k, r, t). In the following we refer δn ν (k, r, t) as the Fermi surface degree of freedom, and δn ν (k, r, t) as the valley degree of freedom. Since the intra-valley scattering keeps particle number in each valley invariant, δn ν (k, r, t) can be relaxed only through the inter-valley scattering. On the other hand, δn ν (k, r, t) can be relaxed through both the inter-and intra-valley scatterings. Therefore, δn ν (k, r, t) and δn ν (k, r, t) in general have different relaxation time and the relaxation of the former is always much slower than that of the latter. In general we can approximate the scattering term as As proved in Appendix E, for the simplest case where both inter-and intra-valley scattering cross sections are constants without k-dependence, i.e., the scattering cross section takes the form Such a model scattering cross section is a good approximation for Weyl semimetals where the Fermi surfaces are very small compared to the total volume of the Brillouin zone. In this case, the δn ν (k, r, t) part takes the form Now we argue that, in the chiral limit, where τ 0 τ v , the Fermi surface degrees and the valley degrees are decoupled, and the collective modes are purely contributed by the valley degrees. To zeroth order of τ 0 , any nonzero δn ν (k, r, t) would be immediately quenched in an infinitely short time, thus the Fermi surface degrees are always in equilibrium, i.e., δn ν (k, r, t) = 0. Therefore, to get the dynamic equation in the chiral limit, we can simply assume δn ν (k, r, t) = δn ν (k, r, t) = c ν e iq(·r−ωt) δ( 0 ν (k)−µ). Here q and ω are the wave-vector and frequency of the corresponding collective modes, respectively. As detailed in Appendix B, by substituting this trial solution into Eq. (3) and integrating over k, we get the following dynamic equation where β ν (B) is the compressibility of the ν-th valley in absence of interactions, χ ν = ±1 is the chirality, and η ν = β ν (B)c ν is the change of quasiparticle density. As discussed in Appendix B, in the semiclassical region there As we only consider the case ω B µ, so there is always To give any credence of Eq. (5), we need to find the upper bound of τ 0 below which the zeroth order discussion is valid. The effect of finite τ 0 is dealt with second-order perturbation theory, as presented in Appendix E. Here we only describe the conclusion: finite τ 0 introduces an effective damping rate ∼ τ 0 v 2 F q 2 for the collective modes. For the collective modes to be stable, the Hermitian part of Eq. (5) must be larger than the non-Hermitian part, or equivalently the eigenfrequency should much larger than the damping. For the gapped CPs, the coulomb interaction dominates in the longwave limit, thus the stable con- 2 , which can be always satisfied in the longwave limit. Therefore, the gapped CPs are stable against τ 0 , which is reasonable because the existence of gapped plasmons does not require the chiral limit. For the CZSs and gapless CPs, the coulomb interaction in Eq. (5) does not contribute (as exampled in the model in the following and proved generically in Appendix B), thus the stable conditions are (i) , which can be satisfied at some q only if For simplicity, here we assume isotropic Fermi surfaces such that β(B) = µ 2 Thus the upper bound of τ 0 for CZSs and gapless CPs to be stable is given by Eq. (6). (See Appendix E for a more general discussion.) Before going to particular models, let us analyze the symmetry of Eq. (5). It should be noticed that the wavevector q enters Eq. (5) only through the q · B term. As a consequence, the symmetry of Eq. (5) is much higher than the little group at q: all the symmetries that preserve q · B and B are kept. In this paper we denote this (magnetic) point group as G c (B). Notice that q · B is a pseudo-scalar, which is invariant under proper rotation and flips sign under inversion and is invariant under time-reversal, and B is a pseudo-vector, which rotates as a vector under proper rotation and is invariant under inversion and flips sign under time-reversal. Therefore, G c (B) consists of proper rotations with axis parallel with B, and time-reversal followed by proper rotations with axis perpendicular with B. On one hand, as discussed at the beginning of this paper and illustrated in the following example, the CPs form the trivial irrep of G c (B), whereas the CZSs form the nontrivial irreps of G c (B). On the other hand, the representation formed by η ν is completely determined by the configuration of WPs. Thus the number of CPs and CZSs can be determined in prior to solving the dynamic equation. As proved in Appendix B, the number of CPs are given by and the number of CZSs are given by Here ν sums over all inequivalent WPs, G 0 is the maximal (magnetic) point group of the space group, G ν is the subgroup of G 0 that leaves the ν-th WP invariant, and |G| represents the number of elements in G. We emphasize that even there can be more than two branches of CPs, only two of them are fully gapped, and others are gapless in the longwave limit but are coupled to the gapped CPs at finite momentum. (See Appendix D for more details.) Here we take TaAs [5] as an example to show the usage of Eq. (7) and (8). TaAs has a space group I4 1 md and time-reversal symmetry, and the maximal (magnetic) point group is G 0 = C 4v + T C 4v . Totally there are 24 different WPs in TaAs, which can be divided into two classes, 8 WPs located at the k z = 0 plane and 16 WPs located off the k z = 0 plane. The WPs within the same class can be related by point group operators and are considered to be equivalent from symmetry point of view. Then the corresponding little groups that leave the WPs unchanged are G 1 = {E} and G 2 = {E, T C 2 }, respectively. Therefore, from equation (5), there are totally 24 independent variables leading to same number of independent modes. Assuming the magnetic field is applied along the C 4 rotation axis, we get G c (B) = C 4 . Therefore, we get N CP = 6 and N CZS = 18 for TaAs.
As discussed above and detailed in Appendix B, in the semiclassical region we always have . Thus, to leading order effect of magnetic field, we can omit the B-dependence in β ν (B). Then Eq. (5) is in first order of B and the corresponding symmetry group becomes higher than G c (B). This higher symmetry group, denoted as G c (0), consists of all the proper rotations, time-reversal (if present), and time-reversal followed by proper rotations (if present) in the original group. Thus G c (0) is nothing but the chiral subgroup of the little group at q = 0. In this approximation, the number of CPs and CZSs (Eq. (7) and (8)) should be calculated with G c (0) instead of G c (B). Now let us consider a model Weyl semimetal system with only two pairs of WPs. The system has C 2v symmetry and a square lattice, as illustrated schematically in Fig. 1. For convenience, here we split the valley index ν as a chirality index s = ±1 and a sub-valley index a = 1, 2 such that the chirality χ sa = s. Then it is direct to write down the representation matrices formed by η ν for the C 2v group as D sa;s a (C 2 ) = τ x a,a σ 0 s,s and D sa;s a (M x ) = τ 0 a,a σ x s,s . Here C 2 and M x represent the 2-fold rotation and the vertical mirror shown in Fig. 1, respectively, and τ 0,x,y,z and σ 0,x,y,z are Pauli-matrices in the chirality space and sub-valley space, respectively. In the following, we would omit the matrix subscripts for brevity. Without loss of generality, here we fix the magnetic field in the z-direction and set β sa (B) as β(B) and to ensure that the short range interaction is positive semi-definite. As discussed above, only the proper rotations are preserved in G c (0). For this model, one can easily verify that the dynamic equation preserves the C 2 symmetry but breaks the M x symmetry. Thus the modes should form the irreps of the point group C 2 . Solving Eq. (5), we get two branches of CPs (9) where ζ 0 (q) = 1+β(B)(f 0 +f 1 +2e 2 /( 0 q 2 )) and ζ 1 (q) = β(B)(f 2 + f 3 + 2e 2 /( 0 q 2 )), and two branches of CZSs where , so the CP modes are gapped and the plasmon frequency is approximately However, the CZS modes have a linear dispersion along the magnetic field direction, with a sound velocity Here we give a rough estimation of the sound velocity for a typical Weyl semimetal system. For simplicity, we where λ 1,2 (q) = (ζ 0 (q) ± ζ 2 0 (q) − ζ 2 1 (q))/ζ 1 (q), and for the CZSs, the eigenvectors are where Here the subscript in η is ordered as (s, a) = (+, 1), (−, 1), (+, 2), (−, 2). Apparently, η (1,2) are invariant under C 2 , so they form the trivial irrep, whereas η (3,4) change sign under C 2 , and so each of η (3,4) forms a nontrivial irrep. The oscillation pattern for both the CP and CZS modes are schematically plotted in Fig. 1 (a) and (c) together with the corresponding contribution to the charge current from the four different WPs in Fig. 1 (b) and (d). We can find clearly from Fig. 1 that CP is such a mode that the quasiparticle densities of Weyl fermions with the same chirality oscillate with the same phase, while those of Weyl fermions with the opposite chirality oscillate with the opposite phase. Therefore, net current will be generated during the oscillation of CP, which couples the motion to the long-range Coulomb interaction and makes it a special plasmon mode in the chiral limit. In contrast, for the CZS modes, the valley density with the same chirality oscillates with the opposite phase, which cancels out the net charge current and makes it decouples to the charge dynamics completely. We would emphasize that both the CP and CZS are entirely different with the normal plasmon and ZS modes because all the normal collective modes here will be damped out in the chiral limit. Moreover, both the CP and CZS can only propagate along the magnetic field. It is insightful to compare the possibility to have zero sound modes in ordinary metals and Weyl semimetals under magnetic field. The collective modes for the former have been discussed in detail in Ref. [41]. Using the description developed above, for a ordinary metal, all the collective modes can be derived from the dynamics of Fermi surface degree of freedom δn ν (r, k, t), which describes the small deviation of the quasi-particle occupation at the Fermi surface. For a system with approximately the sphere symmetry, it can be expanded using the sphere harmonic Y lm (θ k , φ k ). Therefore, the longitudinal mode is formed by the proper linear combination of the sphere harmonics with m = 0 and becomes plasmon mode. The transverse modes are described by the sphere harmonics with m = 0. Among them, the channel with m = ±1 will be absorbed into the Maxwell equation to describe the possible electrical magnetic wave, which contains no solution for frequency less than plasmon edge. Therefore, the only possible channels to have zero sound modes in a ordinary metal system is the channels with |m| ≥ 2 provided that the effective residual interaction in these channels are positively defined to survive the Landau damping. These conditions are difficult to fulfil and so does the zero sound in ordinary metal. There-fore, for Weyl semimetals under magnetic field, the CME provides a unique mechanism to stabilise CZS with any form of residual interaction that doesn't cause instability. At least in the chiral limit, the dynamics of CZS only involves the anomalous current but not the normal current, it is free of Landau damping.
The existence of CP and CZS in the chiral limit can lead to several interesting physical phenomena under a magnetic field. Here we introduce two of them. The first one is the CZS contribution to the specific heat. The CZS can be view as a set of 1D collective modes dispersing only along the magnetic field. As derived in Appendix F, the specific heat contributed by CZS is proportional to T for temperature T Θ CZS , where Θ CZS = c(B)Λ/k B is the corresponding Debye temperature for the CZS. Here c(B) is the sound velocity (Eq. (12)), Λ is the momentum space cutoff, and k B is the Boltzmann's constant. Θ CZS is proportional to the field strength in the weak field case as the velocity of CZS. Although such linear temperature dependence is similar with the quasiparticle contribution to the specific heat, the former can be distinguished from the latter by its unique field dependence. To be specific, in Fig. 2(a) we plot the specific heat as a function of temperature using some typical parameters for the Weyl semimetal systems. As calculated in Appendix F, the specific heat in low temperature region is given by κ(B) ≈ k 2 B T Λ 2 /(6c(B)), and in high temperature region (T Θ CZS ) is given by κ(T ) ≈ k B Λ 3 /(3π 2 ).
Another unusual property caused by CZS is the thermal conductivity. Since the CZS disperses only along the field direction, the thermal current carried by the CZS modes can also only flow in this direction. As a result, the thermal conductivity tensor contributed by CZS modes has only one nonzero entry. As derived in detail in Appendix F, if the magnetic field is applied along the z-direction, this thermal conductivity is given by σ th ij = δ i,z δ j,z τ s (T )c 2 (B)κ(B), where τ s (T ) is the relaxation time for CZS excitations and c(B) is the velocity of the CZS defined in Eq. (12). Since the dynamic equation (Eq. (5)) that leads to CZS depends on the compressibility, which oscillates with the field, if the field is strong enough, as a consequence the velocity of the CZS, as well as the thermal conductivity in general, should also oscillate with the field. To discuss the quantum oscillation of the CZS velocity and thermal conductivity, we need to re-derive the dynamic equation and solve it under the strong field, where the electronic states are already Landau levels. Here we only focus on the case ω B µ so that there are still lots of Landau levels below the chemical potential. As detailed in Appendix G, it turns out that the dynamic equation has the same form of Eq. (5), except that the field dependence of the compressibility is modified. As calculated in Appendix H, the compressibility in strong feild can be expressed as β(B) = β (0) (B) + β (1) (B) + · · · , where the β (l≥1) (B) terms oscillate as the lth harmonics of 1/B. In the presence of finite temperature, the ratio between the first order and zeroth order components is approximately where sinch(x) = (e x − e −x )/(2x), and S ex (µ) is the area enclosed by the extreme circle on the Fermi surface that is perpendicular to B. In deriving this equation, we have assumed eB > 0 and µ > 0. Due to Eq. (5) and (12), the oscillation in β(B) will lead to oscillation in the sound velocity in CZS. Substiting Eq. (15) to Eq. (12), we get the first order oscillation in sound velocity as where c (0) (B) is the non-oscillating component of the sound velocity. As both the specific heat and thermal conductivity are functions of sound velocity, the oscillation in velocity leads to oscillations in specific heat and thermal conductivity as well. As an example, in Fig. 2(b) we plot the thermal conductivity as a function of magnetic field. In normal metals, the thermal conductivity is mainly contributed by electrons and acoustic phonons. The phonon part only couples indirectly to the magnetic field and usually doesn't change much with the field. Therefore the part that oscillates with the field is mainly contributed by the free electrons in the normal metal, which satisfies the Wiedemann-Frantz law. While as we have introduced above, for the Weyl semimetals in the chiral limit, since the CZS can only propagate along the magnetic field, the thermal conductivity along the field will be contributed by both the CZS and free electrons leading to the dramatic violation of Wiedemann-Frantz law, which is absent for thermal conductivity along the perpendicular direction. Such field dependent violation of the Wiedemann-Frantz law had already been seen in the thermal conductivity measurement of TaAs under a magnetic field, indicating the possible contribution from the CZS. We would note that for realistic systems, which are not deeply in the chiral limit, the CZS will also acquire nonzero velocity along the transverse direction of the magnetic field as well, which is caused by the accompany normal current during the oscillation. Therefore, the CZS or the gapless CP can also contribute to the thermal conductivity along the transverse direction but the effect should be much less than the longitudinal direction.
The above mentioned quantum oscillations in specific heat and thermal conductivity can be viewed as strong evidence for the existence of CZS but still indirect. It will be much convincing, if we can also have direct ways to measure it. On this regard, the direct ultrasonic measurement of these materials under magnetic field and low temperature may be difficult but worth trying. Another possible experiment is inelastic neutron scattering spectrum. Although the corresponding scattering cross section for electrons may be very small, the existence of CZS can still be inferred from the spectrum of certain phonon modes, which have the same symmetry representation with the CZS and can hybridise with it at some particular wave vector to form the "polariton mode".
In summary, we have proposed that an exotic collective mode, chiral zero sound, can exist in a Weyl semimetal under magnetic field with the chiral limit, where the inter-valley scattering time is much longer than the intravalley one. The CZS can propagate along the external magnetic field with its velocity being proportional to the field strength in the weak field limit and oscillating in the strong field. The CZS can lead to several interesting phenomena, among which the giant quantum oscillation in thermal conductivity is the most striking and can be viewed as the "smoking gun" evidence for the existence of it.

Appendix A: Boltzmann's equation and collective modes
Let us first derive the Boltzmann's equation, which applies when the Landau level splitting, i.e., ω B = v F √ eB, is smaller than the imaginary part of the quasiparticle self-energy and the chemical potential, µ. The semiclassical equations of motion of Weyl fermion are [49,50] where e = −|e| (|e|) for electron-like (hole-like) quasiparticle, and is the Berry's curvature. The decoupled equations are is the phase space measure. Now we denote the distribution function over phase space as ρ(k, r, t), due to particle number conservation, we have and hence Here we have neglected the scattering term in Boltzmann's equation.
From now on, we assume there are a few valleys and label quantities in different valleys with a subscript ν. For each valley, we introduce a weighted distribution function n ν (k, r, t) = ρ ν (k, r, t)/γ(k), then the multi-valley Boltzmann's equation is given by where, again, the scattering is neglected. In derving Eq. (A9) we have made use of the relations ∂ r · (γ ν (k, B)ṙ ν ) = 0 and ∂ k · (γ ν (k, B)k ν ) = 0. Due to Eq. (A4) and (A5), this two relations are satisfied as long as (i) k is not at the Weyl point, where the semiclassical method does not apply, and (ii) ∂ r · ∂ k (k, r, t) = 0, which is automatically satisfied in our approximation for quasiparticle energy (Eq. (A13)).
In presence of collective mode, the single particle energy ν (k, r, t) should be determined self-consistently. With quasiparticle excitation, the total energy is a functional of the distribution function [47] δn(r, t)δn(r , t) where the second and third terms denote the long range Coulomb and residual short range interaction among the quasiparticles around the WPs respectively. Here is the deviation of distribution from equilibrium, n F ( ) = 1/ 1 + exp − kBT is the Fermi-Dirac distribution, is the charge density at position r and time t, f ν,ν is a real matrix due to the Hermitian condition of Hamiltonian.
The k-dependence of f ν,ν is neglected since we consider the case where the Fermi surfaces are very small compared to the Brillouin zone. The quasiparticle energy can then be derived as the functional derivation of the total energy where ϕ is the scalar potential determined by possion equation Now we assume the deviation from equilibrium takes the form of plane wave δn ν (k, r, t) = δn ν (k)e i(q·r−ωt) .
Following this definition, we can rewrite the quasiparticle energy as The equation of motion to first order of δn ν (k) is given by where v ν (k) = ∂ k 0 ν (k), and δ T ( ) = −∂ n F ( ). For convenience, we replace the 3D variable k in Eq. (A17) with an energy and a 2D wave vector σ on the energy surface. The integration over k in the ν-th valley can be rewritten aŝ where´ means σ takes value on the 2D surface with fixed energy . Apparently, the solution of Eq. (A17) takes the form where σ takes value on the Fermi surface. Integrating the energy, Eq. (A17) becomes wherev ν (σ) = v ν (σ)/|v ν (σ)|, and In Eq. (A20) all the quantities are defined on the Fermi surfaces, so we omit the energy dependence of these quantities, e.g., v ν (σ) is a shorthand of v ν (µ, σ).

Appendix B: The chiral limit
We can decompose δn as a part changing particle number in each valley, δn, and a part deforming the shape of Fermi surface but preserving particle number in each valley, δn. We refer δn as the valley degree of freedoms and δn as the Fermi surface degree of freedom. In general case these two degrees of freedoms are strongly coupled. However, as argued below, in the chiral limit, the dynamic of these two degrees of freedoms are decoupled. In presence of scattering term, δn in general damps with time, but the valley degrees and the Fermi surface degrees can have different relaxation time. We denote the relaxation time of δn as τ 0 whereas the relaxation time of δn as τ v . Then the time derivative term in Eq. (A20) should be replaced by The chiral limit refers to the case that τ 0 is much smaller τ v , i.e., This limit can be achieved when the intra-valley scattering is much stronger than the inter-valley scattering. In Appendix E we discuss the relaxation times contributed by impurity scattering. In the simple case in Appendix E, δn ν is defined as where is the compressibility of the ν-th valley, and ρ ν is total particle density of the ν-th valley. Then Eq. (A20) can be rewritten as In the following, we study the physics in zeroth order of τ 0 , and leave the discussion on finite τ 0 effect in Appendix E. To zeroth order of τ 0 , Fermi surface degrees of freedom is always in thermal equilibrium, i.e., δn ν (σ) = 0: any deviation from equilibrium will be immediately killed by the strong scattering. By integrating σ in Eq. (B5), we get a generalized eigenvalue equation Here χ ν = ±1 is the chirality of the ν-th valley, and η ν = β ν (B)δn ν is the inequilibrium quasiparticle number in the ν-th valley. In deriving Eq. (B6), we have applied Now let us discuss the symmetry of Eq. (B6). Apparently, Eq. (B6) has a higher symmetry than the little group group of q: it contains all the symmetries that preserve the chiralities of WPs and the direction of magnetic field. The direction of q is irrelevant to the symmetry. This is because, in the chiral limit, the electric field, proportional to q, enters the equation only through the q · B term, and thus only couples to the chiral degree of freedom. Therefore, finite q only breaks the symmetries changing chiralities. In the following, we denote the symmetry group of Eq. (B6) as G c (B). We emphasize that some anti-unitary symmetry, like time-reversal followed by a crystalline symmetry, can also keep the chiralities and the magnetic field invariant. And, since f ν,ν is a real matrix, these anti-unitary symmetries act on Eq. (B6) as unitary operators. The explicit representation matrix of all these symmetries is given in Eq. (B9).
It should be noticed that, to leading order of magnetic field, i.e., setting β ν (B) = β ν (0), Eq. (B6) even has a symmetry higher than G c (B): the magnetic field enters the equation only through term q · B, thus the direction of magnetic field becomes irrelevant to the symmetry. We denote this higher symmetry group as G c (0), which consists of all the symmetries preserves the chiralities of WPs.
Solutions of Eq. (B6) must form irreducible representations (irreps) of G c (B). As will be shown in next two sections, the trivial irreps of G c (B) always couple to the charge density oscillation, and thus we call the modes forming trivial irreps as chiral plasmons (CPs). As will be proved, only two of the CPs are gapped, whereas other CPs are gapless in the longwave limit. On the other hand, all the nontrivial irreps are decoupled from density oscillation, so we call them the chiral zero sounds (CZSs). Now let us calculate the number of trivial irreps in the solution of Eq. (B6). We first consider a set of symmetry related WPs in the inner of Brillouin zone, and one of them has the little group G W . We denote the maximal (magnetic) point group of the space group as G 0 , then each symmetry-related WP can be represented by a coset representative of G 0 /G W The representation formed by the valley degrees is given by Now we reduce D to irreps of G c (B). The number of trivial irrep is given by Therefore, for a system with a few set of non-equivalent WPs, the number of CP modes and CZS modes are given by and respectively. Here ν sums over all inequivalent WPs, G ν is the little group of the ν-th WP.

Appendix C: Chiral zero sound (CZS)
If η ν is not a trivial irrep of G c (B), there must be ν η ν = 0, implying that it does not cause any charge density oscillation. Thus for nontrivial irreps we can omit the Coulomb term, and the corresponding modes are the CZSs. Now let us solve the equation of motion for CZS. Notice that the matrix in the right hand side of Eq. (B6) is real and symmetric, so we diagonalize it as where β(B) is the averaged β ν (B), O is an orthogonal matrix, and λ a 's are dimensionless numbers. Applying the transformation we can rewrite Eq. (B6) as Applying the transformation we get a regular eigenvalue problem where Ξ has the symmetry of G c (B). The dispersion of CZS is given by Here ξ n is the n-th eigenvalue of Ξ. It should be noticed that Ξ is real and symmetric (such that ξ n 's are real) only if all λ a 's are positive. Thus the number of CZS modes is given by Eq. (B13) only if Eq. (C1) is positive definite. Otherwise, only irreps where all λ a 's are positive correspond to physically observable modes. The irreps having negative λ a 's in general have complex ξ and so are not stable.

Appendix D: Chiral plasmon (CP)
For the trivial irreps, in general we have ν η ν = 0. Therefore the trivial irreps contribute to density oscillation and thus the Coulomb term must be considered. However, as e 2 0 q 2 is a rank-1 operator, in the longwave limit, there should be only one channel that responses to Coulomb interaction. To separate this channel, we define the projection operator P ν,ν = 1 N W , where N W is the number of WPs, and divide the terms in the right hand side of Eq. (B6) to four components: Here e 2 0q 2 represents the matrix where every element is e 2 0q 2 , 1 β(B) represents the diagonal matrix 1 βν (B) δ ν,ν , and Q = I − P , where I is the identity matrix. We apply an orthogonal transformation V = I + S − S T , where S = P SQ, to remove the mixing term between P and Q subspaces. To second order of q, we find that and Using the fact ν V ν ,ν χ ν V ν ,ν = χ ν δ ν,ν + O(q 2 ), we can rewrite Eq. (B6) as To solve this generalized eigenvalue equation, we apply the technique used in Appendix C: diagonalizing the matrix in the right hand side and transforming the equation to a regular eigenvalue problem. Let us write the matrix in the write hand side as 1 we get a regular eigenvalue problem Ξη = e(q · B) The frequencies of CP modes are then given by where ξ n 's are eigenvalues of Ξ. Now let us analyse the spectrum of Ξ. For convenience, we set O ν,1 as the eigenvector of P , so the corresponding eigenvalue is which is singular in the limit q → 0. Then, due to Eq. (D12), the Ξ matrix takes the form where and Ξ a,a = Ξ a,a , a, a = 2, 3 · · · , are submatrices of Ξ. We emphasize that for a ≥ 2, λ a is not singular. Thus in the limit q → 0, Ξ approaches a constant matrix, whereas ζ 1,a ∼ |q|. Therefore, by diagonalizing Ξ , we can rewrite Ξ as where ξ a 's are eigenvalues of Ξ , and U is some orthogonal matrix. Now we prove that one of ξ a=2,3··· is zero. We denote the diagonal matrix δ ν,ν χ ν as [χ]. Then the projected Apparently, η ν = 1 and η ν = χ ν are two zero eigenvectors of Q[χ]Q, wherein η ν = 1 is in subspace P whereas η ν = χ ν is in subspace Q. As Ξ is equivalent to Q[χ]Q up to an invertable transformation, Ξ has one zero eigenvalue in the Q subspace. Therefore, one of ξ a=2,3··· is zero. Here we choose ξ 2 = 0. In the limit q → 0, we have the Ξ eigenvalues as Therefore, due to Eq. (D18), n = 1, 2 correspond to the gapped CP modes, whereas n = 3 · · · N CP (B) correspond to the gapless CP modes. The low energy behavior of the gapless CPs are very similar with the CZSs: both of them have a linear dispersion in the limit q → 0. However, a vital difference is that gapless CPs are coupled to gapped CPs, through the c a≥3 terms in Eq. (D12), whereas the CZSs cannot. As a result, the dispersions of gapless and gapped CPs form anti-crossings, whereas dispersions of CZSs and gapped CPs form symmetry-protected crossings.

Appendix E: Finite intra-valley scattering
In this section, we solve Eq. (B5) to first order of τ 0 and justify the chiral limit approximation using a second order perturbation theory.
Firstly, let us derive τ 0 and τ v explicitly in terms of scattering cross section. We model the scattering cross section as W ν,ν (σ, σ ) = δ ν,ν w 0 + (1 − δ ν,ν )w 1 . (E1) Thus the scattering term is given by and δn ν (σ) = δn ν (σ) − δn ν , we can rewrite the scattering term as The first term relaxes deformation of Fermi surfaces that does not change quasiparticle number in each valley, the second term relaxes the quasiparticle number in each valley, and the last term is a feedback from the change of total quasiparticle number. For simplicity, in the following we will neglect the ν-dependence in τ 0,ν , i.e., setting τ 0,ν = τ 0 .
For simplicity, here we consider isotropic Fermi surfaces, where |v ν (σ)|, β ν (σ), Ω ν (σ) ·v ν (σ) do not depend on σ, and γ ν (σ, B) = 1. According to Eq. (B5), to leading order of τ 0 we get the δn ν (σ) part as δn ν (σ) = −iτ 0 (q · v ν (σ)) ω + i τv χ ν e 4π 2 (q · B) η ν + O(τ 2 0 ). (E7) Now we look at the leading order effect of τ 0 on CZS modes. For a specific branch of CZS, Eq. (E7) gives where ξ is the corresponding eigenvalue of Ξ matrix (Eq. (C4)). Substiting it back into Eq. (B5) and integrating σ, we get where Apparently, finite τ 0 introduce a non-Hermitian term in the dynamic equation of η. This term would lead to a damping rate proportional to τ 0 q 2 v 2 F /ξ. Therefore, for the zero sound to be stable, the following relation should be satisfied Considering the inter-valley scattering, the following relation should also be satisfied eBq βξ The above two inequalities are equivalent to which have solutions only if Eq. (E14) gives the upper limit of τ 0 , above which the CZS modes become unstable. It should be noticed that 1/ξ is in order of 1 + β(B)|f |. Perturbation theory for gapless CP modes is similar with the perturbation theory for the CZS modes, and the stable condition of gapless CP modes is also given by Eq. (E14). Now we consider the gapped CP modes. For the positive branch of gapped CP, the frequency of which is denoted as ω CP , Eq. (E7) gives Following the above analysis, we find this term would lead to a damping rate proportional to β(B)(qv F )(τ 0 ω CP )/eB ∼ (qv F )(τ 0 ω CP ) µ 2 ω 2

B
. Thus, for the gapped CP to be stabel, there should be Therefore, the CP modes in the longwave limit are always stable against the intra-valley scattering.