Many-body fermionic excitations in Weyl semimetals due to elastic gauge fields

TL;DR

This paper investigates how elastic strain-induced gauge fields affect the single-particle spectrum of Weyl semimetals, revealing new anisotropic fermionic excitations caused by electron-phonon interactions.

Contribution

It demonstrates the emergence of well-defined fermionic satellite peaks in Weyl semimetals due to electron-phonon coupling, a phenomenon not seen in Dirac systems with Coulomb interactions.

Findings

Satellite peaks appear at lowest order in perturbation theory.

The new excitations are anisotropic.

Behavior is regulated by the Debye frequency.

Abstract

We study the single-particle spectrum of three-dimensional Weyl semimetals taking into account electron-phonon interactions that are the result of straining the material. We find that a well-defined fermionic excitation appears in addition to the standard peak corresponding to quasiparticle states as suggested by Landau-Fermi liquid theory. Contrary to the case of Dirac systems interacting via the Coulomb interaction, these satellite peaks appear even at lowest order in perturbation theory. The new excitations are anisotropic, as opposed to the single-particle spectrum, and their behavior is dictated by the Debye frequency, which naturally regulates the electron-phonon coupling.

Figures

Equations

$${\mathcal{H}}_{0,\chi }^W(\mathbf{k})=v{\mathbf{\sigma }}_{\perp }\cdot {\mathbf{k}}_{\perp }-\chi v_3(k_3-\chi {\lambda }_3){\sigma }_3,$$

$$\delta \mathcal{H}[u]=v{g}_0{\sigma }_0\text{Tr}[u]+v{g}_1{\sigma }_i{u}_{ij}{\lambda }_j+v{g}_2{\sigma }_i{\delta }_{ij}{\lambda }_j\text{Tr}[u],$$

$$\delta {\mathcal{H}}_2[u]=\frac{i{v}_F{g}_1}{2}{\int }_{\mathbf{k},\mathbf{q}}{\psi }_{\mathbf{k}+\mathbf{q}}^{†}{\sigma }_i{\lambda }_j(q_i{\delta }_{jr}+q_j{\delta }_{ir})u_r(\mathbf{q}){\psi }_{\mathbf{k}},$$

$$g_{\mathbf{q}}^{\alpha }=\frac{i}{2}\frac{v_F{g}_1}{\sqrt{2\rho {\omega }_{\alpha }(\mathbf{q})}}{\sigma }_i{\lambda }_a(q_i{\delta }_{ar}+q_a{\delta }_{ir}){\hat{e}}_r^{\alpha },$$

$$D_0^{\alpha {\alpha }^{\prime }}(i{\omega }_m,\mathbf{q})=\frac{2{\omega }_{\alpha }(\mathbf{q})}{(i{\omega }_m{)}^2-{\omega }_{\alpha }^2(\mathbf{q})}\vartheta (k_D^{\alpha }-|\mathbf{q}|){\delta }_{\alpha {\alpha }^{\prime }},$$

$$\Sigma (\omega ,\mathbf{k})$$

$$\displaystyle\phantom{=}\qquad\qquad+\frac{u\Lambda_{+}(\bm{k}-\bm{q})\vartheta(\mu-\varepsilon_{\bm{k}-\bm{q}})}{\omega^{+}-\varepsilon_{\bm{k}-\bm{q}}+u\omega_{\alpha}({\bm{q}})}\bigg{]}g^{\alpha}_{\bm{q}},$$

$$\frac{\Sigma(\omega,{\bm{k}})}{v_{F}|{\bm{\lambda}}|}=G_{1}\big{[}\Sigma_{0}(\omega,{\bm{k}})\sigma_{0}+\Sigma_{1}(\omega,{\bm{k}})\bm{\sigma}\cdot\hat{{\bm{k}}}\big{]},$$

$${\Sigma }_1^L(\omega ,\mathbf{k})={\hat{\lambda }}^i{\hat{\lambda }}^j{\int }_{k_D>|\mathbf{q}|}\frac{{\hat{q}}^i{\hat{q}}^j}{|\mathbf{k}-\mathbf{q}|}{H}_{-}(\omega ,\mathbf{k}-\mathbf{q},\mathbf{q})$$

$$\displaystyle\phantom{=}\qquad\qquad\qquad\times\big{[}(2{\bm{k}}\cdot{\bm{q}}-|{\bm{q}}|^{2})(\hat{{\bm{k}}}\cdot\hat{{\bm{q}}})-|{\bm{k}}||{\bm{q}}|\big{]},$$

$${\Sigma }_n^{\theta }(\omega ,\mathbf{k})=\frac{1}{8}{\int }_{k_D>|\mathbf{k}+\mathbf{q}|}\frac{f_n^{\theta }(\mathbf{k},\mathbf{q})H_{{\eta }_n}(\omega ,\mathbf{q},\mathbf{k}+\mathbf{q})}{|\mathbf{k}+\mathbf{q}{|}^{2n+1}},$$

$$H_{\pm }(\omega ,\mathbf{x},\mathbf{y})\equiv \frac{\vartheta (\mu -|\mathbf{x}|)}{{\omega }^{+}-|\mathbf{x}|+\epsilon |\mathbf{y}|}-\frac{\vartheta (\mu -|\mathbf{x}|)}{{\omega }^{+}-|\mathbf{x}|-\epsilon |\mathbf{y}|}$$

$$+\frac{1}{{\omega }^{+}-|\mathbf{x}|-\epsilon |\mathbf{y}|}\pm \frac{1}{{\omega }^{+}+|\mathbf{x}|+\epsilon |\mathbf{y}|}.$$

