Optical properties and electromagnetic modes of Weyl semimetals

TL;DR

This paper provides a comprehensive theoretical analysis of electromagnetic modes in Weyl semimetals, revealing how their unique electronic structures influence optical responses and can be probed through various electromagnetic measurements.

Contribution

It introduces detailed conductivity tensors for bulk and surface states, enabling extraction of Weyl semimetal electronic properties from optical experiments.

Findings

Derived conductivity tensors for bulk and surface states.

Demonstrated how optical measurements reveal Weyl node positions and separations.

Showed the influence of Fermi energy and Fermi arcs on electromagnetic modes.

Abstract

We present systematic theoretical studies of both bulk and surface electromagnetic eigenmodes, or polaritons, in Weyl semimetals. We derive the tensors of bulk and surface conductivity taking into account all possible combinations of the optical transitions involving bulk and surface electron states. We show how information about electronic structure of Weyl semimetals, such as position and separation of Weyl nodes, Fermi energy, and Fermi arc surface states, can be unambiguously extracted from measurements of the dispersion, transmission, reflection, and polarization of electromagnetic waves.

Figures

a and 2b. The union of these $k_{y}=$ const lines and the bulk-state constant energy surface is a set of bulk and surface energy states with the same energy. In particular, at the energy equal to the Fermi energy $E_{F}$ the $k_{y}=E_{F}/(\hbar v_{F})$ line forms a Fermi arc.

a and 2b. The union of these $k_{y}=$ const lines and the bulk-state constant energy surface is a set of bulk and surface energy states with the same energy. In particular, at the energy equal to the Fermi energy $E_{F}$ the $k_{y}=E_{F}/(\hbar v_{F})$ line forms a Fermi arc.

Equations

$$\hat{H}=v_F\left(\frac{{\hat{Q}}^2-{\hslash }^{2}m(z)}{2\hslash b}{\hat{\sigma }}_x+{\hat{p}}_y{\hat{\sigma }}_y+{\hat{p}}_z{\hat{\sigma }}_z\right),$$

$$(1){\hat{Q}}^2$$

$$(2){\hat{Q}}^2$$

$$(3){\hat{Q}}^2$$

$$\hat{H}=v_F\left(\frac{{\hat{Q}}^2-{\hslash }^{2}m(z)}{2\hslash b}{\hat{\sigma }}_x+{\hat{p}}_z{\hat{\sigma }}_y-{\hat{p}}_y{\hat{\sigma }}_z\right),$$

$$\widehat{H_{\mathbf{k}}}=\hslash v_F\left(K_x\left(\mathbf{k}\right){\hat{\sigma }}_x+k_z{\hat{\sigma }}_y-k_y{\hat{\sigma }}_z\right),$$

$$(1)K_x$$

$$(2)K_x$$

$$(3)K_x$$

$$\left|\Psi_{\boldsymbol{k}}\right\rangle=\left(\begin{array}[]{c}\Psi_{1}\\ \Psi_{2}\end{array}\right)e^{i\boldsymbol{kr}-i\frac{E}{\hbar}t},$$

$$\begin{pmatrix}-k_{y}-\frac{E}{\hbar v_{F}}&K_{x}\left(\boldsymbol{k}\right)-ik_{z}\\ K_{x}\left(\boldsymbol{k}\right)+ik_{z}&k_{y}-\frac{E}{\hbar v_{F}}\end{pmatrix}\left(\begin{array}[]{c}\Psi_{1}\\ \Psi_{2}\end{array}\right)=0.,$$

$$E=s\hslash v_F\sqrt{K_x^2+k_y^2+k_z^2},$$

$$\left(\begin{array}[]{c}\Psi_{1}\\ \Psi_{2}\end{array}\right)=\frac{1}{\sqrt{2V}}\left(\begin{array}[]{c}\sqrt{1-s\cos\theta_{\boldsymbol{k}}}e^{-i\phi_{\boldsymbol{k}}}\\ s\sqrt{1+s\cos\theta_{\boldsymbol{k}}}\end{array}\right),$$

$$\left|\Psi_{B}\right\rangle=\frac{e^{ik_{x}x+ik_{y}y}}{2\sqrt{V}}\left[\left(\begin{array}[]{c}\sqrt{1-s\cos\theta_{\boldsymbol{k}}}e^{-i\phi_{\boldsymbol{k}}}\\ s\sqrt{1+s\cos\theta_{\boldsymbol{k}}}\end{array}\right)e^{ik_{z}z}-\left(\begin{array}[]{c}\sqrt{1-s\cos\theta_{\boldsymbol{k}}}e^{i\phi_{\boldsymbol{k}}}\\ s\sqrt{1+s\cos\theta_{\boldsymbol{k}}}\end{array}\right)e^{-ik_{z}z}\right],$$

$$\frac{E}{\hbar v_{F}}=-k_{y},\ \left|\Psi_{S}\right\rangle=\sqrt{\frac{2\kappa}{S}}\left(\begin{array}[]{c}1\\ 0\end{array}\right)\Theta\left(-z\right)e^{\kappa z+ik_{x}x+ik_{y}y},$$

$$\frac{E}{\hbar v_{F}}=+k_{y},\ \left|\Psi_{S}\right\rangle=\sqrt{\frac{2\kappa}{S}}\left(\begin{array}[]{c}0\\ 1\end{array}\right)\Theta\left(z\right)e^{-|\kappa|z+ik_{x}x+ik_{y}y},$$

