{'fraction_non_alphanumeric': 0.04886511738122239, 'fraction_numerical': 0.013760769113589817, 'mean_word_length': 3.744017986189176, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 50, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Phonon helicity and Nieh-Yan Anomaly in the Kramers-Weyl semimetals of Chiral Crystals

Nieh-Yan anomaly describes the non-conservation of chiral charges induced by the coupling between Dirac fermions and torsion fields. Since the torsion field is beyond general relativity, this effect remains hypothetical and its relevance to our universe is unclear in the context of high-energy physics. In this work, we propose that the phonons can induce a torsion field for the Kramers-Weyl fermions through electron-phonon interaction in a non-magnetic chiral crystal, thus leading to the occurrence of the Nieh-Yan anomaly in this condensed matter system. As a consequence, the Nieh-Yan term can strongly influence the phonon dynamics and lead to the helicity of acoustic phonons, namely, two transverse phonon modes mix with each other to form a circular polarization with a non-zero angular momentum at a finite phonon momentum and the phonon angular momentum reverses its sign for opposite momenta due to time reversal symmetry. The phonon helicity can be probed through measuring the total phonon angular momentum driven by a temperature gradient.

In Weyl semimetals, position-dependent and time-dependent perturbations, such as magnetic fluctuations 27 , strain 28,29 and structure inhomogeneity 30,31 , can shift the position of Weyl nodes and thus act as an emergent gauge field, dubbed the "pseudo-gauge field" 32 , which provides a unified description of various physical phenomena in Weyl semimetals. The chiral zero Landau levels induced by pseudo-gauge fields have been observed experimentally in photonic and acoustic Weyl metamaterials 30,31 . Furthermore, these perturbations can also give corrections to the Fermi velocity of Weyl fermions and thus play the role of the frame fields (also called tetrad or vierbein) 33,34 , allowing to mimic Weyl fermions in the curved spacetime 26,[35][36][37][38][39] . An exciting theoretical proposal is to realize the Nieh-Yan anomaly 10-13 due to the torsion field in Weyl semimetals [40][41][42][43][44][45][46][47][48] , which may be probed through topological magnetotorsional effect 49 , anomalous thermal Hall conductance 42 , or the sound-wave-induced current oscillations in tilted Weyl semimetal interfaces under magnetic fields 50 . However, one challenging still remains: since the perturbation terms normally shift the Weyl node location and modify the Fermi velocity simultaneously, both the chiral anomaly and Nieh-Yan anomaly can be induced by these perturbation terms, and thus their roles as the pseudo-gauge fields and the frame fields are difficult to distinguish. Furthermore, it is still unclear how the quasiparticle excitations due to these perturbations, such as elastic wave, phonon dynamics and magnetic dynamics, are influenced by the Nieh-Yan anomaly of Weyl fermions.
To address these challenges, we turn to the so-called "Kramers-Weyl" (KW) fermions in chiral crystals. Chiral crystals are the crystal structures with a well-defined handedness and non-magnetic chiral crystals can generally host the KW fermions when taking into account spinorbit coupling (SOC) [51][52][53][54][55] . The location of KW fermions is pinned at time-reversal (TR) invariant momenta by TR symmetry. Consequently, any TR-preserving perturbations cannot shift the KW fermions in the momentums space and thus will not induce the vector potential of pseudo-gauge field for KW fermions. In this work, we will consider the role of acoustic phonons in the KW semimetal phase of chiral crystals and demonstrate that the phonons will induce a torsion field for the KW fermions. By integrating out the KW fermions, we prove that the phonon self-energy contains an off-diagonal term that originates from the Nieh-Yan anomaly of the torsion fields. While this self-energy correction has no influence on the longitudinal phonon mode, it will mix two transverse modes and give rise to the phonon angular momentum (PAM) at a finite phonon momentum. In particular, the induced PAM reverse its sign for opposite phonon momentum, as shown in Fig. 1a, and form a hedgehog texture in the momentum space for one phonon branch (Fig. 1b), analogous to the helical spin texture in spinorbit coupled electronic band structures. Thus, we term it as "phonon helicity" 56 . The phonon helicity can be probed through measuring the total PAM induced by a temperature gradient 57 .
Electron-phonon Coupling in Chiral Crystals -In nonmagnetic chiral crystals, all the energy bands are at least doubly degenerate (spin degeneracy) at high symmetry momenta, labelled as Γ i (i = 1, 2..., 8), and SOC can lift this spin degeneracy for the momentum away from Γ i , thus giving rise to the KW fermions. If the full ro- tation symmetry is assumed in the low-energy sector of the KW fermions (see Appendix A and B), the low energy physics can be described by an isotropic Hamiltonian H 0 = 2 k 2 2m * 0 + v f k · σ − µ expanded around Γ i up to the k 2 terms, where σ labels the Pauli matrix for spin, m * 0 is the effective mass, v f is the Fermi velocity and µ is the chemical potential. As discussed in Ref [58], this isotropic Hamiltonian remains valid for the chiral point groups T and O, while for other chiral point groups, the Fermi velocity v f becomes anisotropic and should be replaced by a tensor. We only focus on the isotropic Hamiltonian in this work. In chiral crystals, there are normally multiple Fermi surfaces, and here we only focus on the parameter regime where the linear-k terms dominate over the k 2 terms at the Fermi surfaces around different Γ i , as depicted in Fig. 1c. We call this parameter regime as the "KW semimetal phase", in which the k 2 term only limit the valid regime in the momentum space for Weyl physics and thus can be replaced by a momentum cut-off Λ = 2m * 0 v f (assuming µ 2 Λ 2 2m * 0 ).
The electron-phonon coupling in this system can be derived from the standard k · p theory up to the second order perturbation terms with the help of Schrieffer-Wolf transformation, as discussed in Appendix A. The TR symmetryT = iσ y K, where K is the complex conjugate, gives a strong constraint on the form of the resulting effective Hamiltonian. In this work, we focus on the acoustic phonons, which couple to electrons through an internal strain, described by the strain tensor u ij = 1 2 (∂ i u j + ∂ j u i ) with the displacement vector u i (i, j = x, y, z). Since the strain tensor is even underT , it cannot directly couple to either electron momentum k or spin σ, both of which are odd underT . This implies that the strain tensor cannot play the role of the vector potential of pseudo-gauge field. Based on the symmetry construction in Appendix B, the symmetry-allowed strain-electron coupling Hamiltonian for an isotropic KW fermion is written as H ep = (C 1 + g 0 k j σ j )u ii + g 1 u ij k i σ j up to the order of k i σ j , since both u ij and T ij = k i σ j are rank-2 tensors. We have assumed the summation over the duplicated indices. C 1 , g 0 , g 1 are three independent parameters. The C 1 term provides a correction to the chemical potential, while both g 0 and g 1 terms give the corrections to the Fermi velocity tensor. Microscopically, both g 0 and g 1 terms are from the second order perturbation, combining the k · p term and the bare electronphonon coupling term, as shown in Appendix A.
Putting H 0 and H ep together, one can derive the effective action for the KW semimetals around Γ i , given by whereψ Γ i is the field operator for the KW fermions, τ is the imaginary time, A 0 = C 1 u ii , χ = sign( v f ) and e j a = δ j a + ∆ j a with the ∆ field given by where i, j, a = x, y, z. Here the imaginary-time formalism is used for the study of finite temperature effect.
The Fermi velocity v f is absorbed by re-scaling the spatial coordinate, while we still keep track of the chirality χ (the sign of Fermi velocity). The e j a describes the frame field 33,34 and it only involves non-trivial spatial part due to TR symmetry. The effective Lagrangian (1) describes the Weyl fermions in a space with non-trivial frame field but zero spin connection, which is known as the "Weitzenböck spacetime" studied in "teleparallel gravity theory" 59 . Zero spin connection means vanishing Riemann curvature while non-trivial frame field implies non-zero torsion field, given by T a ij = ∂ iē a j − ∂ jē a i , wherē e a i is the coframe field define byē a i e j a = δ j i (See Appendix D). Thus, strain or phonon can induce a torsion field for the KW fermions in chiral crystals.
Stress-stress correlation function and Nieh-Yan anomaly -Our next step is to integrate out Weyl fermions to obtain an effective action for the ∆ field that is directly related to the strain field. To do that, we separate the effective Lagrangian into Here both the A 0 and ∆ fields in L 1 are related to the strain field induced by acoustic phonons and thus treated as perturbations. By integrating out the KW fermions in L ef f (see Appendix C), we find the effective action W [A 0 , ∆] for the A 0 and ∆ fields (Eq. C15 in Appendix C), which provide the corrections to the phonon effective action. Here we focus on the term W N Y for the ∆ field whereq = (q, iν m ) and q = 1 βV q,iνm . The stressstress correlation function Φ ab ij is defined as (4) wherek = (k, iω n ), G 0 = (iω n + µ − χ(k · σ)) −1 and the stress tensor T a i (k, k ) = χ 2 (k i + k i )σ a with i, a = x, y, z. Our next task is to show that Φ ab ij includes a contribution from the Nieh-Yan anomaly. To see that, we consider iν m = 0 and treat both the momentum q and the chemical potential µ as the perturbations. In this limit, we can first expand Φ ab ij (q) ≈ Φ ab ij (0) + (∂ q l Φ ab ij ) q=0 q l up to the linear order in q. Since ∆ is directly proportional to the strain tensor, the zero order term Φ ab ij (0) just provides the corrections to the elastic moduli (See Appendix G). Thus, below we focus on the linear-q term with its coefficient denoted as Φ ab ij,l = (∂ q l Φ ab ij ) q=0 , which can be further expanded as With these approximations, the effective action is after the Fourier transform into the real space, where µ 5 = χµ 2 has the physical meaning of chiral chemical potential. Now we will show that W N Y is a manifestation of the Nieh-Yan anomaly. To see that, let's recall the Nieh-Yan anomaly equation ∂ µ j 5ξ = F 4 η ab ξνλρ T a ξν T b λρ with the axial 4-current j 5ξ for one Dirac fermion, the torsion field T a ξν = ∂ ξē a ν − ∂ νē b ξ , and certain parameter F (See Appendix D for a short review). This anomaly equation can be obtained from j 5ξ = δS N Y /δA 5ξ with the effective action S N Y = F d 4 xη ab ξνλρ A 5ξē a ν ∂ λē b ρ 49 , where A 5ξ is the axial gauge potential. Now let's consider two KW fermions with opposite chiralities at two TR-invariant momenta, which together form one massless Dirac fermion. For this case, the chiral chemical potential µ 5 in Eq. (5) should be given by µ 5 = 1 2 (µ + − µ − ) with ± for the chiralities of two KW fermions. Since µ 5 is identical to A 50 and the ∆ field gives the spatial part ofē a ν , one finds that W N Y takes the exactly same form as the terms with ξ = 0 in S N Y .
Next let's discuss the coefficient F 0 in W N Y , which is found to take form of F 0 = F 0 + F 1 (k B T ) 2 for the finite temperature (See Appendix E). A simple dimension counting suggests that F 0 has the dimension of 1 L 2 where L is the length dimension. Thus, F 0 is not universal and generally depends on the momentum cut-off F 0 ∝ Λ 2 . The value of F 0 depends on the normalization scheme and our calculation with the momentum cut-off Λ gives F 0 = 0 for the lowest order perturbation term. On the other hand, for the temperature-dependent term, (k B T ) 2 exactly carries the dimension of 1 L 2 . Here we follow the convention in high energy physics with e = = v f = 1, so that energy and momentum have the same dimension. Consequently, the parameter F 1 is dimensionless. It was argued that this parameter is universal and reflects the central charge of Dirac fermions 42 . Our calculations gives F 1 = − 1 12 , consistent with the results in literature obtained from different approaches 42,43,45 , and this T 2dependent term is called "thermal Nieh-Yan anomaly".  (6). The length is in unit of the lattice constant a0, so that the momentum is in unit of 1/a0. The energy is in unit of v f /a0 and Φ xy xx is in unit of v f /a 4 0 .
The above discussion is for a single Dirac fermion with µ → 0. In realistic systems, multiple KW fermions can exist at different TR-invariant momenta Γ i near the Fermi energy. Since our formalism here treats each KW fermion separately, we can simply introduce a single chiral chemical potential µ 5 = 1 2 Γi χ Γi µ Γi for different Γ i into Eq. (5). Thus, below we just focus on a single KW fermion. For each KW fermion, the chemical potential µ with respect to the Weyl point can be far away from zero. This requires us to numerically evaluate the full expression (4) for Φ ab ij (q, iω n ). Here we consider the momentum q along the z-direction (q = qê z ) and numerically compute the component Φ xy xx (See Appendix F for the numerical methods), which is shown in Fig. 2. Φ xy xx is found to be purely imaginary and Im(Φ xy xx ) as a function of q is shown in Fig. 2a for different µ. Im(Φ xy xx ) is linearly proportional to q for q 1. Furthermore, Im(Φ xy xx ) vanishes at µ = 0 and its sign reverses for opposite µ. Fig. 2b shows the dependence of Im(Φ xy xx ) on µ for different q. Im(Φ xy xx ) depends on µ asymmetrically, and for µ close to zero, one can see that Im(Φ xy xx ) is flat and almost zero. This suggests that the linear µ term in Im(Φ xy xx ) should almost vanish, consistent with our analytical result F 0 = 0. Fig. 2c shows Im(Φ xy xx ) as a function of µ for different temperatures k B T . With increasing temperature, the linear-µ term in Im(Φ xy xx ) gradually appears. Finally, we plot Im(Φ xy xx ) as a function of k B T for different µ in Fig. 2d, which shows a T 2 -dependence, as required by thermal Nieh-Yan anomaly. At small µ and low k B T , this temperature dependence can be well captured by the expression from the perturbation expansion (See Appendix F), as shown by the blue curve in Fig.2d for µ = 0.05. When increasing µ, a derivation from Eq. (6) is found for the red, black and green curves in Fig.2d for µ = −0.2, 0.2, 0.4, respectively, implying that the terms with higher orders of µ come into play.
Helicity of acoustic phonons -Since the ∆ field in the effective action W N Y is determined by the strain field u in Eq.
(2), one may naturally expect this term will influence the acoustic phonon dynamics (elastic wave). As discussed in Appendix G, this term will provide a correction d 3 rdtξ ijklm (∂ i u jk )u lm to the effective action of acoustic phonons. In a uniform isotropic system with TR symmetry, there is only a single independent parameter, labelled by ξ 0 , for the rank-5 tensor ξ ijklm . Thus, the effective phonon action is taken the form where ξ 0 = µ 5 F 0 g 2 1 . Here we have changed from imaginary time to real time. This term gives rise to the correction to the equation of motion for acoustic phonons as where u is the displacement field, and c t(l) is the velocity of transverse (longitudinal) modes. We can decompose u into the longitudinal and transverse components, u = u l + u t , and this new term only contributes to the transverse part and has no influence on the longitudinal component. With the wave solution ansatz u t = u 0 e i(q·r−ωt) , the equation of motion of the transverse modes is changed to ω 2 u 0 = c 2 t q 2 u 0 + ξq 2 2 iq × u 0 , which possess the eigen-frequencies ω t s = c 2 t q 2 + s |ξ0| 2 q 3 for the s-mode (s = ±). More importantly, these two transverse modes carry angular momentum l s (q), defined as l s,i (q) = u † 0,s (q)M i u 0,s (q) 57,60-63 , where u 0,s (q) is the polarization vector for the mode s = ± and the M matrix is given by (M i ) jk = (−i) ijk with i, j, k = x, y, z. Direct calculation gives l s (q) = s q q , which satisfies the relation l s (q) = −l s (−q) from the TR symmetry. Thus, the PAM l s forms a hedgehog texture in the momentum space, as shown in Fig. 1b, resembling the helical spin texture of spin-orbit coupled electronic systems. Thus, we term this PAM texture l s (q) as "phonon helicity" 56 .
The total PAM is defined as I ph = s,q l s (q)(f (ω t s ) + 1 2 ), where f (ω t s ) describes the distribution function for phonons. At the equilibrium state, the total PAM vanishes due to the TR symmetry. The phonon helicity l s (q) can give rise to a total phonon angular momentum I ph by driving a thermal current with a temperature gradient 57 , as shown in Fig. 3a. This is in analog to the Edelstein effect of electrons due to spin-orbit coupling 64 . In the linear response regime, the response tensor α ij , defined as I ph i = α ij ∂T ∂xj , can be evaluated from the formula , whereτ is the phonon relaxation time, v ph s = ∂ω t s /∂q is the group velocity of s-phonon mode and f 0 (ω t s ) is the Bose-Einstein distribution function. Since our model is isotropic, the non-zero component α xx = α yy = α zz is calculated as Thus, by plotting α zz /T 3 as a function of T 2 , our theory predicts a linear line dependence, as shown in Fig.3b, and the intercept and the slope of the curve determine the contribution of the normal and thermal Nieh-Yan anomaly. Furthermore, we notice that the chiral chemical potential µ 5 will determine the sign of α zz , which provides a testable signature in experiments. Discussion and Conclusion -To summarize, we have demonstrated that acoustic phonons can induce a torsion field for the KW fermions and consequently, the Nieh-Yan anomaly of the KW fermions can lead to the phonon helicity, which can be probed through measuring phonon angular momentum driven by a temperature gradient. Phonon dynamics in Weyl/Dirac semimetals was previously studied in the context of chiral anomaly [65][66][67][68][69][70][71] and Kohn anomaly 72,73 . The concept of PAM has recently been explored in the studies of the conversion between electron spin and mechanical rotation 60 , selection rules for phonon-assisted electronic inter-valley scattering 63,74 , phonon-induced spin relaxation [75][76][77] , phonon-magnon interaction 78,79 , thermal Hall effect 80 , phonon rotoelectric effect 61 and phonon orbital magnetic moments 81 , and experimental evidence of PAM has been found 62,78 . The PAM due to temperature gradient can be measured through the rigid-body rotation or the phonon orbital magnetic moments, and the theoretical estimate for several realistic materials suggests these phenomena are detectable in experiments 57 . Temperature gradient can also drive electron spin polarization in spin-orbit-coupled systems [82][83][84] , which shows a different temperature dependence from Eq. (9) for PAM. In Appendix H, we also demonstrate that Coulomb interaction will not affect the Nieh-Yan term for the relevant phonon modes. A large number of chiral crystal materials for the KW fermions have been proposed in Ref 51, and to identify the best candidate materials requires more sophisticated first principles calculations of both electronic band structure and phonon spectrum in these compounds, which will be left for the future study.
Appendix A: k · p Hamiltonian and electron-phonon interaction of the Kramers-Weyl fermions In this section, we will describe the k · p theory for the KW fermions in chiral crystals. To be general, we start from the Schrödinger equation where m 0 is the electron mass,p = −i ∇ and the periodic potentialV (r + R n ) =V (r) including both the lattice potential and spin-orbit coupling. The Bloch theorem requires Ψ(r) = e iK·r u(r) where u(r + R n ) = u(r) and K is the crystal momentum. Let's consider the perturbation expansion around a time-reversal (TR) invariant momentum Γ i (i = 1, 2, ... labels different TR invariant momenta) and assume that the Schrödinger equation (A1) has been solved at Γ i with the eigenenergy E Γi,n and eigen-wave function ψ Γi,n (r) = e iΓ i ·r u Γi,n (r). We are interested in the momentum K = Γ i + k where k is considered as a perturbation, and expand the wave function u K (r) = n C n u Γ i (r). Thus, ψ K (r) = e iK·r u K (r) = e i(Γ i +k)·r n C n u Γ i (r) = e ik·r n C n ψ Γ i ,n (r). The Schrödinger equation at K = Γ i + k can be expanded as and P mn = m0 ψ Γ i ,m |p|ψ Γ i ,n , and the corresponding Hamiltonian can be written as Since we focus on TR symmetry, we denote the quantum number m = (α, s), where s labels different spin states, and α labels other band indices. With the TR operatorT , we define the bands of a Kramers' pair asT ψ Γi,α,s = sψ Γi,α,s and E Γi,α,s = E Γi,α,s , wheres = −s. We may assume around Γ i , only one Kramers' pair of bands, denoted as ψ Γi,0,s (α = 0), contribute to the low energy physics and thus the low energy effective Hamiltonian can be written as with s, s = ±. Here we only keep the correction to the electron mass (m 0 → m * 0 ) for the terms of the second order in k. In general, the effective mass m * 0 should be a tensor, but we only consider the quadratic term to provide a cut-off for the valid regime of the KW fermion physics and thus simply treat m * 0 as a scalar. We can decompose the parameter matrix P 0 as P j 0,ss = a P j 0,a (σ a ) ss where σ a=x,y,z are three Pauli matrices. Here TR symmetry requires that the linear k term cannot couples to identity matrix. With T r(σ a σ b ) = 2δ a b , we have P j 0,a = 1 2 T r(P j 0 σ a ). Correspondingly, the effective Hamiltonian can be written as The parameters P j 0,a will be constrained by the crystal symmetry of the system. Since the chiral crystals only involve rotation symmetry, we consider the full rotation symmetry here. As discussed in B, one can show that up to a unitary transformation, the effective Hamiltonian can be written as where the Fermi velocity is given by v f = 2m0 T r(P x 0 σ x ) = 2m0 T r(P y 0 σ y ) = 2m0 T r(P z 0 σ z ) (with full rotation symmetry) and its sign determines the chirality of the KW fermions. It should be mentioned that the above Hamiltonian is also valid for the crystals with the O and T chiral point groups 58 . In the second quantization language, we can In the real space, the Hamiltonian is written as The Fourier transform is given byψ with V the volume of the system, where we have suppressed the lower indices.

Deformation potential for electron-acoustic phonon interaction
The electron-phonon coupling for acoustic phonons should vanish in the long-wave length limit with the phonon momentum q → 0, while there is no such constraint for optical phonons. Thus, our starting deformation Hamiltonians for acoustic phonons and optical phonons are slightly different. For the acoustic phonons, we will start from the continuous model with the electron field described byψ Γi (r) while the acoustic phonon modes can be described by the displacement field u(r). In the case of more than one atoms in one unit cell, the acoustic phonon modes correspond to the polarization vector with the same amount of shift in the same direction for all the atoms in one unit cell and thus we can still use one vector field. The displacement field u can give rise to the polarization P which couples to electron density through its divergence ∇ · P. Let's assume this coupling is local as the lowest order perturbation. Thus, our starting point will be the volume deformation Hamiltonian with the coupling constant g a . Here the spin index s does not appear because volume deformation potential does not involve spin. ∇ · u = i u ii is the trace part of the strain tensor (volume dilation), which is defined asū below. In the momentum space, we have Here we only consider the lowest order bands, and thus only the volume dilation enters into the effective Hamiltonian. If we consider all the other bands, other components of the strain tensor can also appear to couple different bands and the most general form of the electron-acoustic phonon coupling is given by where the form of the coupling parameter g ij a,αβ can be constrained by the crystal symmetry, which will be discussed in Sec. B.

Deformation potential for electron-optical phonon interaction
The optical phonons involve the relative motions of atoms within one unit cell and thus we cannot directly start from the continuous model. Let's next consider the deformation potential where R nτ = R n + r τ with the lattice vector R n and the position r τ for the atom τ in one unit cell. Here the electron field operator can be expanded asΨ(r) = 1 V αsK e iK·r u αsK (r)ĉ αsK , where s labels spin and α labels other band index, while u nτ is the displacement field, which can be expanded as u nτ = 1 √ N Mτ λq Q λ,q e λ τ e iq·Rn . Here M τ is the atom mass, Q λ,q labels the normal mode of the λ-phonon and the polarization vector e λ τ satisfies the equation of motion where i, j = x, y, z and D k (τ i; ξj) is the dynamical matrix for phonons in the momentum space. Since the displacement u nτ is real, Q λ,q = Q † λ,−q . By substituting the expansion ofΨ(r) and u nτ into Eq. (A15) and after some straightforward simplifications, we obtain where We can again expand g λ αs,βt (K, q) in terms of the 2 by 2 matrices as g λ αs,βt (K, q) = µ g λ αβ,µ (K, q)(σ µ ) st where µ = 0, x, y, z. The electron-phonon coupling Hamiltonian is written as The hermitian property of the Hamiltonian H † ep = H ep requires (g λ αβ,µ (K, q)) * = g λ βα,µ (K, −q). The time reversal symmetry provides a constraint on the form of g λ αs,βt (K, q).

Second order perturbation terms from Schrieffer-Wolf transformation
We next consider the second order perturbation contribution to the electron-phonon interaction through the Schrieffer-Wolf transformation. We consider the full k · p Hamiltonian including the electron-phonon interaction around Γ i with the form H Γ i = H 0 + H 1 with the perturbation Hamiltonian H 1 = H k·p + H ep , where H k·p is given by Eq. (A5) with m = (α, s), n = (β, t) and H ep is given by the Eq. (A14) for acoustic phonons and the Eq. (A17) for optical phonons, in which K should be replaced by k since all the momenta are expanded around Γ i .
To get the second order perturbation terms, we consider the Schrieffer-Wolf transformation and we require the first order term to vanish, H 1 − [S, H 0 ] = 0. The corresponding second order term is given by Since and Here we choose optical phonons as the example and keep the lowest order terms, so we can drop the index K, q for the coupling constant g λ αs,βt . For acoustic phonons, we just need to replace Q λ (q) byū ij (q). With the decomposition of S, the effective Hamiltonian includes The first perturbation term is given by which provides a correction to the effective mass. The second and the third term together give is of order g 2 and we neglect this term. Now we apply the above second order perturbation formalism to the Kramers-Weyl fermions to project into the lowest-energy subspace with α = β = 0, and obtain the effective Hamiltonian where and Here we have made the approximation of isotropic effective mass where ... k means the average over all different directions of the momentum k.
and H (2) ef f,ep can be re-written as Transforming back to the real space, we find that H ef f,0 is given by Eq. (A9), H ef f,ep is given by Eq. (A23) and Now let's impose the TR symmetry on the above expression and we find C λ,j should be real. This means that for the C λ,j 0 term, the Hamiltonian will only depend on ∂ j Q λ (r), while the C λ,j a terms (a = x, y, z) depend on Q λ (r). In the long wavelength limit, the C λ,j 0 term will be much smaller than the C λ,j a terms. Thus, we only consider the C λ,j a terms below. The Hamiltonians (A9), (A23) and (A39) together form the starting point for our study on the electron-phonon coupling in the Kramers-Weyl semimetals, and we can rewrite them in a more compact form where the chemical potential µ has included the energy E Γi,0 , A 0 = g λ o Q λ (r) is from the zero order term of electronphonon coupling and behaves as a chemical potential fluctuation (or scalar potential), and the frame field e j a = δ j a + ∆ j a (r) with ∆ j a (r) = 2 v f λ C λ,j a Q λ (r) comes from the second order perturbation involving electron-phonon coupling. In deriving the above Hamiltonian, we only used the TR symmetry and the spatial rotation symmetry will give additional constraints on the parameters C λ,j a and reduce the number of independent parameters, as discussed in the Sec. B. From the above form of the Hamiltonian, it is clear that the phonon provides a background frame field of the curved space for the KW fermions.
The above derivation is directly applied to optical phonons and for acoustic phonons (elastic wave), we just need to change the normal mode Q λ (r) to the strain tensor u ij = 1 2 (∂ i u j + ∂ j u i ), as discussed above. In the next section, we will construct the Hamiltonian coupled to the strain tensor by analyzing the symmetry of the Hamiltonian.

Appendix B: Symmetry construction of the Hamiltonian
In this section, we will consider the construction of the effective Hamiltonian for the KW fermions from the symmetry principle 85,86 . Here we will focus on the acoustic phonons, which couple to electrons through the strain tensor u ij . Similar discussion can also be applied to the optical phonons, which is described by the normal modes Q λ .
The chiral crystals only involve rotation symmetries and to simplify the problem, let's first consider the isotropic systems with the full rotation symmetry (SO(3)) and TR symmetry. In this case, we can classify all the physical operators in terms of their angular momentum. The Hamiltonian should be invariant under the full rotation and thus should carry the angular momentum 0. Both the momentum k and the spin σ carry the angular momentum 1, and thus allow to construct two terms k 2 and k · σ. As discussed above, the acoustic phonons couple to electrons through the strain tensor u ij , which is a rank-2 tensor and symmetric with respect to i and j (u ij = u ji ). For the momentum k and spin σ, one can define another rank-2 tensor T ij = k i σ j , which is actually the stress tensor of Weyl fermions. Now we need to construct the invariant terms based on two rank-2 tensors u ij and T ij and there are two ways: (1) ij u ii T jj and ij u ij T ij . Thus, we can write down the Hamiltonian as up to the order of k 2 and u ij k l , whereū = i u ii is the trace of strain tensor. For TR symmetryT , k and σ are TRodd while u ij and T ij is TR-even. Therefore, all the terms above are allowed for a TR-invariant Hamiltonian. We may further consider the strain as a field with spatial and temporal variations and thus require to change The corresponding Hamiltonian can be written as . Since we only focus on the spatial component, the upper and lower indices do not have specific meaning as in general relativity.
Here we notice that there are two independent parameters that characterize the coupling between strain tensor u ij and the stress tensor T ij . This conclusion can also be obtained by classifying the rank-2 tensors according to their angular momentum. The strain tensor can be decomposed into two parts, u J=0;M =0 = i u ii =ū with the angular momentum 0 and u J=2;M =±2,±1,0 with the angular momentum 2. The explicit form of u J=2;M is given by Similarly, the stress tensor component with the angular momentum 2 can also be written as The term m u 2,m T 2,−m is invariant under the rotation. Collecting all the invariant terms, we obtain the form of the Hamiltonian One can easily check that this Hamiltonian is the same as Eq. (B1). When reducing the symmetry group from SO(3) to chiral point group symmetry, the number of independent parameters will increase. Here we take the O group as an example and other chiral symmetry groups can be worked out in a similar manner. According to the character table of O group 86 , there are five irreducible representations (IRRs), labelled by A 1 , A 2 , E, T 1 and T 2 . The Hamiltonian should be invariant and thus belongs to the A 1 IRR. Both k and σ belong to the T 1 and since T 1 ⊗ T 1 = A 1 + E + T 1 + T 2 , we have one term k · σ can be constructed to be invariant. For the strain tensor, u 0,0 = i u ii belongs to the A 1 IRR, (u xx − u yy , 2u zz − u xx − u yy ) belongs to the E IRR, and (u xy , u yz , u zx ) belongs to the T 2 IRR. This allows us to construct the following invariant terms: (1) u 0,0 k · σ and u 0,0 k 2 ; (2) (2u zz − u xx − u y y)k z σ z + (2u xx − u yy − u zz )k x σ x + (2u yy − u xx − u zz )k y σ y ; (3) u xy (k x σ y + k y σ x ) + u yz (k y σ z + k z σ y ) + u zx (k z σ x + k x σ z ). Therefore, the first four terms in the Hamiltonian (B1) remain the same while the last term is changed to which possesses two independent parameters g 1 and g 2 .

Appendix C: Effective action and Correlation function
In this section, we will show the derivation of the formalism for the effective action for the phonon dynamics and the corresponding correlation functions.
From the Hamiltonian (A40), the effective action can be written as where we have used the imaginary time τ = it and e j a = δ j a + ∆ j a . To implement the perturbation calculations, we may separate the full action into two parts S ef f = S 0 + S 1 , where and Now let's define the Fourier transform aŝ Then, we have and Here iω n = i (2n+1)π β and iν m = i 2mπ β . We have dropped the quadratic term 2 k 2 2m * 0 in S 0 and instead, we limit the momentum summation within the range set by the cut-off Λ = is at the same order as the linear term v f Λ.
Below, we introduce the notationk = (iω n , k) for short so that the summation over both the frequency and the momentum can be simplified as k = 1 βV iωn,k . One should keep in mind the value of iω n is for boson or for fermion operators. The full action can be written as where with the stress tensor operator T a i (k,k ) = v f 2 (k i + k i )σ a for Weyl fermions. Since the full action is quadratic in fermion operators, one can directly integrate out the Weyl fermions to get the effective action for the A 0 and ∆ fields. Let's consider the partition function and the zero order partition function where Det is the determinant. The effective action W [A 0 , ∆] is defined as Z = Z 0 e −W and thus which gives us the perturbation expansion. The first order term is given by T r(XG 0 ) = k X(k,k)G 0 (k). This involve only the terms A 0 (q = 0) and ∆ i a (q = 0) in X. Since we only consider the fluctuation for X and thus can choose A 0 (q = 0) = 0 and ∆ i a (q = 0) = 0 when defining S 1 . Thus, the first order term will be zero. Next we consider the second order term, given by where three types of correlation functions are defined as Here Π 0 is the density-density correlation function, Φ ab ij is the stress-stress correlation function, while Θ a i is the stress-density correlation function.

Appendix D: A brief review of Nieh-Yan anomaly
In this section, we will first review the Weyl/Dirac fermions in the curved space and the Nieh-Yan anomaly. Then we will discuss the connection of the stress-stress correlation function Φ ab ij and Nieh-Yan anomaly in our formalism. In a curved space, the action for Dirac fermions can be written as 33,34 where γ a s are the standard γ-matrices with {γ a , γ b } = 2η ab and η ab is the metric in the Minkowski space, and e µ a is the frame field. Here the Einstein summation rule has been assumed. The frame fields satisfy η ab e µ a e ν b = g µν with the metric g µν in the curved space, so that {γ µ , γ ν } = 2g µν where γ µ = γ a e µ a . ∇ µ = ∂ µ + Γ µ is the covariant derivative with the spin connection Γ µ = i 2 η ac Γ c bµ J ab and the generator J ab = − i 4 [γ a , γ b ] of the Lorentz group. It is convenient to introduce the so-called spin connection 1-form Γ c b = Γ c bµ dx µ and the coframe field 1-form e a =ē a µ dx µ , where the coframe field is the inverse of the frame field, e ν aē a µ = δ ν µ . The Cartan's structure equations of the Riemann-Cartan spacetime 87 define the curvature 2-form R a b = dΓ a b +Γ a c ∧Γ c b and the torsion 2-form T a = de a +Γ a c ∧e c . In the standard general relativity, the torsion is assumed to be zero (zero torsion constraint), leading to the relation between the coframe field and the spin connection (de a + Γ a c ∧ e c = 0). However, in the Einstein-Cartan theory 87,88 , the torsion can be non-zero and consequently, the frame field (or coframe field) and the spin connection should be treated as two independent fields.
By comparing the effective action (A40) for our system with the action of Dirac field in the curved space, it is clear that the electron-phonon coupling can create a non-trivial coframe field e a but the spin connection Γ c b is still zero. Consequently, the Weyl fermions in our system will feel nonzero torsion, given by T a = de a , but zero curvature R a b = 0. This is in sharp contrast to the standard general relativity with non-zero curvature but zero torsion, and is known as the "Weitzenböck spacetime", which was studied in the "teleparallel gravity theory" 59 .
When the Dirac fermion is coupled to an electromagnetic field, quantum correction can give rise to the anomalous non-conservation of the chiral charge, known as the "chiral anomaly", given by where F µν = ∂ µ A ν − ∂ ν A µ is the strength of the gauge field. The chiral anomaly effect plays an essential role in predicting and understanding a number of physical phenomena. In our system, the phonon will induce the torsion field for the KW fermions, instead of the pseudo-gauge field. It turns out that the torsion field can also contribute to the non-conservation of the chiral charge, which is known as the Nieh-Yan anomaly 10-13 and given by Since the spin connection vanishes, the torsion field T a = de a gives the field strength of the coframe field e a and its components are given by T a µν = ∂ µē a ν − ∂ νē a µ and the Nieh-Yan anomaly equation is also written as Unlike the chiral anomaly with the universal dimensionless coefficient 1 16π 2 (Here we follow the convention of high energy physics and set e = = 1), the coefficient F in the Nieh-Yan anomaly has the dimension [1/L] 2 , where [L] labels the dimension of length. This can be easily obtained from the dimension counting as follows. In the Dirac action, the fermion field has the dimension [1/L] 3/2 and the chiral current j 5µ = ψ γ 5 ψ has the dimension [1/L] 3 . Thus, the left hand side of the anomaly equation (∂ µ j 5µ ) has the dimension [1/L] 4 . On the other hand, the coframe fieldē a ν is dimensionless and thus ∂ µē a ν ∂ λē b ρ has the dimension [1/L] 2 . This immediately means that the coefficient F has the dimension [1/L] 2 . It turns out that F depends on the cut-off Λ in the momentum space and is given by F = Λ 2 4π 2 12,42 (the factor 4π 2 is not important since one can always re-scale the cut-off Λ). In contrast, it was recently proposed 42,43,45 that at a finite temperature, the coefficient takes the form F = F 0 + F 1 (k B T ) 2 with F 0 = Λ 2 4π 2 and F 1 = − 1 12 . The temperature dependent term is proportional to (k B T ) 2 , which absorbs the dimension [1/L] 2 , and thus, the coefficient F 1 becomes dimensionless and universal in the sense that it is only proportional to the central charge of (1+1)-dimensional Dirac fermion 42 . Thus, the anomaly induced by this temperature dependent term is dubbed "thermal Nieh-Yan anomaly".
In order to compare with our results, it is convenient to derive the effective action that corresponds to the Nieh-Yan anomaly. Since the right hand side of Eq. (D4) is actually a total derivative, the chiral current can be given by where A 5µ is the chiral gauge potential. Since the phonon can only induce chiral chemical potential in our case, we only keep the A 50 term in the above action and thus obtain where ν, λ, ρ = x, y, z. In our effective action (1) of the main text, the frame field is given by e j a = δ j a + ∆ j a and the corresponding coframe field isē a j = δ a j − ∆ j a with a, j = x, y, z. Since we only concern the spatial coordinate here, the upper and lower indices do not have much meaning. Furthermore, ∆ j a is proportional to the strain tensor u ja which is symmetric with respect to j and a. We next focus on the stress-stress correlation function term (C17) in the effective action (C15). We may expand Φ ab ij (q) as a function ofq, Φ ab Here we only focus on the term that is linearly proportional to q and denote Φ ab ij,l = (∂ q l Φ ab ij )q =0 . Then, the corresponding term in the action is given by As demonstrated in the next Sec. E, the coefficient Φ ab ij,l is purely imaginary and proportional to δ ij abl , where is the Levi-civita symbol. Moreover, this term is linearly proportional to the chemical potential µ for a small µ. Therefore, by replacingē a ν with ∆ i a , W N Y takes the exact form as S N Y for the Nieh-Yan anomaly, and we call this term as the Nieh-Yan term below. Actually, the symmetry analysis of the possible terms for the strain tensors, given in the Sec. G, suggests that the Nieh-Yan term is the only term that is allowed by symmetry at this order of q in a uniform isotropic system.

Appendix E: Analytical evaluation of the Nieh-Yan term
In this section, we will evaluate the coefficient Φ ab ij,l analytically. With Eq. (C17), we have The Green's function G 0 is given in Eq. (C10) and here we focus on the situation with a small chemical potential µ, so we expand G 0 up to the linear order in µ as G 0 (µ) = G 0 (µ = 0) + (∂ µ G 0 ) µ=0 µ with Below we absorb v f into the definition of the momentum k, but keep track on the chirality of the Weyl fermions. Thus, we define the chirality as χ = sign( v f ), and obtain Thus, the Green's function is given by Let's first look at the first term in Eq. (E1) and the direct calculation gives In the second step of the above derivation, we collect all the terms that have even order in the momentum k since these are the nonzero term under the integral of the momentum angle. Similar calculation also shows the second term in Eq. (E1) is also zero. Thus, only the last term is left non-zero. To evaluate this term, we first need to calculate ∂ k l G 0 (k), which is given by 2k l D 4 − µ D 6 (2k l + 2iω n χσ l )((iω n ) 2 − k 2 ) + 4k l ((iω n ) 2 + k 2 + 2iω n χk · σ) By substituting these terms into the last term in Eq. (E1), we find In the above, since only the terms with even number of the momentum k can be non-zero, we thus obtain two terms, and let's calculate them separately. For the first term, we have In the last step above, we have dropped all the terms with odd number of iω n since the frequency summation will make them vanishing. For the second term, we have Again here we only pick up the terms with even order of iω n . Finally, we can add these two terms together and obtain In the above expression, the integral over the momentum k can only be non-zero when i = j. Furthermore, the integral over the polar and azimuthal angles of the momentum gives rise to with k = (k, θ, ϕ) in the spherical coordinate. Thus, let's define the function where a cut-off Λ in the momentum space has been assumed, and from which one can see that Φ ab ij,l is indeed pure imaginary (F 0 is real, as shown below), and proportional to alb δ ij . Thus, we demonstrate the effective action (H9) together with Eq. (E15) indeed takes the form of Nieh-Yan term.
Next we will evaluate the coefficient F 0 in Eq. (E14) analytically. Let's separate F 0 into F 01 and F 02 with and so that Here we have substitute the expression for D 2 . We perform the frequency summation in F 01 and consider the contour integral in the complex-z plane with the radius R of z. In the limit R → ∞, since the integrand f (z)n F (z) decays to zero fast enough, this integral should be zero. On the other hand, this integral can re-expressed as the integral around the poles of the integrand in the complex-z plane, and thus is determined by the residuals of these poles according to the residual theorem. We next need to figure out all the poles and the corresponding residuals for the integrand f (z)n F (z). For the Fermi function n F (z) = 1 e βz +1 , the poles are z n = iω n = i (2n+1)π β and the corresponding residuals for f (z)n F (z) are Res[f (z)n F (z), z = z n ] = − 1 β f (iω n ). There are another two poles z ± = ±k from the function f (z). Since they are not simple poles, we need a bit more work to extract the residuals. For z = z + = k, let's denote f (z)n F (z) = g(z) (z−k) 2 with g(z) = n F (z) (z+k) 2 . We need to expand g(z) around z + = k and pick up the first order term ∂g(z) ∂z z=k (z − k) since the coefficient in this term will give the residual of f (z)n F (z). Direct calculations give rise to and thus the corresponding residual is we obtain Putting all the results together, we have and We may consider the coefficient F 01 at the zero temperature and its correction from the finite temperature, F 01 = F 01 (T = 0) + δF 01 (T ). At zero temperature (β → ∞), the first term in the above integral diverges as Λ 2 since tanh(βΛ/2) → 1 while the second term converges (actually vanishes). Thus, we have and the finite temperature correction is given by where we have used x = βk/2 and Λ → ∞.

Appendix F: Numerical Methods for Evaluating Stress-stress Correlation Function
In this section, we will describe our numerical methods to evaluate the stress-stress correlation functions. Here we only focus on the stress-stress correlation function Φ ab ij (q, iν m ) in Eq. (C17). One can re-write the Green's function in Eq. (C10) as where ξ s,k = s v f k − µ and the projection operator P s (k) = 1 2 1 + s k·σ k . Then after the Matsubara frequency summation, Eq. (C17) can be simplified as Direct evaluation of G ab ij gives Up to now, the whole formalism is general and can be calculated numerically in principles. As our purpose here is to extract the term that is related to the Nieh-Yan anomaly, we want to look for the terms Φ ab ij (q, iν m ) ∝ δ ij . Furthermore, let's choose the momentum along the z-direction, q = (0, 0, q z ), and then we should expect that the relevant component should be Φ xy ii (q z , iν m ). Thus, our numerical simulations focus on the component Φ xy xx (q z , iν m ) below (Φ xy xx (q z , iν m ) = Φ xy yy (q z , iν m )). For this component, we have Now let's consider the momentum integral in Eq. (F2), and choose the spherical coordinate for the momentum, k = (k sin θ k cos ϕ k , k sin θ k sin ϕ k , k cos θ k ).
Furthermore, we denote k = k − q = (k sin θ k cos ϕ k , k sin θ k sin ϕ k , k cos θ k ) = (k sin θ k cos ϕ k , k sin θ k sin ϕ k , k cos θ k − q z ). Thus, we have ϕ k = ϕ k , k sin θ k = k sin θ k and k cos θ k = k cos θ k − q z . The latter two equalities give rise to k = k 2 + q 2 z − 2kq z cos θ k and cos θ k = k cos θ k −qz k . By choosing q along the z direction, ξ s,k−q only involves k and θ k and does not depend on ϕ k . Thus, we can definē and then We may decompose Φ xy xx (q z , iν m ) into several different components Φ xy xx (q z , iν m ) = Φ xy,− xx (q z , iν m ) + Φ xy,+ xx (q z , iν m ), where ± labels the contribution from different s bands. We further decompose Φ xy,± xx (q z , iν m ) = Φ xy,± xx,1 (q z , iν m ) + Φ xy,± xx,2 (q z , iν m ), where 1 and 2 here are for the intra-band and inter-band contribution. Explicitly, . (F8) Since the above integral involves the singular points, we need to treat it as the Cauchy principal value integral. To numerically evaluate the integral, we need some more simplification. First, the integral over θ k can be replaced by the integral over k . All the above integrals takes the form of With d 3 k δ(k − k + q) = d 3 k δ(k − k 2 + q 2 z − 2kq z cos θ k ) = 1, we have 89 Here the integral range of k is determined by requiring | cos θ k | ≤ 1. Next one may absorb the velocity v f into the definition of the momentum,k = v f k,k = v f k , q = v f q z . Thus, one can see that G does not give any additional factors while the integral J 1 should contain an additional pre-factor 1 ( v f ) 3 . Now G should be a function ofk,k , q, denoted as G(k,k , q). Finally, one can take the transformationk = 1 2 (x + y) andk = 1 2 (x − y) and the integral can be transformed into Using the equation, we can rewrite the integrals as whereḠ xy xx (ss ; x, y, q) = and we have replaced iν m by ω + iη. Now one can see that the poles of the integral are determined by x = ±ω or y = ±ω, and the corresponding Cauchy integral is easy to deal with numerically. These expressions form the basis for our numerical calculations.
To recover the unit, one notices that q, x, y in the above expression are in the unit of energy. v f is in the unit of [E · L], and thus we need to further choose certain length unit of our system, which is denoted as a 0 . Then, the unit of energy can be chosen as v f /a 0 . From the expression of Φ xy xx , one can easily find that Φ xy xx is in the unit of energy density or equivalently v f /a 4 0 . For the numerical calculations, we simply choose v f = 1 and a 0 = 1, and this is also convenient to compare with analytical results in Sec. E. The main results of our numerical calculations are summarized in the Fig. 2 in the main text.
Below we provide some more understanding on the analytical aspects of the formalism. We consider the small q expansion of Eq. (F7) at iν m = 0, and perform the perturbation expansion. The calculation here is equivalent to the calculations in Sec. E, except that we below keep the chemical potential up to µ 3 . For a small q = q zêz , we have the perturbation expansion up to the first order in q z . Now let's denote Φ xy xx = s,s Φ xy xx,ss and evaluate Φ xy xx,ss separately. First, let's consider s = s and we haveḠ which only possesses the linear-q z term. Thus, With ∂n F ∂ξ s,k = − βe βξ s,k (1+e βξ s,k ) 2 and x = βξ s,k , we find we have which contains both q z -independent and linear-q z terms. In addition, up to the linear order in q z . Putting all these results together, we have Φ xy xx,s,−s = where we have used In the above expression, since s (n F (ξ −s,k ) − n F (ξ s,k )) = 0, the first term will vanish after the summation over s. Thus, Now let's define x = βξ −s,k and find Putting all the results together, we have The second term reproduces the result of Eqs. (E15) and (E39) for the thermal Nieh-Yan anomaly (choosing v f = 1), while the first term is of order µ 3 and thus not included in the calculations of Sec. E, which only keep the terms up to linear µ. Therefore, this analytical result is consistent with the results given in Sec. E. The expression (F33) can fit well with our numerical results at small µ, as shown in Fig. 2d in the main text.
where Λ 0 (q) = C 2 1 Π 0 (q) + i,j We can expand Λ 0 (q), Λ a 2,i (q) and Λ ab 3,ij (q) as a function ofq, Λ(q) = Λ(0) + (∂ q Λ)q =0 · q + (∂ ωn Λ)q =0 ω n + .... At the zero order term, we have Here all the terms are quadratic in the strain tensor u ij and let's consider the general form of the effective action for acoustic phonons (elastic wave), which is defined as 90,91 where the rank-4 tensor λ ijkl is the elastic modulus. It is clear that W (0) just provides the correction of the elastic modulus λ ijkl . Next we focus on the first order term with the linear q dependence, which takes the form is just a number. From Eq. (G4), one can see that the last term is just the Nieh-Yan term. Although the above expression seems quite complex, the symmetry gives a strong constraint on the form of the allowed effective action. It is clear that the above action contributes to the following general form of the effective action of acoustic phonons where ∂ i = ∂ ∂ri . Now the question is what is the most general form of the coefficient ξ ijklm , a rank-5 tensor, for a uniform isotropic system. Since the strain tensor u ij is symmetric with respect to the indices i and j, the coefficient ξ ijklm should satisfy the relations ξ ijklm = ξ ikjlm = ξ ijkml . Furthermore, it also satisfies the anti-symmetric relation ξ ijklm = −ξ ilmjk . These relations together with the full rotation symmetry will fix the form of ξ ijklm .
It is also equivalent to ask how to construct an invariant term with the form of (∂ i u jk )u lm . We analyze this problem based on the irreducible representation of the angular momentum. As discussed in Sec. B, the symmetric rank-2 tensor u ij can be decomposed to two parts,ū with the angular momentum 0 and u J=2,M with the angular momentum 2. Since the expression involves two strain tensors, the corresponding angular momentum decomposition can be given by where the bold number represents the total angular momentum. Furthermore, the expression (∂ i u jk )u lm involves one ∂ i in addition to two strain tensors, and ∂ i possesses the angular momentum 1. To make the term (∂ i u jk )u lm invariant (with the angular momentum 0), only the irreducible representation with the angular momentum 1 in 2 ⊗ 2 can lead to 1 ⊗ 1 = 2 ⊕ 1 ⊕ 0. From such angular momentum combination, only a single term is allowed for ξ ijklm in a uniform isotropic system. This term can be explicitly constructed as so that the parameter ξ ijklm given by We can apply the above analysis to the action W (1) , and it is clear that it is impossible to construct an invariant term from the first two terms in Eq. (G7), and thus the coefficient ∂ j Λ 0 and ∂ j Λ a 2,i must vanish. The effective action then takes the form with the coefficient Here we have used Eq. E15. Substituting the coefficient and also changing from imaginary time to real time, we have Comparing with S ph,1 in Eq. (G12), one can see that the Nieh-Yan term can contribute to a new term in the phonon effective action and the corresponding coefficient is given by In the above discussion, we have dropped the term with the form of u ij ω n u kl and in the time-domain, this term should have the form of u ij (∂ τ u kl ). In the elastic theory, this term is nothing but the viscosity, generally defined as [90][91][92] However, the existence of the viscosity requires the breaking of TR symmetry, while TR symmetry exists in our system. Thus, we will expect the viscosity coefficient η ijkl extracted from our model should vanish. Now we hope to solve the acoustic phonon dispersion (or elastic wave) in our system based on the actions (G6) and (G12). For the uniform isotropic system, there are two independent parameters for the elastic moduli λ ijkl and the corresponding equation of motion for the displacement field u is given by where c l and c t are the velocities of longitudinal and transverse waves. The Nieh-Yan term S ph,1 can add an additional term into the above elastic wave equation, as given by By including the Nieh-Yan term, the equation of motion for the displacement field is now written as We can decompose the displacement field into the longitudinal and transverse components, u = u l + u t . One can show that the longitudinal and transverse components are still decoupled from each other. The longitudinal part does not have any correction from the Nieh-Yan term, while the equation of motion for the transverse part is revised as Let's take u t = u 0 e iq·r−iωt and the equation of motion is changed to Now let's choose q = qê z and u 0 = (u 0x , u 0y , 0) T (the superscript T here represents transpose) and then the equation of motion (G23) can be written as This eigen equation can be easily solved. Two branches of dispersion relations for the transverse modes are given by ω t ± = c 2 t q 2 ± |ξ0| 2 q 3 = q c 2 t ± |ξ0| 2 q, from which one can see that the Nieh-Yan term mainly corrects the phonon velocity and gives rise to the different velocities for two transverse phonon modes. The corresponding normalized eigen vectors are given by where s = ±, i = x, y, z, and (M i ) jk = (−i) ijk . One can easily show that only l s,z is non-zero and l ±,z = ± for the momentum q = qê z . Due to the full rotation symmetry of our model, the general form of angular momentum for the s-mode is given by l s = s q withq = q q . Next let's consider the phonon total angular momentum induced by a temperature gradient. The total phonon angular momentum per volume can be related to the temperature gradient by In the above expression, the integral of x is convergent, and thus just gives a number. Therefore, all the temperature dependence is given by the coefficient before the x-integral. Since ξ 0 = g 2 1 µ 5 F 0 with F 0 = F 0 + F 1 (k B T ) 2 , we expect α zz has T 3 dependence for the standard Nieh-Yan term and T 5 dependence for the thermal Nieh-Yan term. Let's denote a 0 = x 4 dx e x (e x −1) 2 + 1 2 e x (e x −1) 3 (e x (1 − x) − x − 1) = −4π 4 /15, and then we have

(G37)
Appendix H: Electron-electron Interaction Since the electron-phonon interaction may be influenced by electron-electron interaction, one may wonder how electron-electron interaction affects the proposed phenomena here. In this section, we will demonstrate that electronelectron interaction is actually decoupled from the phonon modes discussed in this work and thus all the results remain valid even taking into account electron-electron interaction.
Next let's consider the path integral and we can perform the transformatioñ ϕ(q) = ϕ(q) + ieΘ a i (−q)∆ i a (q)D(q) (H11) ϕ(−q) = ϕ(−q) + ieΘ a i (q)∆ i a (−q)D(−q). Correspondingly, and W 1 [ϕ, ∆] = 1 2 q D −1 (q)φ(q)φ(−q) + Φ ab ij (q) + e 2 Θ a i (−q)D(q)Θ b j (q) ∆ i a (q)∆ j b (−q) . From this expression, we conclude that the full stress-stress correlation function, denoted asΦ ab ij (q), should generally acquire a correction from the electron-electron interaction through Φ ab ij (q) = Φ ab ij (q) + e 2 Θ a i (−q)D(q)Θ b j (q), where the stress-density correlation function Θ a i (q) is defined in Eq. (C18). One can see that the Θ a i (q) function couples electrons to acoustic phonons and this function involves one density vertex and one stress tensor vertex from Eq. (C18). Next we will show that for the phonon modes that we are interested in, the Θ a i (q) function is always zero if the system has full rotation symmetry. In the above expressions, a, b, i, j = x, y, z which is not convenient for the symmetry analysis. It is more convenient to relabel the strain and stress tensor with the angular momentum, as shown in Eq. (B3)-(B12). The rotationally symmetric form of the Hamiltonian (B13) is useful for our symmetry analysis below.
To make the problem more concrete, let's consider the phonon momentum along the z-direction, q = q zêz , and we first figure out which strain tensors, as well as the stress tensors, are involved in the Nieh-Yan term. To see that, let's write down explicitly the terms in Eq. (G12) along the q zêz direction, which involves δ jl zkm (∂ z u jk )u lm = zxy (∂ z u jx )u jy + zyx (∂ z u jy )u jx = ((∂ z u xx )u xy + (∂ z u yx )u yy + (∂ z u zx )u zy ) − ((∂ z u xy )u xx + (∂ z u yy )u yx + (∂ z u zy )u zx ) = 2c∂ z (u xx − u yy ) + 2(∂ z u zx )u zy ,