$$\!\!\rho^{\theta}_{\pm}(\omega,{\bm{k}})=\frac{-G_{1}\text{Im}[\Sigma^{\theta}_{\pm}]/\pi}{\big{(}\omega\mp k-G_{1}\text{Re}[\Sigma^{\theta}_{\pm}]\big{)}^{2}+\big{(}G_{1}\text{Im}[\Sigma^{\theta}_{\pm}]\big{)}^{2}},$$

$$\Sigma (i{\omega }_m,\mathbf{k})=-\frac{1}{\beta }\sum _{i{\omega }_n,\alpha ,{\alpha }^{\prime }}{\int }_{|\mathbf{q}|<k_D^{\alpha }}{g}_{-\mathbf{q}}^{\alpha }{G}_0(i{\omega }_m-i{\omega }_n,\mathbf{k}-\mathbf{q})g_{\mathbf{q}}^{{\alpha }^{\prime }}{D}_0^{\alpha {\alpha }^{\prime }}(i{\omega }_n,\mathbf{q}).$$

$$G_0(i{\omega }_n,\mathbf{k})=\frac{1}{2}\sum _{a=\pm }\frac{{\sigma }_0+a\mathbf{\sigma }\cdot \hat{\mathbf{k}}}{i{\omega }_n+\mu -a{\epsilon }_{\mathbf{k}}}.$$

$$\Sigma (i{\omega }_m,\mathbf{k})$$

$$I_{ab}^{\alpha }(\mathbf{k},\mathbf{q},i{\omega }_m)$$

$$=\frac{1}{2\pi i}{\oint }_C\text{d}z\frac{N_{\text{B}}(z)}{(i{\omega }_m+\mu -z-a{\epsilon }_{\mathbf{k}-\mathbf{q}})(z-b{\omega }_{\mathbf{q}}^{\alpha })},$$

$$\frac{N_{\text{B}}(z)}{(i{\omega }_m+\mu -z-a{\epsilon }_{\mathbf{k}-\mathbf{q}})(z-b{\omega }_{\mathbf{q}}^{\alpha })}=\frac{N_{\text{B}}(z)}{i{\omega }_m+\mu -(a{\epsilon }_{\mathbf{k}-\mathbf{q}}+b{\omega }_{\mathbf{q}}^{\alpha })}\left(\frac{1}{z-b{\omega }_{\mathbf{q}}^{\alpha }}-\frac{1}{z-i{\omega }_m-\mu +a{\epsilon }_{\mathbf{k}-\mathbf{q}}}\right).$$

$$I_{ab}^{\alpha }(\mathbf{k},\mathbf{q},i{\omega }_m)=\frac{N_{\text{B}}(b{\omega }_{\mathbf{q}}^{\alpha })-N_{\text{B}}(i{\omega }_m+\mu -a{\epsilon }_{\mathbf{k}-\mathbf{q}})}{i{\omega }_m+\mu -a{\epsilon }_{\mathbf{k}-\mathbf{q}}-b{\omega }_{\mathbf{q}}^{\alpha }}.$$

$$N_{\text{B}}(i{\omega }_m+\mu -a{\epsilon }_{\mathbf{k}-\mathbf{q}})=-N_{\text{F}}(\mu -a{\epsilon }_{\mathbf{k}-\mathbf{q}}),$$

$$N_{\text{B}}(-{\omega }_{\mathbf{q}}^{\alpha })$$

$$N_{\text{F}}(\mu -{\epsilon }_{\mathbf{k}})$$

$$I_{++}^{\alpha }(\mathbf{k},\mathbf{q},i{\omega }_m)$$

$$I_{+-}^{\alpha }(\mathbf{k},\mathbf{q},i{\omega }_m)$$

$$I_{-+}^{\alpha }(\mathbf{k},\mathbf{q},i{\omega }_m)$$

$$I_{--}^{\alpha }(\mathbf{k},\mathbf{q},i{\omega }_m)$$

$$\Sigma (\omega ,\mathbf{k})$$

$$g_{-\mathbf{q}}^L{\Lambda }_{\pm }(\mathbf{k}-\mathbf{q})g_{\mathbf{q}}^L$$

$$\displaystyle=\frac{\tilde{g}_{1}^{2}({\bm{\lambda}}\cdot{\bm{q}})^{2}}{\omega_{L}({\bm{q}})}\bigg{[}\sigma^{0}\pm\frac{\big{[}2{\bm{k}}\cdot\bm{q}-\bm{q}^{2}\big{]}(\bm{q}\cdot\bm{\sigma})}{|{\bm{k}}-{\bm{q}}||\bm{q}|^{2}}\mp\frac{({\bm{k}}\cdot{\bm{\sigma}})}{|\bm{k}-\bm{q}|}\bigg{]},$$

$$\sum _{\alpha =T_1,T_2}{g}_{-\mathbf{q}}^{\alpha }{\Lambda }_{\pm }(\mathbf{k}-\mathbf{q})g_{\mathbf{q}}^{\alpha }$$

$$\displaystyle=\frac{\tilde{g}_{1}^{2}}{4\omega_{T}(\bm{q})}\big{[}{\bm{\lambda}}^{2}\bm{q}^{2}+({\bm{\lambda}}\cdot\bm{q})^{2}\big{]}\sigma^{0}\pm\frac{\tilde{g}_{1}^{2}}{4\omega_{T}(\bm{q})|\bm{k}-\bm{q}|}\big{[}({\bm{\lambda}}\cdot\bm{q})^{2}-{\bm{\lambda}}^{2}\bm{q}^{2}\big{]}(\bm{k}\cdot\bm{\sigma})$$