$$\frac{E}{\hbar v_{F}}=-k_{y},\ \left|\Psi_{S}\right\rangle=N\left(\begin{array}[]{c}1\\ 0\end{array}\right)e^{ik_{x}x+ik_{y}y}\left\{\begin{array}[]{ll}e^{\int_{0}^{z}\frac{m\left(z\right)-k_{x}^{2}}{2b}dz}&{\rm for\;Hamiltonian\;1}\\ e^{\int_{0}^{z}~{}\frac{m\left(z\right)-k_{x}^{2}-k_{y}^{2}}{2b}~{}dz}&{\rm for\;Hamiltonian\;2,}\end{array}\right.$$

$$|{\Psi }_B⟩\approx \frac{e^{i{k}_{x}x+i{k}_{y}y}}{2\sqrt{V}}$$

$$\times$$

$$k_z=\sqrt{2b\sqrt{\frac{E^2}{{\hslash }^2{v}_F^2}+k_x^2}-\left(k_x^2+k_y^2+b^2\right)},\kappa =\sqrt{2b\sqrt{\frac{E^2}{{\hslash }^2{v}_F^2}+k_x^2}+\left(k_x^2+k_y^2+b^2\right)},$$

$$r=-\frac{e^{{\alpha }_{\kappa }}+e^{-i{\varphi }_{\mathbf{k}}}}{e^{{\alpha }_{\kappa }}+e^{i{\varphi }_{\mathbf{k}}}},\sinh {\alpha }_{\kappa }=\frac{\kappa }{\sqrt{\frac{E^2}{{\hslash }^2{v}_F^2}-k_y^2}},l=2i\frac{\sin {\varphi }_{\mathbf{k}}}{e^{{\alpha }_{\kappa }}+e^{i{\varphi }_{\mathbf{k}}}}.$$

$$\begin{pmatrix}-k_{y}-\frac{E}{\hbar v_{F}}&\frac{k_{x}^{2}+k_{y}^{2}-\kappa^{2}-b^{2}}{2b}-\kappa\\ \frac{k_{x}^{2}+k_{y}^{2}-\kappa^{2}-b^{2}}{2b}+\kappa&k_{y}-\frac{E}{\hbar v_{F}}\end{pmatrix}\left(\begin{array}[]{c}\Psi_{1}\\ \Psi_{2}\end{array}\right)=0$$

$$\frac{E}{\hbar v_{F}}=-k_{y},\ \left|\Psi_{S}\right\rangle=\sqrt{\frac{2}{S\left(\frac{1}{\kappa_{1}}+\frac{1}{\kappa_{2}}-\frac{4}{\kappa_{1}+\kappa_{2}}\right)}}\left(\begin{array}[]{c}1\\ 0\end{array}\right)\Theta\left(-z\right)\left(e^{\kappa_{1}z}-e^{\kappa_{2}z}\right)\cdot e^{ik_{x}x+ik_{y}y};$$

$$\frac{E}{\hbar v_{F}}=k_{y},\ \left|\Psi_{S}\right\rangle=\sqrt{\frac{2}{S\left(\frac{1}{\kappa_{1}}+\frac{1}{\kappa_{2}}-\frac{4}{\kappa_{1}+\kappa_{2}}\right)}}\left(\begin{array}[]{c}1\\ 0\end{array}\right)\Theta\left(z\right)\left(e^{-\kappa_{1}z}-e^{-\kappa_{2}z}\right)\cdot e^{ik_{x}x+ik_{y}y}.$$

$$\left|\Psi_{S;z<0}\right\rangle\propto\left(\begin{array}[]{c}1\\ 0\end{array}\right)\left(e^{\kappa_{1}z}-\zeta e^{\kappa_{2}z}\right)e^{ik_{x}x+ik_{y}y},$$

$$\left|\Psi_{S;z>0}\right\rangle\propto\left(\begin{array}[]{c}1\\ 0\end{array}\right)\frac{\kappa_{2}-\kappa_{1}}{\kappa_{2}+\sqrt{m_{\infty}}}e^{-\sqrt{m_{\infty}}z}e^{ik_{x}x+ik_{y}y},$$

$${\left({\hat{p}}_{x,y,z}-\frac{e}{c}{A}_{x,y,z}\right)}^2⟹{\hat{p}}_{x,y,z}^2+\frac{e^2}{c^2}{A}_{x,y,z}^2-\frac{e}{c}\left({\hat{p}}_{x,y,z}{A}_{x,y,z}+A_{x,y,z}{\hat{p}}_{x,y,z}\right),$$

$$-i\frac{\partial }{\partial x,\partial y,\partial z}⟹-i\frac{\partial }{\partial x,\partial y,\partial z}-\frac{e}{c\hslash }{A}_{x,y,z}.$$

$$\varphi =\frac{1}{2}\varphi (\omega )e^{-i\omega t+i\mathbf{q}\cdot \mathbf{r}}+c.c.,$$

$$\mathbf{A}=\frac{1}{2}[{\mathbf{x}}_0{A}_x(\omega )+{\mathbf{y}}_0{A}_y(\omega )+{\mathbf{z}}_0{A}_z(\omega )]e^{-i\omega t+i\mathbf{q}\cdot \mathbf{r}}+c.c.$$

$$n\left(\mathbf{r}\right)=\sum _{\alpha \beta }{n}_{\beta \alpha }\left(\mathbf{r}\right){\rho }_{\alpha \beta },\mathbf{j}\left(\mathbf{r}\right)=\sum _{\alpha \beta }{\mathbf{j}}_{\beta \alpha }\left(\mathbf{r}\right){\rho }_{\alpha \beta },$$

$$n_{\beta \alpha }={\Psi }_{\beta }^{\ast }{\Psi }_{\alpha },{\mathbf{j}}_{\beta \alpha }=\frac{1}{2}\left[{\Psi }_{\beta }^{\ast }\left(\hat{\mathbf{j}}{\Psi }_{\alpha }\right)+\left({\hat{\mathbf{j}}}^{\ast }{\Psi }_{\beta }^{\ast }\right){\Psi }_{\alpha }\right],$$

$$\mathbf{j}\left(\mathbf{r}\right)=\frac{1}{2}\sum _{\mathbf{q}}{\mathbf{j}}^{\left(\mathbf{q}\right)}{e}^{i\mathbf{qr}}+c.c.,n\left(\mathbf{r}\right)=\frac{1}{2}\sum _{\mathbf{q}}{n}^{\left(\mathbf{q}\right)}{e}^{i\mathbf{qr}}+c.c.,$$

Full text

I Introduction

Weyl semimetals (WSMs) have attracted a lot of interest as a new class of gapless three-dimensional topological materials. Their Brillouin zone contains an even number of band-touching points, or Weyl nodes, that can be described by topological invariants defined as integrals over the two-dimensional Fermi surface. For each pair of Weyl nodes, these invariants can be viewed as topological chiral charges of opposite sign of chirality volovik2007 . The electron dispersion near each Weyl node corresponds to three-dimensional massless Weyl fermions. For crystals with broken time-reversal or inversion symmetry (or both), the Weyl nodes of opposite chirality are separated in momentum space. The separation makes them stable against small perturbations and also gives rise to surface states with Fermi arcs. For reviews of WSMs discovered so far and their properties, see hosur2013 ; yan2017 ; hasan2017 ; armitage2018 ; burkov2018 ; kotov2018 .

So far, the bulk of the research has been focused on measuring and modeling the electronic structure of WSMs and topological signatures in electron transport. However, it is becoming increasingly clear that optical studies in the terahertz to mid-infrared spectral regions (e.g. shushkov2015 ) can provide a sensitive and sometimes more selective probe into the unique properties of these materials as compared to other methods. For a WSM in a magnetic field several proposals explored the signatures of the chiral anomaly in the interband optical absorption and plasmon mode properties; see e.g. the calculations of the magnetooptical conductivity in the quasiclassical limit spivak2016 ; zhou2015 ; pell2015 ; tabert2016 ; ashby2014 ; ma2015 ; gorbar2017 and the quantum-mechanical theory in a strong magnetic field ashby2013 ; long2018 . Note that these studies did not include finite separation of Weyl nodes in a microscopic Hamiltonian.

Here we study electromagnetic eigenmodes of WSMs in the presence of finite separation between Weyl nodes in momentum space and without an external magnetic field. To calculate the optical response, one needs to determine a realistic low-energy Hamiltonian that captures the essential topological structure of WSMs. While many WSMs discovered in experiment have a complicated arrangement of several pairs of Weyl nodes, essential physics and electronic properties of WSMs are already revealed in a model containing only two Weyl nodes separated in momentum space. Such models serve as a usual starting point for theoretical studies of transport and optical phenomena. Probably the simplest approach is to add a Zeeman-like constant shift term to the Hamiltonian for a Dirac semimetal, which preserves the linear form of the Hamiltonian with respect to momentum operators chen2015 . The bulk optical conductivity for this model was calculated in tabert2016-2 . In another approach, developed in zyuzin2012 and used in many optical response studies to date, a phenomenological axion $\theta$-term is introduced in the action for the electromagnetic field. This gives rise to the gyrotropic terms in the dielectric permittivity tensor and associated effects of Faraday and Kerr rotation, linear dichroism, modification of surface plasmon dispersion etc.; see e.g. kargarian2015 ; hofmann2016 ; ukhtary2017 ; kotov2018 .

In yet another approach, Okugawa and Murakami okugawa2014 derived a minimal 2x2 Hamiltonian (one conduction and one valence band) containing one parameter which describes the transition from the normal insulator to the WSM with two Weyl nodes separated in momentum space and eventually to a topological insulator in the bulk. In the WSM phase, this Hamiltonian allows for surface state solutions with Fermi arcs. Therefore, a single microscopic Hamiltonian can be used to describe optical transitions between the bulk states, surface states, and surface-to-bulk states. As a result, both bulk and surface tensors of the optical conductivity can be derived. The Hamiltonian of okugawa2014 has been recently used to develop a quantum-mechanical theory of surface plasmons (Fermi arc plasmons) and their dissipation andolina2018 .

Here we use a slightly more general Hamiltonian, which is free of certain surface state pathologies, to perform quantum-mechanical derivation of the tensors of both bulk and surface conductivity. We take into account all possible combinations of transitions between bulk and surface electron states. We then proceed to determine the properties of bulk and surface electromagnetic eigenmodes, or polaritons. We show how information about the electronic structure of WSMs, such as position and separation of Weyl nodes, Fermi energy, surface states, Fermi arcs, etc. can be extracted from the transmission, dispersion, reflection, and polarization of electromagnetic modes. We identify the most sensitive optical signatures of the electronic properties of WSMs and discuss the potential use of WSM thin films for optoelectronic applications.

Section II describes the effective Hamiltonian, or rather a family of Hamiltonians used in this study and derives the properties of corresponding bulk and surface electron states. Section III gives the classification of possible optical transitions and outlines all steps in the derivation of tensors of bulk and surface optical conductivity. The explicit expressions for the tensor elements are given in the Appendix. Section IV provides a detailed description of the electromagnetic normal modes (polaritons) in bulk WSMs. Section V provides boundary conditions which are then used in Sec. VI to calculate the reflection of incident radiation from the surface of a WSM. Section VII describes surface electromagnetic eigenmodes, i.e. surface plasmon-polaritons. Conclusions are in Sec. VIII. The Appendix contains matrix elements of the current density operator, general expressions for elements of the bulk and surface conductivity tensor, their low-frequency limit and the limit of small Weyl node separation.

II Effective Hamiltonian

In this section we describe the family of Hamiltonians that serve as a microscopic basis in this study. We derive the properties of bulk and surface electron states and use them to calculate the optical conductivity. Consider a family of Hamiltonians of the type

[TABLE]

where the function $m(z)$ takes into account that the system may be nonuniform along $z$ and, in particular, has boundaries. Here $\hat{\sigma}_{x,y,z}$ are Pauli matrices and the operator $\hat{Q}^{2}\textbf{ }$ is defined by one of the following three expressions:

[TABLE]

The first option is the original Hamiltonian in okugawa2014 .

To make the derivation of surface states more convenient okugawa2014 , we apply the unitary transformation $\hat{H}\Longrightarrow\hat{S}^{-1}\hat{H}\hat{S}$ to Eq. (1), where $\hat{S}=\frac{1}{\sqrt{2}}\left(1-i\hat{\sigma}_{x}\right).$ This gives

[TABLE]

One can check that this Hamiltonian violates time-reversal symmetry due to the term proportional to $\hat{\sigma}_{x}$. The gyrotropy axis is the $x$-axis. In $\boldsymbol{k}$-representation the Hamiltonian of Eq. (2) becomes

[TABLE]

where $K_{x}\left(\boldsymbol{k}\right)$ for the same three Hamiltonians is given by

[TABLE]

In all three cases the Weyl nodes are located at $k_{x}=\pm\sqrt{m}$ assuming that $m>0$. We have found bulk and surface eigenstates for all three Hamiltonians. Below is a summary of main results related to electron states.

II.1 Hamiltonians 1 and 2

II.1.1 Bulk states

The stationary spinor eigenstate of the Hamiltonian in Eq. (3) is

[TABLE]

where the components are determined from

[TABLE]

From Eq. (5) one can get the eigenenergy of the bulk states $E\left(\boldsymbol{k}\right)$

[TABLE]

and corresponding components of the spinor eigenstate in Eq. (4):

[TABLE]

where $\cos\theta_{\boldsymbol{k}}=\frac{k_{y}}{\sqrt{K_{x}^{2}+k_{y}^{2}+k_{z}^{2}}}$, $e^{i\phi_{\boldsymbol{k}}}=\frac{K_{x}+ik_{z}}{\sqrt{K_{x}^{2}+k_{z}^{2}}}$ ; $s=\pm 1$ denotes the conduction and valence bands, and $V$ is the quantization volume.

To visualize the dispersion of electron states, we take for simplicity $m=b^{2}$. The 3D plot for one projection of 3D dispersion of the Hamiltonian 2 is shown in Fig. 1. For small energies $|\frac{E}{\hbar v_{F}}|\ll b$ the constant energy surface consists of two disconnected spheres, each of them enclosing a corresponding Weyl point; see Fig. 2. At $|\frac{E}{\hbar v_{F}}|=\frac{b}{2}$ a separatrix isoenergy surface is a 3D “figure of eight”. For $|\frac{E}{\hbar v_{F}}|>\frac{b}{2}$ the constant energy surface is simply connected and encloses both Weyl points. Figures 2a and 2b shows contours of constant energy surfaces on the plane $k_{z}=0$ for the Hamiltonians 2 and 1, respectively. The electron dispersion is strongly anisotropic. This will result in different values for the diagonal elements of the bulk dielectric permittivity tensor, as in two-axial crystals. The dotted circle in Fig. 2a is the boundary of a region that contains surface states. For Hamiltonian 1 in Fig. 2b the surface states exist between the dotted lines.

II.1.2 Reflection from the boundary. Surface states and Fermi arcs

Following okugawa2014 , we define the boundary as a jump in the parameter $m$, so that $m=b^{2}$ inside the WSM and $m=-m_{\infty}$ outside. Then Eqs. (3) and (5) will contain the parameter m as a function of the coordinate $r_{j}$ orthogonal to the boundary, and the corresponding component of the quasimomentum $k_{j}$ is replaced by $k_{j}\Longrightarrow-i\frac{\partial}{\partial r_{j}}$.

For the boundary parallel to the gyrotropy axis x, we assume that it coincides with the surface $\textit{z}=0$ and the WSM fills the halfspace $z<0$. In this case $m=b^{2}$ for $z<0$ and $m=-m_{\infty}$, $m_{\infty}\rightarrow\infty$ for $z>0$.

For Hamiltonian 3, the Schrödinger equation given by Eq. (5) is a 4th order differential equation, since its matrix elements contain $\frac{\partial^{2}}{\partial z^{2}}$ . For Hamiltonians 1 and 2 we get a 2nd order set of equations. The velocity operator $\hat{v}_{z}=\frac{i}{\hbar}\left[H,z\right]$ for Hamiltonian 3 is $\hat{v}_{z}=-i\frac{v_{F}}{b}\hat{\sigma}_{x}\frac{\partial}{\partial z}+v_{F}\hat{\sigma}_{y}$, i.e. it depends on the coordinate derivative. In contrast, the velocity operator $\hat{v}_{z}=v_{F}\hat{\sigma}_{y}$ for Hamiltonians 1 and 2 does not depend on the coordinate derivative. Therefore, for Hamiltonian 3 at $z=0$, the continuity of both the eigenstate and its derivative is required, whereas one only needs the continuity of the eigenstates for Hamiltonians 1 and 2.

Using Eq. (5) one can find that the eigenstate of Hamiltonians 1 and 2 in the region $\textit{z}>0$ at $m_{\infty}\rightarrow\infty$ is \left|\Psi_{\infty}\right\rangle\propto\left(\begin{array}[]{c}1\\ 0\end{array}\right)e^{ik_{x}x+ik_{y}y-\frac{m_{\infty}}{2b}z}. In the region $\textit{z}<0$ we take the eigenstate $\left|\Psi_{B}\right\rangle$ which is given by Eq. (7). Stitching together these two eigenstates $\left|\Psi_{\infty}\right\rangle$ and $\left|\Psi_{B}\right\rangle$ at the boundary yields the following expression for the bulk state:

[TABLE]

where the quantization volume is limited from one side by the $\textit{z}=0$ plane. The eigenenergy is still given by Eq. (6), and the angles $\theta_{\boldsymbol{k}}$ and $\phi_{\boldsymbol{k}}$ are defined below Eq. (7).

If $\frac{E^{2}}{\hbar^{2}v_{F}^{2}}<k_{y}^{2}+K_{x}^{2}$ the value of $k_{z}$ in Eq. (6) is imaginary: $k_{z}=\pm i\kappa$. In order to connect the eigenstate \left|\Psi_{\infty}\right\rangle\propto\left(\begin{array}[]{c}1\\ 0\end{array}\right) in $\textit{z}>0$ with the eigenstate localized at $\textit{z}<0$ which is $e^{\kappa z}$ , the localized eigenstate should be also a spinor \left(\begin{array}[]{c}1\\ 0\end{array}\right). After replacing $k_{z}\Rightarrow-i\kappa$ in Eq. (5), we obtain the following eigenenergies and eigenvectors for surface states in the limit $m_{\infty}\rightarrow\infty$:

[TABLE]

where $\Theta$ is a step function, S is the quantization area, $\kappa=-K_{x}>0$. For Hamiltonian 2 the surface states exist inside a dashed circle $b^{2}>k_{x}^{2}+k_{y}^{2}$ in Fig. 2a. For Hamiltonian 1 the surface states exist in the region $b^{2}>k_{x}^{2}$ in Fig. 2b.

If a WSM occupies the region $\textit{z}>0$, instead of Eqs. (9) we obtain

[TABLE]

where $\kappa=+K_{x}<0$. Equations (9),(10) can be easily generalized to the case of a parameter $m(z)$ which varies continuously between the values $b^{2}$ and $-m_{\infty}$ okugawa2014 . For example, instead of Eqs. (9) we get

[TABLE]

where $N$ is a normalization factor.

Note that the constant surface energy lines $k_{y}=$ const are tangent to the points where the bulk-state constant energy surface intersects the boundary of the surface states, shown as dotted lines in Fig. 2a and 2b. The union of these $k_{y}=$ const lines and the bulk-state constant energy surface is a set of bulk and surface energy states with the same energy. In particular, at the energy equal to the Fermi energy $E_{F}$ the $k_{y}=E_{F}/(\hbar v_{F})$ line forms a Fermi arc.

II.2 Hamiltonian 3

For a 4th order set of differential equations the construction of electron states including their interaction with a boundary is more complicated.

First, we use Eq. (6) to find the value of $k_{z}$ for given $k_{x,y}$ and $E$. Consider the parameter range $m\leq b^{2}$, including both positive and negative values of $m$. If $\frac{E^{2}}{\hbar^{2}v_{F}^{2}}>k_{y}^{2}+\frac{\left(k_{x}^{2}+k_{y}^{2}-m\right)^{2}}{4b^{2}}$, one always has two real solutions $k_{z1}=-k_{z2}>0$ together with two imaginary solutions corresponding to evanescent states: $k_{z3,4}=i\kappa_{3,4}$, where $0<\kappa_{3}=-\kappa_{4}$ . If $\frac{E^{2}}{\hbar^{2}v_{F}^{2}}<k_{y}^{2}+\frac{\left(k_{x}^{2}+k_{y}^{2}-m\right)^{2}}{4b^{2}}$ , all four solutions are imaginary and correspond to evanescent states: $k_{z1,2,3,4}=i\kappa_{1,2,3,4}$, where $0<\kappa_{1}=-\kappa_{3}$, $0<\kappa_{2}=-\kappa_{4}$. In the region $z>0$ (i.e. outside the sample, where $m=-m_{\infty}$ ) it is reasonable to take the solution as a superposition of two localized modes $e^{-|\kappa_{3,4}|z}$. In this case for $z<0$, i.e. inside the sample where $m=b^{2}$, there could be two options:

(i) A superposition of two counterpropagating waves with quasimomenta $k_{z1}=-k_{z2}$ together with a localized wave $e^{\kappa_{3}z}$. The localized solution cannot be discarded, since without it the number of constants becomes smaller than the number of the boundary conditions.

(ii) A superposition of two localized waves i.e. the surface state. In this option the number of constants is always smaller than the number of the boundary conditions, so such a state can exist only at certain values of energy.

The procedure of stitching the spinor components and their derivatives is simplified if $m_{\infty}\rightarrow\infty$ since in this limit the continuity of the derivative is equivalent to setting both components of a spinor $\Psi_{1,2}$ equal to zero in the cross section $z=0$.

II.2.1 Bulk states near the boundary

In case (i) we obtain

[TABLE]

where

[TABLE]

[TABLE]

Clearly, $|r|^{2}=1$, which corresponds, as expected, to the total reflection from the boundary. The quantization volume in Eq. (18) is chosen in such a way that its length along the $z$ axis is much larger than $k_{z}^{-1}>\kappa^{-1}$. Therefore, the last term in the brackets in Eq. (18) is unimportant in a sense that it does not affect the eigenstate normalization or the matrix elements.

II.2.2 Surface states

To construct the surface states (option (ii)) it is convenient to to go back to Eq. (5), use $m=b^{2}$, and make the substitution $k_{z}=-i\kappa$:

[TABLE]

Consider the solution of Eq. (19), corresponding to different positive values of $\kappa_{1,2}$ but the same spinor constant \left(\begin{array}[]{c}a\\ b\end{array}\right). One can build a nontrivial localized solution \left|\Psi_{S}\right\rangle\propto\left(\begin{array}[]{c}a\\ b\end{array}\right)\Theta\left(-z\right)\left(e^{\kappa_{1}z}-e^{\kappa_{2}z}\right), which corresponds to the null boundary conditions at the surface $\textit{z}=0$. Such a solution of Eq. (19) is possible under the following conditions:

$-k_{y}-\frac{E}{\hbar v_{F}}=\frac{k_{x}^{2}+k_{y}^{2}-\kappa^{2}-b^{2}}{2b}+\kappa=0$, or $k_{y}-\frac{E}{\hbar v_{F}}=\frac{k_{x}^{2}+k_{y}^{2}-\kappa^{2}-b^{2}}{2b}-\kappa=0$, or $k_{y}-\frac{E}{\hbar v_{F}}=\frac{k_{x}^{2}+k_{y}^{2}-\kappa^{2}-b^{2}}{2b}-\kappa=0$, where\left(\begin{array}[]{c}a\\ b\end{array}\right)=\left(\begin{array}[]{c}1\\ 0\end{array}\right) or \left(\begin{array}[]{c}a\\ b\end{array}\right)=\left(\begin{array}[]{c}0\\ 1\end{array}\right)respectively. It is easy to see that the first option corresponds to the surface state when the WSM occupies the halfspace $z<0$, whereas the second option corresponds to the WSM in the region $z>0$, since in this case the values of $\kappa_{1,2}$ are negative. The resulting states are as follows.

(i) WSM in $z<0$:

[TABLE]

(ii) WSM in $z>0$:

[TABLE]

Here $\kappa_{1,2}=b\mp\sqrt{k_{x}^{2}+k_{y}^{2}}$ .

In the region $b^{2}<k_{x}^{2}+k_{y}^{2}$ there is only one localized evanescent solution for any fixed value of energy, which is not enough to satisfy the boundary conditions. Therefore, the region $b^{2}>k_{x}^{2}+k_{y}^{2}$, where the surface states exist, is the same in the models described by the Hamiltonian 2 and Hamiltonian 3 (see the dotted circle in Fig. 2a).

Taking into account a finite value of $m_{\infty}$ modifies the above expression, but their general structure remains the same. For example, when a WSM fills the halfspace $z<0$, then the eigenstate in Eq. (20) is replaced by

[TABLE]

[TABLE]

where $\displaystyle\zeta=\frac{\kappa_{1}+\sqrt{m_{\infty}}}{\kappa_{2}+\sqrt{m_{\infty}}}$.

II.3 The boundary orthogonal to the gyrotropy axis

Any Hamiltonian, 1, 2, or 3, contains the second derivative $\frac{\partial^{2}}{\partial x^{2}}$. Therefore, the analysis of the bulk and surface states near the boundary orthogonal to the gyrotropy axis is similar to the one for the boundary parallel to the gyrotropy axis when the Hamiltonian contains the second derivative $\frac{\partial^{2}}{\partial z^{2}}$. Repeating the same arguments as in the previous section, we obtain that the orthogonal boundary is trivial and does not contain surface states.

II.4 Comparison of Hamiltonians 1, 2, and 3

The only principal difference between the eigenstates of the effective Hamiltonians considered above is the region where the surface states exist. Such a region is determined by the inequality $b>\sqrt{k_{x}^{2}+k_{y}^{2}}$ for Hamiltonians 2 and 3, and the inequality $b>|k_{x}|$ for Hamiltonian 1. The latter condition leads to an infinite density of surface states, which is unphysical and would have to be restricted artificially. Therefore, it is better to work with Hamiltonian 2 or 3. Hamiltonian 2 leads to a simpler $z$-component of the velocity operator: $\hat{v}_{z}=v_{F}\hat{\sigma}_{y}$ instead of $\hat{v}_{z}=-i\frac{v_{F}}{b}\hat{\sigma}_{x}\frac{\partial}{\partial z}+v_{F}\hat{\sigma}_{y}$, which corresponds to Hamiltonian 3. The velocity operator of Hamiltonian 2 makes calculations of the surface current easier without losing any essential physics. Therefore, we will use Hamiltonian 2 for subsequent calculations of the optical properties.

III Optical transitions and the tensors of bulk and surface conductivity

In the presence of external fields one should replace $\boldsymbol{\hat{\boldsymbol{p}}}\Longrightarrow\boldsymbol{\hat{p}}-\frac{e}{c}\boldsymbol{A}$, and also add the electrostatic potential $\hat{H}\Longrightarrow\hat{H}+e\varphi\hat{1}$ in Eq. (2). Particles are assumed to have charge $e$ where $-e$ is the magnitude of the electron charge. If the potential has a coordinate dependence $\boldsymbol{A}(\boldsymbol{r})$ we assume symmetrized operators

[TABLE]

and in the expressions for the velocity operator we need to replace

[TABLE]

Throughout the paper, we will consider the potentials corresponding to a monochromatic electromagnetic field propagating in the arbitrary direction $\boldsymbol{r}$ with angular frequency $\omega$ and wavevector $\boldsymbol{q}$, i.e.

[TABLE]

Bulk-to-bulk and surface-to-surface transitions contribute to the bulk and surface conductivity tensors, respectively. The contributions are detailed in the Appendix. Surface-to-bulk transitions contribute to the surface conductivity tensor only. They have to be handled with more care, as we briefly describe below.

Generally, the electron and current densities expressed in terms of the density matrix are given by

[TABLE]

[TABLE]

where $\boldsymbol{\hat{j}}=e\boldsymbol{\hat{v}}\textit{.}$

The Fourier harmonics of the the electron and current densities are

[TABLE]

where

[TABLE]

For their matrix elements we have

[TABLE]

where

[TABLE]

To find the current without the effect of a boundary we can use the states given by Eq. (7).

Now consider the states near the surface. We will denote the bulk states with latin indices and surface states with greek ones. For convenience we rewrite them, having in mind an upper boundary $z=0$ with the WSM located at $z<0:$

[TABLE]

where $E_{m}=s\hbar v_{F}\sqrt{\left(\frac{k_{x}^{2}+k_{y}^{2}-b^{2}}{2b}\right)^{2}+k_{y}^{2}+k_{z}^{2}}$ is the eigenenergy, $s=\pm 1$ is the band index, the values $k_{x,y}$ can be of either sign whereas $k_{z}>0$; $\cos\theta_{\boldsymbol{\boldsymbol{k}_{\parallel}}}=\frac{k_{z}}{\frac{|E|}{\hbar v_{F}}}$.

[TABLE]

where $S$ is the area; the energy of the state is $E_{\alpha}=-\hbar v_{F}k_{y}$, $\kappa=\frac{b^{2}-k_{x}^{2}-k_{y}^{2}}{2b}$, $\sqrt{k_{x}^{2}+k_{y}^{2}}<b$.

Let us limit the surface states by the condition $\kappa>\kappa_{min}$, where the latter could be a typical scattering length $\sim\kappa_{min}^{-1}$. We will assume that $\kappa_{min}^{-1}$ is much smaller than $L$, which enters the quantization volume $V=SL$ in Eq. (29). When we calculate the matrix elements of the interaction Hamiltonian in the von Neumann equation, the matrix elements $V_{mn}^{\left(int\right)}$ ,$V_{\alpha\beta}^{\left(int\right)}$ and $V_{m\alpha}^{\left(int\right)}$ have no peculiarities: the integration is carried out over the whole volume. However when we calculate the matrix elements of the density and current, and if at least one of the indices belongs to the surface state, we will perform the integration over $dz$:

[TABLE]

[TABLE]

These quantities will depend only on x and y, and therefore determine the surface current and density.

The following current component is nontrivial: $\sum_{\alpha\beta}\left(j_{z}\right)_{\beta\alpha}\rho_{\alpha\beta}+\sum_{m\alpha}\left(j_{z}\right)_{m\alpha}\rho_{\alpha m}$. It determines the polarization of a thin double layer:

[TABLE]

This layer radiates, but not normally to the layer, and it cannot affect the surface density of carriers.

With properly defined matrix elements of the current and density the continuity equation is satisfied automatically, so we can consider the volume current flowing toward the boundary $\left(\sum_{mn}\left(j_{z}\right)_{nm}\rho_{mn}\right)_{z=0}$ as a source in the surface continuity equation.

III.1 Tensors of bulk and surface conductivity

The matrix elements of the Fourier components of the current density operator are evaluated in Appendix A. After evaluating them, in Appendix B and C we used the Kubo-Greenwood formula to calculate the bulk and surface conductivity tensors, respectively; e.g.

[TABLE]

for the bulk conductivity, where $g=2$ is the spin degeneracy factor and $\alpha,\beta$ denote Cartesian coordinate components. The surface conductivity tensor has a similar structure, but the contribution is summed over surface-to-surface and surface-to-bulk transitions, and the normalization is over the surface area $S$ instead of a volume $V$. Both interband and intraband transitions are included. Three-dimensional integrals over electron momenta in Appendix B and C cannot be evaluated analytically, except limiting cases of small frequencies or small $b$ (see Appendix D and E). Therefore, integrals were calculated numerically at zero temperature for all plots below.

The bulk (3D) conductivity tensor due to low-energy electrons near Weyl points is

[TABLE]

where $\sigma_{zy}^{B}=-\sigma_{yz}^{B}$. The surface conductivity tensor at $z=0$ has a similar structure, with superscript $B$ replaced by $S$ and $\sigma_{zy}^{S}=-\sigma_{yz}^{S}$.

The background bulk dielectric tensor in the most general form which corresponds to the one for a two-axial nongyrotropic crystal is

[TABLE]

so that the total dielectric permittivity tensor is

[TABLE]

where

[TABLE]

Note that for Hamiltonian 3 we would have $\sigma_{yy}^{B}=\sigma_{zz}^{B}$, whereas for Hamiltonian 2 (used in all calculations of the conductivity tensors in this paper) we have $\sigma_{yy}^{B}\neq\sigma_{zz}^{B}$. Therefore, even if the background dielectric tensor is isotropic, the contribution of massless Weyl electrons will create a two-axial anisotropy. In the numerical plots below we will take an isotropic background dielectric tensor and neglect its frequency dependence at low frequencies, $\varepsilon_{xx}^{(0)}=\varepsilon_{yy}^{(0)}=\varepsilon_{zz}^{(0)}=10$, so that all nontrivial effects of anisotropy and gyrotropy are due to Weyl fermions.

The salient feature of both bulk and surface conductivity tensor is the presence of nonzero off-diagonal (gyrotropic) components due to time-reversal symmetry breaking in the Hamiltonian. These terms originate from the finite separation of the Weyl nodes in momentum space and the existence of surface states (Fermi arcs). The gyrotropic effects in the propagation, reflection, and transmission of bulk and surface modes can serve as a definitive diagnostic of Weyl nodes, surface states, and Fermi surface. They could also find applications in optoelectronic devices such as Faraday isolators, modulators etc.

Figures 3-6 show spectra of $\varepsilon_{xx}(\omega)$, $\varepsilon_{yy}(\omega)$, $\varepsilon_{zz}(\omega)$, and $g(\omega)$ for several values of the Fermi momentum $k_{F}$ (at zero temperature), when the Weyl node separation $2\hbar v_{F}b=200$ meV. The characteristic feature in all plots is strong absorption and dispersion at the onset of interband transitions, when $\omega=2v_{F}k_{F}$. Another common feature is a Drude-like increase in the absolute value of all tensor components at low frequencies. Indeed, as shown in Appendix D, in the limit $\omega\ll v_{F}k_{F}\ll v_{F}b$ when only the intraband transitions in the vicinity of each Weyl point are important, the off-diagonal components are equal to zero and the diagonal conductivity components are reduced to the same Drude form:

[TABLE]

Note an absorption peak at $\omega=100$ meV at low Fermi momenta, which corresponds to a Van Hove singularity at the interband transitions between saddle points of conduction and valence bands at $k=0$, i.e. in the middle between the Weyl points.

Note also that diagonal and off-diagonal parts of the conductivity tensor are of the same order at low frequencies comparable to the Weyl node separation, which indicates that gyrotropic effects should be quite prominent.

All figures in this paper are plotted for a relatively high dephasing rate $\gamma=10$ meV, which smoothes out all spectral features and introduces strong losses for electromagnetic eigenmodes even below the interband transition edge. The dephasing rate originates from electron scattering and obviously depends on the temperature and material quality in realistic materials. Its derivation is beyond the scope of the present paper.

IV Bulk polaritons in Weyl semimetals

Consider first the propagation of plane monochromatic waves in a bulk Weyl semimetal. For complex amplitudes of the electric field and induction, $(\boldsymbol{D,E})e^{i\boldsymbol{kr}-i\omega t}$, where $\boldsymbol{D}=\hat{\varepsilon}\boldsymbol{E}$ and $\hat{\varepsilon}$ is a bulk dielectric tensor, Maxwell’s equations give $\boldsymbol{n\cdot D}=0$, where $\boldsymbol{n}=\frac{c\boldsymbol{k}}{\omega}$. The resulting dispersion equations are

[TABLE]

or

[TABLE]

The structure of these equations indicate strongly anisotropic and gyrotropic properties of bulk polaritons. The dispersion is drastically different for normal modes propagating perpendicular to the $x$-axis and to the $y$-axis. For each direction, there are furthermore two normal modes with different refractive indices. We will consider each case separately.

IV.1 Propagation perpendicular to the anisotropy x-axis

In this case we have $n_{x}=0$, $n^{2}=n_{y}^{2}+n_{z}^{2}$, $n_{z}=n\cos\theta$, $n_{y}=n\sin\theta$, where $\theta$ is the angle between the wave vector and $z$-axis. From Eqs. (41) we obtain two normal modes that can be called an ordinary (O) and extraordinary (X) wave. An O-wave has an electric field along $x$ and the refractive index

[TABLE]

Therefore, its dispersion and absorption are completely described by the spectrum of $\varepsilon_{xx}(\omega)$. As shown in Fig. 7, at low frequencies the O-mode experiences strong metallic absorption and at $\omega=2E_{F}=160$ meV there is an onset of interband transitions.

An X-wave have an electric field in the $(y,z)$ plane and the refractive index showing strong $\theta$-dependence and resonances:

[TABLE]

For normal incidence $\theta=0$,

[TABLE]

It is obvious from Eq. (43) that the refractive index for an X-wave is strongly enhanced (singular in the absence of losses) when

[TABLE]

which corresponds to the bulk plasmon excitation. Indeed, from Maxwell’s equations in the Coulomb gauge one can show that $|\frac{1}{c}\frac{\partial A}{\partial t}|/|\nabla\varphi|\sim|\frac{\omega^{2}}{\omega^{2}-c^{2}k^{2}}||\frac{j_{\perp}}{j_{\parallel}}|$, where $\mathbold{j}=\mathbold{j}_{\perp}+\mathbold{j}_{\parallel}$, $\nabla\times\mathbold{j}_{\parallel}=0$, $\nabla\cdot\mathbold{j}_{\perp}=0$. Therefore, if $|\mathbold{j}_{\perp}|\sim|\mathbold{j}_{\parallel}|$, which corresponds to a general oblique propagation in an anisotropic medium, the wave is quasi-electrostatic at $n^{2}\gg 1$. Eq. (45) corresponds to the condition $\mathbold{n\cdot D}=0$ for $\mathbold{E}=-\nabla\varphi\parallel\mathbold{n}$ . If $\varepsilon_{yy}=\varepsilon_{zz}=\varepsilon_{\perp}$ the dispersion equation for a plasmon propagating in the plane orthogonal to the $x$-axis has a simple form $\varepsilon_{\perp}=0$.

Figure 8 shows real and imaginary parts of the refractive index $n_{X}$ of an X-wave as a function of frequency for different values of the propagation angle $\theta$. Near the bulk plasmon resonance, i.e. around 100 meV for normal incidence, the value of $n_{X}^{2}$ becomes negative in the absence of losses according to Eq. (44)This corresponds to a non-propagating photonic gap. Since we include significant loss rate $\gamma=10$ meV in all simulations, the real part of $n_{X}$ does not go all the way to zero, but there is a strong absorption peak in the imaginary part of $n_{X}$. We will later see that this spectral region leads to a telltale change of phase in reflection. The second feature in all plots is an onset of interband transitions at $2E_{F}=160$ meV.

The real part of the bulk plasmon resonance frequency at normal incidence as a function of the Fermi energy is shown in Fig. 9. Note that according to Eq. (44) the magnitude of the refractive index at frequencies around plasmon resonance is determined by the value of the off-diagonal component of the dielectric tensor $g$. Therefore, measurements of the transmission and reflection provide a sensitive measure of the Weyl node separation.

The same is true about the polarization effects. From the third row of Eqs. (41) one can get the expression for the polarization coefficient:

[TABLE]

Substituting Eq. (43) into Eq. (46) we get

[TABLE]

At the resonant plasmon frequency defined by $\cos^{2}\theta\varepsilon_{zz}+\sin^{2}\theta\varepsilon_{yy}=0$ we obtain $K_{X}=\frac{1}{\tan\theta}$, which is expected. If we set $\theta=0$, which corresponds to normal incidence, $K_{X}=\frac{ig}{\varepsilon_{zz}}$, i.e. again proportional to $g$. In this case, the plasmon frequency is given by $\varepsilon_{zz}=0$, and $K_{X}\rightarrow\infty$ in the absence of losses. If $\varepsilon_{yy}=\varepsilon_{zz}=\varepsilon_{\perp}$, Eq. (47) gives

[TABLE]

For an isotropic medium, when $g^{2}=0$, the last expression gives $K_{X}=-\tan\theta$, as it should be for a transverse wave in an isotropic medium.

IV.2 Propagation transverse to the $y$-axis

In this case $n_{y}=0$, $n^{2}=n_{x}^{2}+n_{z}^{2}$, $n_{x}=n\cos\phi$, $n_{z}=n\sin\phi$;

[TABLE]

[TABLE]

Note that the solution of Eq. (50) at $\phi=\frac{\pi}{2}$ corresponds to the normal incidence propagation along $z$ and therefore should coincide with Eqs. (42), (43) at $\theta=0$. Indeed, from Eq. (50) for $\phi=\frac{\pi}{2}$ we obtain

[TABLE]

from which $n_{O}^{2}=\varepsilon_{xx}$, $n_{X}^{2}=\varepsilon_{yy}-\frac{g^{2}}{\varepsilon_{zz}}$, as expected.

The case $n^{2}\rightarrow\infty$ in the absence of losses, when

[TABLE]

corresponds to the condition $\mathbold{n\cdot D}=0$ where $\mathbold{E}=-\nabla\varphi\parallel\mathbold{n}$. From Eq. (50) we obtain

[TABLE]

In Eq. (53) the signs $\pm$ are chosen for $n_{O,X}^{2}$ according to the limiting case $\phi=\frac{\pi}{2}$.

For the propagation along the $x$-axis of anisotropy, when $\phi=0$, Eq. (50) gives

[TABLE]

Note that the $x$-axis is also a gyrotropy axis related to the Weyl node separation along $x$. Therefore, the propagation along $x$ is similar to the Faraday geometry in a magnetic field. In our case the normal modes are elliptically polarized, and an incident linearly polarized wave experiences Faraday rotation and gains ellipticity after traversing a sample in $x$-direction. To quantify the effect, Fig. 10 shows the polarization coefficient $K_{X}=E_{z}/E_{y}$ after traversing a 1-$\mu$m thick film for a wave initially linearly polarized in $y$-direction. The real part of $K_{X}$ is a measure of the polarization rotation whereas its imaginary part is a measure of ellipticity. Clearly, a rotation by $\sim\pi/2$ by very thin (0.5-1 $\mu$m) Weyl semimetal films is possible at frequencies near the interband absorption edge. This is a giant Faraday rotation, comparable to the one observed at THz frequencies in narrow-gap semiconductors in the vicinity of a cyclotron resonance in Tesla-strength magnetic fields; see e.g. arikawa for the review. Note that in our case no magnetic field is needed and the effect is controlled by the Weyl node separation and by the Fermi level. Previously Faraday rotation and nonreciprocity in light propagation associated with it was studied in kargarian2015 ; kotov2018 using the model with an axion $\theta$-term in the electromagnetic field action.

IV.3 Oblique propagation of bulk polaritons

In the general case the direction of the wave vector is determined by two angles $\theta$ and $\phi$:

[TABLE]

The general expression for $n_{O,X}^{2}$ is quite cumbersome. At the same time, in the particular case of $\varepsilon_{yy}=\varepsilon_{zz}=\varepsilon_{\perp}$, the result should not depend on the angle $\theta$ and should coincide with the one for a magnetized plasma:

[TABLE]

The condition $\mathbold{n\cdot D}=0$ at $\mathbold{E}=-\nabla\varphi\parallel\mathbold{n}$ in the case of an oblique propagation gives

[TABLE]

Therefore, Eq. (56) determines the frequencies of bulk plasmons in the general case. Under the condition $\varepsilon_{yy}=\varepsilon_{zz}=\varepsilon_{\perp}$ the plasmon dispersion equation takes a form similar to plasmons in a magnetized plasma:

[TABLE]

V Boundary conditions

So far we considered propagation and transmission of electromagnetic waves in bulk samples. Now we turn to effects of reflection and surface wave propagation that are equally sensitive to the electronic structure of WSMs. Moreover, in many situations they are easier to observe than bulk propagation effects.

We start with the derivation of the boundary conditions at $z=0$ surface. Assume that there is an isotropic dielectric medium with dielectric constant $n_{up}^{2}=\varepsilon_{up}$ above a WSM. The boundary conditions include:

(i) Gauss’ law for the normal components of the electric induction vector:

[TABLE]

where $\rho^{S}$, $j_{x}^{S}$ and $j_{y}^{S}$ are the surface charge and components of the surface current that are connected by the continuity equation. For the wave field we have $\frac{\partial}{\partial x,\partial y}\rightarrow ik_{x,y}$.

(ii) Equations for the magnetic field components:

[TABLE]

Due to the presence of the components of the surface conductivity $\sigma_{zz}^{S}$ and $\sigma_{zy}^{S}=-\sigma_{yz}^{S}$ a surface dipole layer is formed at the boundary between the two media. Its dipole moment is

[TABLE]

Note that when dealing with a surface response, we will always choose the fields at $z=-0$ in Eq. (62) and similar relationships. The presence of the dipole layer changes the boundary conditions for the tangential field components of $\mathbold{E}$. Consider Maxwell’s equations

[TABLE]

For convenience, let’s assume that the dipole layer has a small but finite thickness $L$:

[TABLE]

Using $\frac{\partial}{\partial x,\partial y}\rightarrow ik_{x,y}$ and integrating $\int_{-\frac{L}{2}}^{\frac{L}{2}}\ldots dz$ , we obtain

[TABLE]

We neglect the integral over the magnetic field components assuming that $\frac{\omega}{c}L\rightarrow 0$. Next we use Gauss’ law under the condition $|k_{x,y}|L\rightarrow 0$, which will yield in the region of the dipole layer:

[TABLE]

Here $P_{z}$ is a component of the volume polarization whereas $p_{z}$ describes the distribution of the polarization in the dipole layer, so that

[TABLE]

Substituting $E_{z}=-4\pi\left(P_{z}+p_{z}\right)$ into Eq. (63) and integrating over $dz$ at $|k_{x,y}|L\rightarrow 0$ and finite $P_{z}$ , we obtain

[TABLE]

The boundary condition Eq. (64) looks unusual but it can be easily deduced from the radiation field of an individual dipole.

Figures 11-14 show spectra of the surface conductivity components for different values of the Fermi momentum. Note that the surface conductivity in Gaussian units has a dimension of velocity and its value is normalized by $e^{2}/(2\pi\hbar)\simeq 3.5\times 10^{7}$ cm/s in all plots. In contrast with the bulk conductivity, the surface conductivity had a Drude-like behavior at low frequencies only for the $yy$-component because of the surface state dispersion $E=-\hbar v_{F}k_{y}$. The surface optical response decreases with increasing Fermi energy and vanishes when all surface states within $k_{x}^{2}+k_{y}^{2}<b^{2}$ are occupied.

VI Reflection from the surface of a Weyl semimetal

Consider radiation incident from a medium with refractive index $n_{up}$ on a WSM at an angle $\theta$ between the wavevector of the wave and the normal to a WSM. For simplicity consider the propagation transverse to the $x$-axis. The reflection spectra provide information about both bulk and surface conductivity components. Here we will pay particular attention to the case when the contribution of the surface states becomes significant or dominant, thus allowing one to probe surface states by optical means.

VI.1 Reflection with excitation of an $O$-mode

In this geometry, the complex amplitudes of the electric field of the incident $E_{1}$, reflected $E_{2}$, and transmitted $E_{O}$ wave are parallel to the $x$-axis. The refractive index of the transmitted wave is $n_{O}^{2}=\varepsilon_{xx}=\varepsilon_{xx}^{\left(0\right)}+i\frac{4\pi}{\omega}\sigma_{xx}^{B}$ (see Eq. (42)).

Applying Maxwell’s equations with standard boundary conditions including the surface current, we arrive at

[TABLE]

where $n_{up}\sin\theta=n_{O}\sin\theta_{O}$. Assuming $\sigma_{xx}^{S}=0$ we obtain $R=\frac{E_{2}}{E_{1}}=\frac{\cos\theta n_{up}-\cos\theta_{O}n_{O}}{\cos\theta_{O}n_{O}+\cos\theta n_{up}}$, which is a standard Fresnel formula.

For the same magnitude of $\sigma_{xx}^{S}$, the relative contribution of surface states to the reflected field depends on the parameter $\frac{|\varepsilon_{xx}^{\left(0\right)}|}{4\pi|\sigma_{xx}^{B}|/\omega}$. If $\frac{\omega|\varepsilon_{xx}^{\left(0\right)}|}{4\pi|\sigma_{xx}^{B}|}\gg 1$, the relative contribution of surface states is determined by the expression: $\frac{2\omega|\sigma_{xx}^{S}|/c}{|\sigma_{xx}^{B}|/|\varepsilon_{xx}^{\left(0\right)}|}$. If $\frac{\omega|\varepsilon_{xx}^{\left(0\right)}|}{4\pi|\sigma_{xx}^{B}|}\ll 1$, one needs to evaluate the ratio $\frac{2\sqrt{\pi}\sigma_{xx}^{S}/c}{\sqrt{\sigma_{xx}^{B}/\omega}}$.

VI.2 Reflection with excitation of an $X$-mode

In this geometry, the complex Fourier harmonics for the incident and reflected waves are

[TABLE]

The transmitted wave is

[TABLE]

where $n_{X}^{2}$ and $K_{X}$ are given by Eqs. (43) and (46), in which one should substitute $\theta\rightarrow\theta_{X}$. The corresponding complex amplitudes of the magnetic field are $B_{1x}=\frac{n_{up}}{\cos\theta}E_{1}$, $B_{2x}=-\frac{n_{up}}{\cos\theta}E_{2}$, $B_{\left(X\right)x}=n_{X}\left(\cos\theta_{X}-\sin\theta_{X}K_{X}\right)E_{X}$.

At the plasmon frequency, when $K_{X}=\frac{1}{\tan\theta_{X}}$, the last equation gives $B_{\left(X\right)x}=0$, as should be expected. For an isotropic medium, when $K_{X}=-\tan\theta_{X}$, we obtain $B_{\left(X\right)x}=\frac{n_{X}}{\cos\theta_{X}}E_{X}$ which is also expected for a transverse wave (note that $E_{X}$ is an amplitude of the y-component of the extraordinary (X-)mode).

We will use the boundary conditions

[TABLE]

[TABLE]

to obtain

[TABLE]

where $n_{up}\sin\theta=n_{X}\sin\theta_{X}$. In the limit of an isotropic medium, where $K_{X}=-\tan\theta_{X}$, $\sigma_{ij}^{S}=0$, we obtain $R=\frac{E_{2}}{E_{1}}=\frac{n_{up}\cos\theta_{X}-n_{X}\cos\theta}{n_{X}\cos\theta+n_{up}\cos\theta_{X}}$ which is a standard Fresnel equation.

For the normal incidence the expressions are simplified:

[TABLE]

which gives

[TABLE]

The contribution of surface states is less trivial for X-mode excitation as compared to the excitation of an $O$-mode. For normal incidence (see Eq. (69)) one can see that at the plasmon resonance frequency, when $\varepsilon_{zz}\rightarrow 0$ in the absence of losses, the contribution of the surface conductivity can become dominant. Indeed, in Eq. (69) the term $\sigma_{yz}^{S}\frac{g}{\varepsilon_{zz}}$ diverges as $\frac{1}{\varepsilon_{zz}}$, whereas the refractive index $n_{X}$ diverges weaker, as $\frac{1}{\sqrt{\varepsilon_{zz}}}$. When $\sigma_{ij}^{S}=0$ while $n_{X}\gg n_{up}$ we have $R=-1$ (we take into account that the magnitude of $n_{X}$ is large at the plasmon frequency). In the opposite case, when the contribution of the surface conductivity dominates, i.e. $\frac{4\pi}{c}|\sigma_{yz}^{S}\frac{g}{\varepsilon_{zz}}|\gg|n_{X}|\approx\frac{g}{\sqrt{|\varepsilon_{zz}|}}$ , we obtain $R=+1$ , i.e. the phase of the reflected field is rotated by 180 degrees.

The enhanced contribution of the surface conductivity at normal incidence in the vicinity of the bulk plasmon resonance is expected. Indeed, at plasmon resonance the $z$-component $E_{z}$ of the field in the medium becomes very large, which leads to a dominant contribution of the surface current $j_{y}^{S}=\sigma_{yz}^{S}E_{z}$.

For oblique incidence $\theta\neq 0$ and small losses the calculations of the reflection in the vicinity of plasmon resonance have a technical subtlety, related to the presence of the term $n_{X}\cos\theta\left(\cos\theta_{X}-\sin\theta_{X}K_{X}\right)$ in Eq. (68). Indeed, at the plasmon frequency $n_{X}\rightarrow\infty$ as losses $\gamma\rightarrow 0$; however, for a plasmon we also have $K_{X}\rightarrow\frac{1}{\tan\theta_{X}}$, i.e. $\left(\cos\theta_{X}-\sin\theta_{X}K_{X}\right)\rightarrow 0$. One needs to treat the resulting uncertainty of the product with caution. The details are presented in Appendix F.

The main result is that the contribution of surface states to the reflected wave is determined by the ratio

[TABLE]

and therefore becomes significant or dominant at the plasmon resonance frequency, when $\varepsilon_{zz}=\varepsilon_{zz}^{\left(0\right)}+i\frac{4\pi}{\omega}\sigma_{zz}^{B}\rightarrow 0$. When the bulk contribution dominates the reflection coefficient $R$ is close to $-1$. When the surface contribution dominates, $R$ is close to $+1$ i.e. the phase of the reflected field flips.

VII Surface plasmon-polaritons

Surface plasmon-polaritons can be supported by both bulk and surface electron states. Here we derive dispersion relations for surface waves including both bulk and surface conductivity for several specific cases. Emphasis is placed on the situations where the dispersion is significantly affected or dominated by surface states and can therefore be used for diagnostics of surface states and Fermi arcs. Previously, surface plasmons in WSMs have been considered in the low-frequency limit within a semiclassical description of particle motion with added ad hoc anomalous Hall term song2017 and with a quantum-mechanical description andolina2018 based on the Hamiltonian in okugawa2014 . Both studies indicated strong anisotropy and dispersion of surface plasmons.

VII.1 Quasielectrostatic approximation

Within the quasielectrostatic approximation the electric field can be defined through the scalar potential:

[TABLE]

We introduce the vector of electric induction, $\mathcal{\vec{D}}=\real\left[\vec{D}\left(z\right)e^{ik_{x}x+ik_{y}y-i\omega t}\right]=\hat{\varepsilon}\mathcal{\vec{E}}$ and use Gauss’ law for each halfspace:

[TABLE]

In general, there can be an electric dipole layer at the boundary between the two media. The dipole layer has a jump in the scalar potential $\Phi\left(z\right)$,

[TABLE]

where $d_{z}$ is determined by Eqs. (62).

Next, we define the potential $\Phi\left(z\right)$ for the surface mode as

[TABLE]

Using Eq. (70) in each halfspace, we obtain

[TABLE]

Using the boundary condition Eq. (58) we get

[TABLE]

which gives

[TABLE]

where $\varepsilon_{yz}=-\varepsilon_{zy}=ig=i\frac{4\pi\sigma_{yz}^{B}}{\omega}$. Using also the boundary condition Eq. (71) together with Eqs. (62), we obtain

[TABLE]

From these relationships one can get the dispersion equation for surface waves. Note that the confinement constants $\kappa_{W}$ and $\kappa_{up}$ are generally complex-valued. Their imaginary parts give rise to a Poynting flux away from the surface which contributes to surface wave attenuation.

VII.1.1 Neglecting surface states

First, we neglect the surface conductivity to consider surface plasmons supported by bulk carriers only. In this case from Eqs. (72), (75) we get $\kappa_{up}=\sqrt{k_{x}^{2}+k_{y}^{2}}$, $\Phi_{up}=\Phi_{W}$. Denoting $k_{x}^{2}+k_{y}^{2}=k^{2}$, $k_{x}=k\cos\phi$, $k_{y}=k\sin\phi$, we obtain from Eq. (73)

[TABLE]

Furthermore, from Eq. (74) for $\kappa_{up}=k$ and $\Phi_{up}=\Phi_{W}$ we have

[TABLE]

where $\varepsilon_{yz}=ig=i\frac{4\pi\sigma_{yz}^{B}}{\omega}$. Substituting Eq. (76) into Eq. (77), we obtain the dispersion relation

[TABLE]

The dispersion equation Eq. (78) gives the dependence $\omega\left(\phi\right)$, but does not have any dependence on the magnitude of $k$. This situation is similar to the dispersion relation for bulk plasmons in the quasielectrostatic approximation, Eq. (56). It is also similar to waves in classical magnetized plasmas. Of course the range of values of $k$ is constrained by the validity of the quasielectrostatic approximation.

VII.1.2 Including surface states

If we now include the surface conductivity, Eqs. (72)-(75) give

[TABLE]

where the function $D\left(\omega,\phi\right)$ is determined by Eq. (78). As we see, taking the surface conductivity into account brings the dependence on the magnitude of the wave vector $k$ into the dispersion relation. Therefore, measuring the frequency dispersion of the surface plasmon resonance provides a direct characterization of surface states.

Figure 15 shows the surface plasmon dispersion for propagation along $y$, i.e. transverse to the gyrotropy $x$-axis, for two values of the Fermi momentum. The real part of the surface plasmon frequency ignoring the contribution of the surface conductivity is shown as a dashed horizontal line for each value of $k_{F}$. Clearly, the contribution of surface electron states is important everywhere, except maybe in a narrow region of small wavenumbers $k$ where the quasistatic approximation breaks down. The plot has a horizontal axis $ck$ in units of meV in order to directly compare with frequencies. The inequality $ck\gg\omega$ is satisfied almost everywhere.

The fact that the contribution of the surface current is so important can be understood from the structure of Eq. (79). Clearly, the relative contribution of the bulk and surface terms can be estimated by comparing the magnitudes of $|\sigma^{B}|$ and $|k\sigma^{S}|$ where $\sigma^{B}$ and $\sigma^{S}$ are appropriate components of bulk and surface conductivity tensors and $k$ is a wavenumber of a given electromagnetic mode. This is true not only for surface modes but also for other electromagnetic wave processes at the boundary such as reflection. In the mid/far-infrared spectral region of interest to us, $|k\sigma^{S}|\ll|\sigma^{B}|$ for vacuum wavelengths $ck\sim\omega$. However, for large surface plasmon wavenumbers shown in Fig. 15 the opposite condition $|k\sigma^{S}|\geq|\sigma^{B}|$ is satisfied.

Note the dispersion in Fig. 15 is stronger (the slope is steeper) at frequencies corresponding to Re$[\epsilon_{zz}]\approx 0$, i.e. near the resonance for bulk plasmons propagating along $z$. This follows from Eq. (79) where the surface terms contain a factor $1/\sqrt{\epsilon_{zz}}$. Physically, this is expected: indeed, as we already commented, at the plasmon resonance the $z$-component $E_{z}$ of the field in the medium becomes very large, which leads to an enhanced contribution of the surface current $j_{y}^{S}=\sigma_{yz}^{S}E_{z}$.

VII.2 Surface waves beyond the quasielectrostatic approximation

For small wavenumbers the quasielectrostatic approximation is no longer valid. On the other hand, in this case one can neglect the surface conductivity as we pointed out in the previous paragraph. This is not an interesting limit as far as the spectroscopy of surface states is concerned, but we still derive the resulting dispersion relation for completeness. For the electric field of a surface mode in the upper halfspace with the refractive index $n_{up}$,

[TABLE]

the Maxwell’s equation for $\nabla\times\vec{\cal{E}}$ gives

[TABLE]

For the field in the Weyl semimetal,

[TABLE]

the same equation gives, after replacing $\kappa_{up}\rightarrow-\kappa_{W}$ in Eq. (80),

[TABLE]

The inverse decay length for the field in the upper halfspace is given by $\kappa_{up}^{2}=k^{2}-n_{up}^{2}\frac{\omega^{2}}{c^{2}}$.

In a WSM we can use a version of Eq. (41) after replacing $k_{z}\rightarrow-i\kappa_{W}$:

[TABLE]

where $k^{2}=k_{x}^{2}+k_{y}^{2}$.

Consider again a surface wave propagating transverse to the anisotropy axis ($k_{x}=0$). In this case, there are two solutions to the dispersion equation Eq. (82), an O-wave and an X-wave. However, one can show that an O-wave with $E_{x}\neq 0$ does not exist as a surface wave. Moreover, this statement remains true even with the surface current taken into account. Only the X-wave with $E_{y,z}\neq 0$ can exist as a surface wave. Its inverse confinement length in the Weyl semimetal is given by

[TABLE]

where

[TABLE]

is the refractive index of an extraordinary wave propagating in the volume in the $y$-direction (see Eq. (43) for $\theta=\frac{\pi}{2}$). The polarization of an extraordinary wave is determined by

[TABLE]

which follows from Eq. (82). After some straightforward algebra, we obtain the dispersion relation for a surface wave:

[TABLE]

In the limit of large wavenumbers $k$ this equation becomes the quasielectrostatic dispersion relation Eq. (78) at $\phi=\frac{\pi}{2}$.

For the propagation in $x$-direction, one can repeat the above analysis for the case $k_{y}=0$ and obtain that there are no surface wave solutions when the surface conductivity is neglected.

One interesting solution of the dispersion equation Eq. (85) is a strongly nonelectrostatic case when the surface mode is weakly localized in a medium above the WSM surface, e.g. in the air. The energy of this wave is mostly contained in an ambient medium above the WSM surface where there is no absorption. Therefore, such surface waves can have a long propagation length; see e.g. barnes ; homola ; agran .

To find this solution we assume $n_{up}^{2}=1$ and introduce the notation $\frac{\omega}{c}=k_{0}$. A weak localization outside a WSM means that $|\kappa_{up}|\ll k_{0}$. Then, assuming $k\simeq k_{0}+\delta k$, where $k_{0}\gg|\delta k|$, we obtain $\kappa_{up}\simeq\sqrt{2k_{0}\delta k}$. From Eqs. (85) and (83) in the first order with respect to $\sqrt{\frac{\delta k}{k_{0}}}$ we get

[TABLE]

[TABLE]

This solution describes surface waves if Re$[\kappa_{W}]>0$ and Re$[\kappa_{up}]>0$. In addition, $|\delta k|\ll k_{0}$ has to be satisfied. We checked that all three inequalities are satisfied for the numerical parameters chosen to calculate the conductivity tensor. As an example, Fig. 16 shows normalized confinement constants Re$[\kappa_{W}]/k_{0}$ and Re$[\kappa_{up}]/k_{0}\simeq{\rm Re}[\sqrt{2\delta k/k_{0}}]$ as functions of frequency, for the Fermi momentum $k_{F}=0.5b$. Clearly, the solution describes a surface wave which is weakly confined in the air and strongly confined in the WSM. The spectra remain qualitatively the same with increasing Fermi momentum, but the oscillating feature moves to higher energies, roughly following the spectral region where the real parts of $\varepsilon_{zz}$ and $\varepsilon_{yy}$ cross zero. We note again that the confinement constants $\kappa_{W}$ and $\kappa_{up}$ are complex-valued. Their imaginary parts give rise to a Poynting flux away from the surface which contributes to surface wave attenuation.

VIII Summary and conclusions

We presented systematic studies of the optical properties and electromagnetic modes of Weyl semimetals. Both bulk and surface conductivity tensors are derived from a single microscopic Hamiltonian. The presence of separated Weyl nodes and associated surface states gives rise to distinct signatures in the transmission, reflection, and polarization of bulk and surface electromagnetic waves. These signatures can be used for quantitative characterization of electronic structure of Weyl semimetals. Particularly sensitive spectroscopic probes of bulk electronic properties include strong anisotropy in propagation of both bulk and surface modes, birefringent dispersion and absorption spectra of ordinary and extraordinary normal modes, the frequency of bulk plasmon resonance as a function of incidence angle and doping level, and the polarization rotation and ellipticity for incident linearly polarized light. The sensitive characterization of surface electronic states can be achieved by measuring the phase change of the reflection coefficient of incident plane waves, the frequency dispersion of surface plasmon-polariton modes, and strong anisotropy of surface plasmon-polaritons with respect to their propagation direction and polarization.

Potential optoelectronic applications of Weyl semimetal films in the mid-infrared and THz spectral regions will benefit from the strong anisotropy, gyrotropy, and birefringence of these materials, giant polarization rotation for light transmitted along the gyrotropy axis of submicron films, and strongly localized surface plasmon-polariton modes. All effects are tunable by doping.

Acknowledgements.

The authors are grateful to I. Tokman and M. Erukhimova for useful discussions. This work has been supported by the Air Force Office for Scientific Research through Grants No. FA9550-17-1-0341 and FA9550-14-1-0376. I.O. and M.T. acknowledge the support by the Ministry of Education Science of the Russian Federation contract No. 14.W03.31.0032.

Appendix A Evaluation of the matrix elements of the current density operator

We denote everywhere the bulk states by Latin letters, and the surface states by Greek letters, i.e. $\left|n\right\rangle=\left|B\right\rangle,\left|\mu\right\rangle=\left|S\right\rangle$. In this section we evaluate the matrix elements of the current density operator that enter Eq. (34) for the components of bulk and surface conductivity tensors.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where we have used $\kappa=\frac{b^{2}-\left(k_{x}^{2}+k_{y}^{2}\right)}{2b}$.

Appendix B Calculation of the bulk optical conductivity tensor

The 3D integrals over electron momenta cannot be evaluated analytically in most cases, even in the zero temperature limit. Whenever the integrals remain in the final expression, they were evaluated numerically for the plots in the main text.

B.1 Contribution of intraband transitions ($s=+1\rightarrow s=+1$)

In this case the matrix elements $\boldsymbol{j}_{nm}^{(\boldsymbol{q})}$ of the current density operator reduce to

[TABLE]

[TABLE]

[TABLE]

Therefore, we obtain

[TABLE]

Similarly,

[TABLE]

[TABLE]

Here $\Theta(k)$ is the step function and we have used $\cos\theta_{\boldsymbol{k}}=\frac{k_{y}}{\sqrt{K_{x}^{2}+k_{y}^{2}+k_{z}^{2}}}$, $e^{i\phi_{\boldsymbol{k}}}=\frac{K_{x}+ik_{z}}{\sqrt{K_{x}^{2}+k_{z}^{2}}}$ , $K_{x}\equiv\frac{\left(k_{x}^{2}+k_{y}^{2}\right)-b^{2}}{2b}$, and $k_{F}\equiv\frac{E_{F}}{\hbar v_{F}}$.

[TABLE]

B.2 Contribution of interband transitions ($s\rightarrow-s,\,\left|B\right\rangle\leftrightarrow\left|S\right\rangle$)

In this case, i.e. $s_{m}=-s_{n}=\pm 1,\,n\neq m,$ the matrix elements $\boldsymbol{j}_{nm}^{(\boldsymbol{q})}$ of the current density operator reduce to

[TABLE]

[TABLE]

[TABLE]

where $n\neq m.$ Therefore, we obtain

[TABLE]

where we have used $K_{x}\equiv\frac{\left(k_{x}^{2}+k_{y}^{2}\right)-b^{2}}{2b}=-\kappa,$ \cos\theta_{\boldsymbol{k}}\left(-k_{x}\right)=\cos\theta_{\boldsymbol{k}}\left(k_{x}\right),$$\sin\theta_{\boldsymbol{k}}\left(-k_{x}\right)=\sin\theta_{\boldsymbol{k}}\left(k_{x}\right) $\cos\phi_{\boldsymbol{k}}\left(-k_{x}\right)=\cos\phi_{\boldsymbol{k}}\left(k_{x}\right),$ and $\sin\phi_{\boldsymbol{k}}\left(-k_{x}\right)=\sin\phi_{\boldsymbol{k}}\left(k_{x}\right).$

Similarly,

[TABLE]

[TABLE]

The only nonzero off-diagonal element is $\sigma_{zy}^{inter}(\omega)=-\sigma_{yz}^{inter}(\omega)$, as expected:

[TABLE]

Here we have introduced a cutoff at $k=K$ in the integration over electron momenta in order to regularize the divergent integral $\int\frac{d^{3}k}{\left(2\pi\right)^{3}}$ which comes from $\frac{1}{V}\sum_{n}\rightarrow\int\frac{d^{3}k}{\left(2\pi\right)^{3}}$. The divergence is an artifact of the effective Hamiltonian Eq. (1) which has a “bottomless” valence band with electrons occupying all states to $k\rightarrow\infty$. The regularization makes the valence band bounded from below. We chose the cutoff at the momentum corresponding to the energy of 2 eV, i.e. much higher than the range of interest to us near the Weyl nodes. In the numerical examples in the paper the value of half-separation between Weyl nodes $\hbar v_{F}b$ is chosen to be 100 meV. We have verified that an exact value of the cutoff has a negligible effect on the low-energy optical response below 350 meV, as long as $K$ is large enough.

Appendix C Calculation of the surface electrical conductivity

C.1 Surface-to-surface states intraband transitions

[TABLE]

All other tensor components are equal to zero.

C.2 Surface-to-bulk states transitions

[TABLE]

Similarly,

[TABLE]

[TABLE]

The only nonzero off-diagonal element is

[TABLE]

In Eqs. (112)-(115) the integral over $k_{z}$ can be carried out analytically in terms of elementary functions, leading however to very lengthy expressions which we do not present here. The remaining integration was carried out numerically. All integrals are finite, i.e. no cutoff is necessary.

Appendix D Drude-like low-frequency limit

In the limit when the frequency and the Fermi energy are much smaller than $\hbar v_{F}b$, only the electron momenta close to the corresponding Weyl point $k_{x}=\pm b$ matter. Therefore, we introduce $\delta k_{x}=k_{x}-b$ for electron states near one Weyl point and replace the degeneracy factor by $2\times g$ to account for the contribution from the second Weyl point. In this case, K_{x}\sim\frac{\left(k_{x}-b\right)\left(k_{x}+b\right)}{2b}\approx\delta k_{x},\,k_{x}$$=b+\delta k_{x}, and all diagonal intraband components have the same Drude form:

[TABLE]

All off-diagonal conductivity elements are zero in this limit.

Appendix E Small $b$ Expansion

In the limit $b\ll 1$, we can expand the conductivity in powers of $b$ to the leading order: $b\ll 1$, $\frac{1}{b}\gg 1,$ $K_{x}=\frac{\left(k_{x}^{2}+k_{y}^{2}\right)-b^{2}}{2b}\sim\frac{\left(k_{x}^{2}+k_{y}^{2}\right)}{2b}\sim\frac{\left(k_{x}^{2}+k_{y}^{2}+k_{z}^{2}\right)}{2b}\gg k_{x,y,z},\frac{\omega}{v_{F}}$ for $k_{x,y,z}\neq 0$. Then we obtain

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

All off-diagonal surface terms are zero.

Appendix F Reflection in the vicinity of plasmon resonance

For oblique incidence $\theta\neq 0$ and small losses the calculations of the reflection in the vicinity of plasmon resonance have a technical subtlety, related to the presence of the term $n_{X}\cos\theta\left(\cos\theta_{X}-\sin\theta_{X}K_{X}\right)$ in Eq. (68). Indeed, at the plasmon frequency $n_{X}\rightarrow\infty$ as losses $\gamma\rightarrow 0$; however, for a plasmon we also have $K_{X}\rightarrow\frac{1}{\tan\theta_{X}}$, i.e. $\left(\cos\theta_{X}-\sin\theta_{X}K_{X}\right)\rightarrow 0$. One needs to treat the resulting uncertainty of the product with caution.

We substitute the relationship $\sin\theta_{X}=\frac{n_{up}\sin\theta}{n_{X}}$ into the expression for the refractive index of an extraordinary wave:

[TABLE]

which gives

[TABLE]

In the case $\varepsilon_{yy}=\varepsilon_{zz}=\varepsilon_{\perp}$, Eq. (122) for an arbitrary angle $\theta$ leads to the familiar expression $n_{X}^{2}=\varepsilon_{\perp}-\frac{g^{2}}{\varepsilon_{\perp}}$. Next we use Eq. (46):

[TABLE]

Consider the expression $n_{X}\cos\theta\left(\cos\theta_{X}-\sin\theta_{X}K_{X}\right)$:

[TABLE]

The condition $\frac{n_{X}}{n_{up}}\gg 1$, which is satisfied at the plasmon frequency, allows one to simplify the above expressions for any angle of incidence $\theta$

[TABLE]

[TABLE]

Since for $\frac{n_{X}}{n_{up}}\gg 1$ we always have $\sin\theta_{X}\ll 1$, the plasmon frequency always corresponds to $|\varepsilon_{zz}|\ll 1$ (at normal incidence, $\varepsilon_{zz}=0$ exactly). Taking into account Eq. (122), we obtain $1\gg|\varepsilon_{zz}|\sim n_{X}^{-2}$.

Now let us consider the range of incidence angles close to normal incidence, when $\sin^{2}\theta\ll 1$. Two cases need to be treated separately: $|\varepsilon_{zz}|\ll\sin^{2}\theta n_{up}^{2}\ll 1$ and $\sin^{2}\theta n_{up}^{2}\ll|\varepsilon_{zz}|\ll 1$.

(i) $|\varepsilon_{zz}|\ll sin^{2}\theta n_{up}^{2}\ll 1$

In this case

[TABLE]

[TABLE]

where $g=\frac{4\pi\sigma_{yz}^{B}}{\omega}$,

[TABLE]

For real $\sigma_{yz}^{\left(B,S\right)}$ we always have $|R|=1$; however, the phase of the reflected field depends on the contribution of surface states. Since in the vicinity of plasmon resonance $n_{X}\sim\frac{1}{\sqrt{|\varepsilon_{zz}|}}\gg 1$, at these frequencies the contribution of surface states may become important. This is especially clear in the limit of small enough angles, when $n_{up}^{2}sin\theta\ll|\frac{4\pi\sigma_{yz}^{B}}{\omega}|$. In this case

[TABLE]

When the bulk contribution dominates we have $R=-1$, whereas if the surface contribution dominates we obtain $R=+1$, i.e. the phase of the reflected field flips.

The relative contribution of surface states is determined by the ratio $\frac{|\sigma_{yz}^{S}n_{X}|}{\frac{c}{\omega}|\sigma_{yz}^{B}|}$. Taking into account that $|n_{X}|\approx\frac{|g|}{\sqrt{|\varepsilon_{zz}|}}$ and $|g|=\frac{4\pi|\sigma_{yz}^{B}|}{\omega}$, the above ratio can be reduced to $\frac{\frac{4\pi|\sigma_{yz}^{S}|}{c}}{\sqrt{|\varepsilon_{zz}|}}$.

(ii) $\sin^{2}\theta n_{up}^{2}\ll|\varepsilon_{zz}|\ll 1$

This case is similar to the one at $\theta=0$. Indeed, for this range of parameters we obtain

[TABLE]

[TABLE]

[TABLE]

Eqs. (129), (130) are the same as for the normal incidence. Eq. (131) can be obtained from the normal incidence formula Eq. (69) if $|\sigma_{yy}^{S}|\ll|\sigma_{yz}^{S}\frac{g}{\varepsilon_{zz}}|$ and $n_{X}\gg n_{up}$; the latter inequalities are valid near the plasmon resonance, where $n_{X}\sim\frac{1}{\sqrt[]{|\varepsilon_{zz}|}}\rightarrow\infty$.

For real values of $\sigma_{yz}^{\left(S\right)}$ we always have $|R|=1$, but the phase of the reflected field depends on the contribution of surface states. Again, when the bulk contribution dominates we have $R=-1$, whereas if the surface contribution dominates we obtain $R=+1$.

The relative contribution of surface states is determined by the ratio $\frac{\frac{4\pi}{c}|\sigma_{yz}^{S}\frac{g}{\varepsilon_{zz}}|}{|n_{X}|}$. Again taking into account $|n_{X}|\approx\frac{|g|}{\sqrt{|\varepsilon_{zz}|}}$ and $|g|=\frac{4\pi|\sigma_{yz}^{B}|}{\omega}$ we obtain that the above ratio is reduced to exactly the same expression as before: $\frac{4\pi|\sigma_{yz}^{S}|/c}{\sqrt{|\varepsilon_{zz}|}}$.

To summarize, the effect of surface states on the reflected wave is determined by the ratio

[TABLE]

and therefore becomes significant or dominant at the plasmon resonance frequency, when $\varepsilon_{zz}=\varepsilon_{zz}^{\left(0\right)}+i\frac{4\pi}{\omega}\sigma_{zz}^{B}\rightarrow 0$.