$$\displaystyle\phantom{=}\pm\frac{\tilde{g}_{1}^{2}}{4\omega_{T}(\bm{q})|\bm{k}-\bm{q}|}\bigg{\{}2\big{[}{\bm{\lambda}}^{2}-4({\bm{\lambda}}\cdot\hat{\bm{q}})^{2}\big{]}(\bm{k}\cdot\bm{q})+2({\bm{\lambda}}\cdot\bm{q})({\bm{\lambda}}\cdot\bm{k})+5({\bm{\lambda}}\cdot\bm{q})^{2}-{\bm{\lambda}}^{2}\bm{q}^{2}\bigg{\}}(\bm{q}\cdot{\bm{\sigma}})$$

$$\displaystyle\phantom{=}\pm\frac{\tilde{g}_{1}^{2}({\bm{\lambda}}\cdot\bm{q})\big{[}\bm{k}\cdot\bm{q}-\bm{q}^{2}\big{]}}{2\omega_{T}(\bm{q})|\bm{k}-\bm{q}|}({\bm{\lambda}}\cdot{\bm{\sigma}}).$$

$${\Sigma }_0^L(\omega ,\mathbf{k})$$

Full text

Introduction.— Quite generically, equilibrium and transport properties of many-body metallic systems are accounted for by the Landau theory of the Fermi liquid. This theory is based on the notion of quasiparticles: single-particle fermionic excitations that are in one-to-one correspondence to the non-interacting fermionic states. In this regard, the single-particle spectral function $\rho(\omega,\bm{k})$ provides information about how good the quasiparticle approximation is in an interacting system. For non-interacting fermions the spectral function is given by a Dirac delta function, $\rho^{0}(\omega,\bm{k})=\delta(\hbar\omega+\mu-\varepsilon_{\bm{k}})$, where $\varepsilon_{\bm{k}}$ is the dispersion relation of the particle and $\mu$ the chemical potential. Upon including interactions, Landau’s Fermi liquid theory states that, because the damping rate of the quasiparticles scales quadratically with the energy at the Fermi level, the spectral function $\rho^{0}(\omega,\bm{k})$ becomes a Dirac delta function that is weighted by a parameter $0<Z<1$.

It was realized a long time ago Lundqvist (1967); Hedin et al. (1967) that the actual spectral function of an interacting fermionic system is richer than the one suggested by Landau-Fermi liquid theory. Apart from the existence of a well-defined quasiparticle peak, extra weakly-damped, i.e. long-lived, fermionic satellite excitations appear due to the coupling of bare fermions with propagating bosonic degrees of freedom. Such satellite excitations, currently called plasminos, were theoretically reported in the context of the quark-gluon plasma Klimov (1981); Weldon (1989), and more generically in hot or dense relativistic plasmas Baym et al. (1992); Blaizot and Ollitrault (1993); Manuel (2000). In condensed matter physics, they are referred to as plasmarons, and they have been theoretically studied in graphene Polini et al. (2008); Sensarma et al. (2011); LeBlanc et al. (2011, 2012) and Bismuth Chudzinski and Giamarchi (2011). Experimental evidence of these emergent modes has been found using angle-resolved photoemission spectroscopy (ARPES) and optical measurements in graphene Bostwick et al. (2010) and Bismuth Tediosi et al. (2007). These satellite fermionic excitations are not just of purely theoretical interest. They modify the transport and optical properties of the electronic system Chudzinski and Giamarchi (2011); Sadooghi and Taghinavaz (2014): in the case of plasminos, it is well known that they have an opposite helicity-to-chirality ratio as compared to the quasiparticle excitations Weldon (1982).

In the context of condensed-matter physics, Weyl semimetals are materials that host low-energy quasiparticles that linearly disperse Armitage et al. (2018). These quasiparticles also have a chirality and have similar kinematics to chiral ultra-relativistic particles. It is therefore interesting to analyze the many-body fermionic spectrum in these systems and see how many-body quantum effects interfere with ultra-relativistic (or pseudo-relativistic) invariance in solid-state systems.

Despite their appealing similarities to ultra-relativistic systems, we have to keep in mind that Weyl semimetals are condensed matter systems in which the transverse electromagnetic degrees of freedom are not important compared to the collective plasmon mode Kotov et al. (2012). Moreover, in three-dimensional systems, longitudinal and transverse plasmonic modes become gapped, with the gap set by the plasma frequency. Effects associated to electron-plasmon interactions are therefore weaker than in the two-dimensional counterpart Throckmorton et al. (2015); Hofmann et al. (2015); Kozii and Fu (2018). Phonons, on the contrary, remain gapless even after including many-body effects in three dimensions. Besides, electron-phonon interactions appear in Weyl semimetals in the form of elastic gauge fields Cortijo et al. (2015); Ilan et al. (2019), so they become a promising alternative to study the physics of emergent fermionic degrees of freedom different from the quasiparticles, even at lowest order in perturbation theory Klimov (1981); Baym et al. (1992); Blaizot and Ollitrault (1993); Vanderheyden and Ollitrault (1997); Le Bellac and Manuel (1997); Manuel (2000).

In the present work we show that novel fermionic excitations appear in Weyl semimetals due to electron-phonon couplings that can be expressed in terms of elastic gauge (vector) fields. A distinctive feature of the electron-phonon coupling is that extra energy/momentum scales appear in the problem. On the one hand there are the longitudinal/transverse phonon velocities, which are typically much smaller than the Fermi velocity of the electrons, and on the other hand there is the Debye frequency, which naturally acts as a cutoff for the electron-phonon interaction. Another interesting feature is that, as we will see, despite the isotropic nature of the original non-interacting electronic spectrum, the emergent fermionic modes will be anisotropic due the particular nature of the elastic vector fields in Weyl metals. Both the anisotropy and the presence of additional energy and momentum scales make the problem quite different from the situation in relativistic plasmas.

Strain in Weyl metals.— We start by considering the minimal model for a time-reversal symmetry breaking Weyl metal. Setting $\hbar\equiv 1$, the Hamiltonian that describes two Weyl cones separated by $2\lambda_{3}\hat{\bm{z}}$ reads

[TABLE]

with ${\bm{k}}=({\bm{k}}_{\perp},k_{3})$, ${\bm{\sigma}}=({\bm{\sigma}}_{\perp},\sigma_{3})$ Pauli matrices, $\chi=\pm$ the chirality of the cone and $v$ and $v_{3}$ Fermi velocities.

It was shown in Ref. Cortijo et al. (2015) that when strain is considered, the lowest-order modification of the Hamiltonian in Eq. (1) due to strain is $\delta\mathcal{H}[u]=v\chi\bm{\sigma}_{\perp}\cdot\bm{A}^{\text{el}}_{\perp}+v_{3}A^{\text{el}}_{3}\sigma_{3}$. From this expression, we see that the presence of the chirality $\chi$ implies that the coupling between electrons and the strain tensor is through a chiral vector potential $\bm{A}^{\text{el}}=({\bm{A}}^{\text{el}}_{\perp},A_{3}^{\text{el}})$. That is, it couples with a different sign to the two Weyl nodes with opposite chirality. This so-called elastic gauge field can be expressed in terms of the microscopic tight-binding parameters leading to the low-energy Hamiltonian in Eq. (1) and the linearized strain tensor $u_{ij}\equiv(\partial_{i}u_{j}+\partial_{j}u_{i})/2$ that itself is defined in terms of the displacement vector $u_{i}$, which measures the deviation from the equilibrium situation Cortijo et al. (2015); Landau et al. (1995).

The influence of strain can also be derived from a more generic symmetry-based approach. In this line of approach, all possible scalars that can be constructed from the available scalars, vectors and tensors are added to the Hamiltonian for the unstrained material. At first order in the strain tensor, the strain-induced Hamiltonian $\delta\mathcal{H}[u]$ for one cone contains three possible terms Arjona and Vozmediano (2018), i.e.,

[TABLE]

with $g_{0,1,2}$ coupling constants. The first term, proportional to the identity matrix $\sigma_{0}$, is the standard electron-phonon coupling coming from the changes induced in the electronic density due to changes in the volume, $\delta V/V=\text{Tr}[u]$. The Hamiltonian given by Eq. (2) gives rise to different types of electron-phonon coupling. These couplings can be derived by writing the strain tensor $u_{ij}$ in terms of the displacement $u_{i}$ and subsequently quantizing the latter to yield phonons. Phonons come in two types: longitudinal and transverse. By construction, only longitudinal phonons contribute to the trace in the first and third term of Eq. (2), whereas the second term receives contributions from both longitudinal and transverse phonons.

The coupling between transverse phononic degrees of freedom and electrons is similar to what happens in QED (quantum electrodynamics), where electrons couple to transverse gauge degrees of freedom. However, there is a crucial difference: in ultrarelativistic QED both the gauge and fermionic degrees of freedom have the same velocity, and the transverse modes are gapped as the longitudinal collective plasmonic mode Le Bellac (2000). In the case of phonons, in contrast, the phonon velocity is much smaller than the Fermi velocity, thereby resembling the case of sound modes in plasmas due to ionic collective modes. These modes do not develop a gap but they do attain a strong renormalization of the sound velocity. These observations place the vector electron-phonon coupling in Eq.(3) in an unique situation not addressed before in condensed matter, nor in relativistic plasmas.

We therefore proceed by focussing on the second term in Eq. (2). The corresponding second-quantized continuum Hamiltonian $\delta\mathcal{H}_{2}[u]$ reads in momentum space

[TABLE]

where $\int_{{\bm{k}}}\equiv\int\text{d}^{3}{\bm{k}}/(2\pi)^{3}$, $\psi_{{\bm{k}}}^{\dagger}$ ($\psi_{{\bm{k}}}$) is the fermionic creation (annihilation) operator and we denoted the Fermi velocity by $v_{F}$. We wrote Eq. (3) explicitly in terms of the displacement field $u_{r}({\bm{q}})$ using the definition of the strain tensor. Quantizing the displacement field, we find the electron-phonon vertex $g^{\alpha}_{\bm{q}}$, i.e.,

[TABLE]

with $\rho$ the mass density of the material under consideration, $\omega_{\alpha}({\bm{q}})=c_{\alpha}|{\bm{q}}|$ the dispersion relation of the longitudinal ($\alpha=L$) and transverse acoustic phonons ($\alpha=T$) and $c_{\alpha}$ the corresponding speeds of sound. Furthermore, the polarization vectors $\hat{e}^{\alpha}_{r}$ form a basis in momentum space and obey the identity $\sum_{\alpha}\hat{e}^{\alpha}_{i}\hat{e}^{\alpha}_{j}=\delta_{ij}$, with $\hat{e}^{L}_{i}=\hat{q}^{i}$.

It is important to stress that the electron-phonon coupling in Eq. (4) follows from very generic elasticity theory. It therefore holds in any Weyl semimetal. Contributions from optical phonons, in contrast, are rather material and model-dependent. Additionally, due to their defining energy-gap at the $\Gamma$-point, we expect the contributions from optical phonons to be negligible. We therefore do not consider the effect of optical phonons here and proceed with the computation of the self-energy due to the coupling of the electrons to the acoustic phonons.

Electronic self-energy.— For simplicity we consider an isotropic Weyl cone with dispersion relation $\varepsilon_{{\bm{k}}}=v_{F}|{\bm{k}}|$. Furthermore, we treat the acoustic phonons within the Debye model such that upon crossing the Debye frequency $\omega_{D}$, the phonon density of states becomes zero. The Matsubara Green’s function for the acoustic phonons therefore reads Abrikosov et al. (2012)

[TABLE]

with $i\omega_{m}$ a bosonic Matsubara frequency and $k^{\alpha}_{D}=\omega_{D}/c_{\alpha}$ the Debye wavenumber. Upon integrating out the phonons, we obtain the one-loop electron self-energy Sadovskii (2006). At $T=0$ the retarded self-energy reads Blaizot and Ollitrault (1993)

[TABLE]

with $\Lambda_{\pm}({\bm{k}})\equiv(\sigma^{0}\pm{\bm{\sigma}}\cdot\hat{{\bm{k}}})/2$ and $\omega^{+}=\omega+i0$. Note that we absorbed the chemical potential in the frequency $\omega$ and that the integral is rendered finite by the high-momentum cut-off provided for by the Debye wavenumber. The sum over the polarizations in Eq. (6) can be performed separately for the longitudinal and transverse phonons. The full self-energy is then given by the sum of the longitudinal and transverse contribution, i.e., $\Sigma(\omega,{\bm{k}})\equiv\Sigma^{L}(\omega,{\bm{k}})+\Sigma^{T}(\omega,{\bm{k}})$. For simplicity we will take $c_{L}=c_{T}\equiv c$ in the following.

Before we compute the self-energy, it is instructive to introduce dimensionless variables by scaling all parameters on $|{\bm{\lambda}}|$, i.e., $\tilde{\omega}\equiv\omega/v_{F}|{\bm{\lambda}}|$, $\tilde{k}\equiv|{\bm{k}}|/|{\bm{\lambda}}|$, $\tilde{\mu}\equiv\mu/v_{F}|{\bm{\lambda}}|$ and $\tilde{k}_{D}\equiv k_{D}/|{\bm{\lambda}}|$. The ratio of the dimensionless chemical potential and Debye wavenumber is given by $\tilde{\mu}/\tilde{k}_{D}=c\mu/v_{F}\omega_{D}\equiv\varepsilon\mu/\omega_{D}$ in terms of the typically small parameter $\varepsilon\equiv c/v_{F}$. This is rather different from the situation in relativistic plasmas, where the constituent particles are photons and electrons, both moving at the speed of light. Another important difference between the case at hand and relativistic plasmas is that there is no natural Debye frequency scale for photons like we have here. It is therefore interesting to investigate the different possible regimes set by $\tilde{\mu}/\tilde{k}_{D}$.

From now on we use these dimensionless variables and omit tildes. The dimensionless self-energy can be decomposed as

[TABLE]

with $G_{1}\equiv g_{1}^{2}|{\bm{\lambda}}|^{4}/4\rho v_{F}\varepsilon$ a new dimensionless coupling constant. Interestingly, the original coupling constant $g_{1}$ is thus increased by a factor of $1/\varepsilon\gg 1$.

After the scaling on $|{\bm{\lambda}}|$, the self-energy still depends on the direction $\hat{\lambda}^{i}$, as can be seen from the vertex in Eq. (4). To illustrate the resulting anisotropy in the self-energy, we consider the contribution of the longitudinal phonons to $\Sigma_{1}(\omega,{\bm{k}})$ in some more detail. From Eq. (6) it can be expressed as

[TABLE]

with $H_{-}(\omega,{\bm{k}}-{\bm{q}},{\bm{q}})$ a scalar function that will be specified later. The tensorial structure in Eq. (8) can be simplified by noting that the integrand is symmetric under $i\leftrightarrow j$. The integral thus has to yield a linear combination of $\delta^{ij}$ and $\hat{k}^{i}\hat{k}^{j}$. The latter leads to a contribution proportional to $(\hat{{\bm{\lambda}}}\cdot\hat{{\bm{k}}})^{2}$, whereas the former does not depend on the angle between ${\bm{\lambda}}$ and ${\bm{k}}$. Clearly the resulting expression depends on the angle between ${\bm{\lambda}}$ and ${\bm{k}}$ in a non-trivial way.

In what follows we investigate the anisotropy of the self-energy in two distinct cases: ${\bm{\lambda}}\perp{\bm{k}}$ and ${\bm{\lambda}}\parallel{\bm{k}}$. The corresponding self-energy contributions are $\Sigma^{\theta}_{n}(\omega,{\bm{k}})$, with $\theta=\perp,\parallel$ and $n=0,1$. They can be expressed as

[TABLE]

with $\eta_{0}=+,\eta_{1}=-$. The function $H_{\pm}(\omega,{\bm{x}},{\bm{y}})$ reads

[TABLE]

Finally, the functions $f_{n}^{\theta}({\bm{k}},{\bm{q}})$ in Eq. (9) capture the angle-dependence and are functions of $|{\bm{k}}|$, $|{\bm{q}}|$ and ${\bm{k}}\cdot{\bm{q}}$. We give explicit derivations and expressions for them in the Supplemental Material.

Results.— Using the expressions following from Eq. (9), we calculate the imaginary part of the self-energy analytically. Subsequently, we obtain the real part by numerically integrating the corresponding Kramers-Kronig relation. To investigate the many-body behavior resulting from the self-energy, we calculate the dimensionless spectral function $\rho^{\theta}(\omega,{\bm{k}})\equiv\rho^{\theta}_{+}(\omega,{\bm{k}})+\rho^{\theta}_{-}(\omega,{\bm{k}})$, with

[TABLE]

where $\Sigma^{\theta}_{\pm}(\omega,{\bm{k}})\equiv\Sigma^{\theta}_{0}(\omega,{\bm{k}})\pm\Sigma^{\theta}_{1}(\omega,{\bm{k}})$. The spectral function obeys the sum rule $\int_{-\infty}^{\infty}\text{d}\omega\,\rho^{\theta}(\omega,{\bm{k}})=2$, as a Weyl fermion has two degrees of freedom. We checked that the sum rule is satisfied for all numerically obtained spectral functions within $1\%$ accuracy. We plot the spectral function in Fig. 1 for several interesting cases and proceed to discuss the results presented there in detail.

We start with the case ${\bm{\lambda}}\perp{\bm{k}}$. In Fig. 1(a) we show $\rho^{\perp}(\omega,{\bm{k}})$ for $k_{D}/\mu=1/2$ and five different values of $|{\bm{k}}|\equiv|{\bm{k}}|/|{\bm{\lambda}}|$, taking small enough values such that the linear approximation around the Weyl nodes remains justified. As expected, we observe high quasiparticle peaks centered around $\omega_{{\bm{k}}}=\tilde{v}_{F}|{\bm{k}}|$, with $\tilde{v}_{F}$ a renormalized Fermi velocity. There are however, no other satellite peaks present in this case. The situation changes when we increase the ratio $k_{D}/\mu$. In Fig. 1(b) we show the case $k_{D}/\mu=2$, for which a broad peak appears around the Debye frequency. As $k_{D}/\mu$ is increased even further, c.f. Fig. 1(c), the peak centers around the new value for $k_{D}/\mu$ and becomes narrower and higher, signaling the emergence of a well-defined excitation. Interestingly, the peak is split in two for small $|{\bm{k}}|/|{\bm{\lambda}}|$ and is damped out upon increasing $|{\bm{k}}|/|{\bm{\lambda}}|$.

The situation is more pronounced when ${\bm{\lambda}}\parallel{\bm{k}}$. In this case, the behavior in the regime $k_{D}/\mu<1$ is very similar to the behavior presented for the perpendicular case in Fig. 1(a), so we refrain from showing a plot of it here. Upon increasing $k_{D}/\mu$, a distinct behavior from the case ${\bm{\lambda}}\perp{\bm{k}}$ becomes clear, as can be seen in Figs. 1(d) and (e). Especially in Fig. 1(e) it is clear that the peaks are even narrower and higher than in Fig. 1(c). Again the peaks are centered around the Debye frequency $k_{D}$ and they disperse as $|{\bm{k}}|/|{\bm{\lambda}}|$ is increased. In contrast to Fig. 1(c) the damping of the peaks as $|{\bm{k}}|/|{\bm{\lambda}}|$ is increased is much smaller, leading to a well-defined satellite peak at higher wavenumbers. The height of the satellite peak is for a large part determined by the ratio between the phonon and electron velocity $\varepsilon\equiv c/v_{F}$. As $\varepsilon$ is increased from, e.g., $1/100$ to $1/10$, the peaks become lower and broader, in addition to being located at higher frequencies. The corresponding plots are shown in the Supplemental Material.

In order to investigate the dispersive behavior, we collect the locations of the peaks for several values of $|{\bm{k}}|/|{\bm{\lambda}}|$ and plot the corresponding dispersion relation $\omega_{{\bm{k}}}$, normalized on the Debye wavenumber $k_{D}$, in Fig. 1(f). We plot the dispersion relation for different values of the ratio between the phonon and electron velocity and the Debye wavenumber $k_{D}$. Clearly the dispersion is linear for larger wavenumbers, but it seems to slightly change slope as the wavenumber goes to zero. From Fig. 1(f) it is clear that the slope of the linear part of the dispersion is set by the phonon velocity via the dimensionless parameter $\varepsilon$. Upon closer inspection of the last two terms in Eq. (10), it appears that the slope is set by $v_{F}+c=v_{F}(1+\varepsilon)$. Scaling the curves in Fig. 1(f) on $(1+\varepsilon)$ does not lead to a perfect collapse due to numerical inaccuracies.

Conclusion.— In summary, we have shown that gauge field-like electron-phonon coupling leads to novel satellite excitations in Weyl metals. The excitations have a rich phenomenology, depending on three parameters: i) the ratio between the Debye wavenumber and the chemical potential $k_{D}/\mu$, ii) the ratio between the phonon and fermion velocity $c/v_{F}$ and iii) the angle between the external wavenumber and the separation between the Weyl nodes, $\hat{{\bm{\lambda}}}\cdot\hat{{\bm{k}}}$. It is also important to note that the transverse phonons contributed significantly to all the spectral functions in Fig. 1 and excluding them leads to very different results indeed.

Regarding the potential impact of higher orders in perturbation theory, we do not expect vertex corrections to significantly change the results presented here, as Migdal’s theorem for massless fermions dictates that vertex corrections are surpressed by a factor of $c/v_{F}\ll 1$ Roy et al. (2014). Finally, we note that the emergent fermionic modes found in this work are well defined (the imaginary part of the self-energy becomes much smaller than the real part at small momenta) when $c/v_{F}\ll 1$, so in principle it is possible to implement schemes based on kinetic theory to compute the impact of these excitations on transport and optical properties Blaizot and Iancu (1993a, b); Sadooghi and Taghinavaz (2014). A combination of optical, transport, and ARPES measurements could therefore be employed to show the existence of the new fermionic excitation experimentally.

Appendix A Supplemental Material: derivation of the self-energy and additional plots

In this Supplemental Material we first derive Eq. (6). Subsequently we derive the expression in Eq. (9) and the corresponding functions $f_{n}^{\theta}({\bm{k}},{\bm{q}})$. Finally, we show some additional plots of the spectral function to support the discussion of Fig. 1 in the main text.

A.1 Matsubara sum

Upon integrating out the phonons, the full expression for the self-energy is given by

[TABLE]

The Matsubara Green’s function for a Weyl fermion with an isotropic Weyl cone and positive chirality can be written as

[TABLE]

Using this expression and the phonon propagator $D^{0}_{\alpha\alpha^{\prime}}(i\omega_{n},\bm{q})$ in Eq. (5), we rewrite the expression in Eq. (12) as

[TABLE]

in terms of $\Lambda_{\pm}(\bm{k})\equiv(\sigma^{0}\pm{\bm{\sigma}}\cdot\hat{\bm{k}})/2$ and

[TABLE]

where we converted the Matsubara sum into a contour integral in the complex plane over the Bose distribution function $N_{\text{B}}(z)\equiv(e^{\beta z}-1)^{-1}$. To make the poles more explicit, we rewrite the integrand in Eq. (15) as

[TABLE]

The integrand thus has poles at $z_{1}=b\omega^{\alpha}_{\bm{q}}$ and $z_{2}=i\omega_{m}+\mu-a\varepsilon_{\bm{k}-\bm{q}}$. Performing the integral in Eq. (15) using the residue theorem, we find

[TABLE]

As $i\omega_{m}$ is a fermionic Matsubara frequency, the second term can be simplified by noting that

[TABLE]

where $N_{\text{F}}(x)\equiv(e^{\beta x}+1)^{-1}$ is the Fermi-Dirac distribution function. Furthermore, we have the identities

[TABLE]

such that the four terms defined by Eq. (17) read

[TABLE]

We now proceed by going to the zero-temperature limit, $T=0$, and we assume a positive chemical potential, $\mu>0$. In this limit there are no thermally activated phonons and fermions. Therefore $N_{\text{B}}(\omega^{\alpha}_{\bm{q}})=0$ and $N_{\text{F}}(\varepsilon_{\bm{k}-\bm{q}}+\mu)=0$ such that $I^{\alpha}_{-+}(\bm{k},\bm{q},i\omega_{m})=0$. Performing the analytical continuation to real frequencies $i\omega_{m}\rightarrow\omega+i\eta$, with $\eta\rightarrow 0$, we obtain for the retarded self-energy

[TABLE]

which is exactly Eq. (6). Note that we shifted $\omega+\mu\rightarrow\omega$ in Eq. (21), such that all frequencies are measured with respect to the chemical potential. The obtained expression for the self-energy contains four contributions: two proportional to $\vartheta(\mu-\varepsilon_{\bm{k}-\bm{q}})$ and two with just a $1$ in the numerator. The former two are associated with the matter content and the latter two, being non-zero even at $\mu=0$, are vacuum contributions.

A.2 Sum over polarizations

The next step is to perform the sum over the polarizations by using the vertex from Eq. (4). For briefness we introduce $\hat{p}^{l}\equiv(k-q)^{l}/|\bm{k}-\bm{q}|$ and $\tilde{g}^{2}_{1}\equiv v_{F}^{2}g_{1}^{2}/4\rho$. For the longitudinal contribution to the self-energy we then find

[TABLE]

where we used $\hat{e}^{L}_{i}=\hat{q}_{i}$ and the identity $\sigma^{i}\sigma^{j}=\delta^{ij}+i\varepsilon^{ijk}\sigma_{k}$ several times. The first term in the expression above is proportional to $\sigma^{0}$ and leads to $\Sigma_{0}^{L}(\omega,\bm{k})$, whereas the other terms yield $\Sigma^{L}_{1,i}(\omega,\bm{k})\sigma^{i}$. Doing the same for the sum over the transverse modes, we find

[TABLE]

A.3 Tensorial decomposition

We now proceed with the computation of the self-energy that has been made dimensionless by scaling all parameters on $|\bm{\lambda}|$, as we discussed in the paper. Additionally, we take for simplicity $\varepsilon_{L}=\varepsilon_{T}\equiv\varepsilon$ and $k_{D}^{L}=k_{D}^{T}\equiv k_{D}$ in the following. Using Eq. (22) and Eq. (21), we find for the longitudinal contributions to the dimensionless self-energy

[TABLE]

in terms of the function $H_{\pm}(\omega,{\bm{x}},{\bm{y}})$ defined in Eq. (10). Note that $\Sigma_{0}(\omega,\bm{k})$ depends on $H_{+}(\omega,\bm{x},\bm{y})$, whereas $\Sigma_{1,i}(\omega,\bm{k})\sigma^{i}$ depends on $H_{-}(\omega,\bm{x},\bm{y})$. This is ultimately due to the fact that the last term in Eq. (21) is proportional to $\Lambda_{-}(\bm{k}-\bm{q})$, whereas the other terms are proportional $\Lambda_{+}(\bm{k}-\bm{q})$.

In order to simplify the expression Eq. (24a), we note that the integrand is symmetric in $i\leftrightarrow j$ and must therefore be a linear combination of $\delta^{ij}$ and $\hat{k^{i}}\hat{k^{j}}$. Making the ansatz $I^{ij}(k)=\alpha(k)\delta^{ij}+\beta(k)\hat{k}^{i}\hat{k}^{j}$ for the tensor resulting from the integral, one finds \alpha(k)=\big{(}\delta_{ij}-\hat{k}_{i}\hat{k}_{j}\big{)}I^{ij}(k)/2 and \beta(k)=\big{(}3\hat{k}_{i}\hat{k}_{j}-\delta_{ij}\big{)}I^{ij}(k)/2. After shifting the integration variable $\bm{q}\rightarrow\bm{q}+\bm{k}$, we find

[TABLE]

Only the terms on the second line depend on the angle between $\bm{\lambda}$ and $\bm{k}$. In the special cases $\bm{\lambda}\perp\bm{k}$ and $\bm{\lambda}\parallel\bm{k}$, we denote the self-energy as $\Sigma^{L,\perp}_{0}(\omega,\bm{k})$ and $\Sigma^{L,\perp}_{0}(\omega,\bm{k})$ and find

[TABLE]

We follow a similar strategy for Eq. (24b). The second integrand is of the same form as the integrand of Eq. (24a) and can thus be computed in the same way. The first integrand is symmetric in $i\leftrightarrow j\leftrightarrow l$. We can thus decompose the result of the corresponding integral as I^{ijl}(k)=\alpha(k)\big{(}\delta^{ij}\hat{k}^{l}+\delta^{il}\hat{k}^{j}+\delta^{jl}\hat{k}^{i}\big{)}+\beta(k)\hat{k}^{i}\hat{k}^{j}\hat{k}^{l}. The coefficients $\alpha(k)$ and $\beta(k)$ follow by contracting $I^{ijl}(k)$ with $\delta^{ij}\hat{k}^{l}$ and $\hat{k}^{i}\hat{k}^{j}\hat{k}^{l}$. They read \alpha(k)=\big{(}\delta_{ij}\hat{k}_{l}-\hat{k}_{i}\hat{k}_{j}\hat{k}_{l}\big{)}I^{ijl}(k)/2 and \beta(k)=\big{(}5\hat{k}_{i}\hat{k}_{j}\hat{k}_{l}-3\delta_{ij}\hat{k}_{l}\big{)}I^{ijl}(k)/2. Combining all contributions we write $\Sigma_{1,i}^{L,\perp}(\omega,\bm{k})\sigma^{i}=\Sigma_{1}^{L,\perp}(\omega,\bm{k})(\hat{\bm{k}}\cdot\bm{\sigma})$ when $\bm{\lambda}\perp\bm{k}$ and find

[TABLE]

In the case that $\bm{\lambda}\parallel\bm{k}$, we write $\Sigma_{1,i}^{L,\parallel}(\omega,\bm{k})\sigma^{i}=\Sigma_{1}^{L,\parallel}(\omega,\bm{k})(\hat{\bm{k}}\cdot\bm{\sigma})$, with

[TABLE]

Following the same procedure as for the longitudinal part, but starting from Eq.(23) instead of Eq.(22), we find for transverse contribution

[TABLE]

In the perpendicular case $\bm{\lambda}\perp\bm{k}$ we write similarly to the longitudinal case $\Sigma_{1,i}^{T,\perp}(\omega,\bm{k})\sigma^{i}=\Sigma_{1}^{T,\perp}(\omega,\bm{k})(\hat{\bm{k}}\cdot\bm{\sigma})$. Following the same steps that led to Eq. (27), we obtain

[TABLE]

Finally, when $\bm{\lambda}\parallel\bm{k}$ we denote $\Sigma_{1,i}^{T,\parallel}(\omega,\bm{k})\sigma^{i}=\Sigma_{1}^{T,\parallel}(\omega,\bm{k})(\hat{\bm{k}}\cdot\bm{\sigma})$ and find

[TABLE]

We have now obtained all longitudinal and transverse contributions to the self-energy in the distinct cases $\bm{\lambda}\perp\bm{k}$ and $\bm{\lambda}\parallel\bm{k}$. Because we took $\varepsilon_{L}=\varepsilon_{T}\equiv\varepsilon$ we can simplify the full expression for the self-energy by adding the longitudinal and transverse contributions and arrive at Eq. (8). Using Eqs. (26a), (26b), (29a) and (29b) we obtain the functions $f_{0}^{\perp,\parallel}(\bm{k},\bm{q})$. They are given by

[TABLE]

Similarly, we obtain from Eqs. (27),(28), (30) and (31) the functions $f_{1}^{\perp,\parallel}(\bm{k},\bm{q})$. They are given by

[TABLE]

where we introduced $x\equiv\hat{{\bm{k}}}\cdot\hat{{\bm{q}}}$.

Appendix B Additional figures

To supplement the discussion in the main text, we show in Fig. 2 plots of the spectral function $\rho^{\parallel}(\omega,{\bm{k}})$ for three different values for the ratio of the phonon and electron velocity $\varepsilon=c/v_{F}$. From the figure it becomes clear that both the height and width, and the location of the satellite excitations is determined by $\varepsilon$. The satellite peak is well-defined when $\varepsilon\ll 1$ and $|{\bm{k}}|/|{\bm{\lambda}}|$ is sufficiently small.