## 2 Theory of Isofrequency Surfaces in Bianisotropic Materials  
Consider bianisotropic materials with \(6\times 6\) effective parameters matrix \(\hat{M}=\left[ \begin{array}{cc}\hat{e}&\hat{X}\\ \hat{Y}&\hat{\mu}\end{array} \right]\), such that the constitutive relations are given by \(\left( \begin{array}{c}\mathbf{D}\\ \mathbf{B}\end{array} \right)=\hat{M}\left( \begin{array}{c}\mathbf{E}\\ \mathbf{H}\end{array} \right)\). It is known that the isofrequency surface of the most generic bianisotropic material with arbitrary material parameters \(\hat{M}\) is a quartic surface in k-space given by [10]. \[f(k,k_{0})=\sum_{i+j+l+m=4}[\alpha_{ijlm}k_{x}^{i}k_{y}^{j}k_{z}^{l}k_{0}^{m}] =0. \tag{1}\]  
As was demonstrated in Ref. [10], the coefficients \(\alpha_{ijlm}\) follow from the solution of [PERSON]'s equation. The topological asymptotic skeleton of the isofrequency surfaces (Equation (1)) can be found in the high-k limit \(k\gg k_{0}\). The high-k states are those tending to the conical k-surfaces given by  
\[h(\mathbf{k})=f(k\rightarrow\infty,k_{0})=\sum_{i+j+l=4}[\alpha_{ij0}k_{x}^{i}k_{y }^{j}k_{z}^{l}]=0 \tag{2}\]  
Below we establish the detailed form and properties of the function \(h(\mathbf{k})\) in terms of the material parameters \(\hat{M}\) and the topological properties of the skeleton asymptotic surfaces of Equation (2). Figure 1: Isofrequency k-surface in k-space with high-k polaritons (HKPs) in a conventional monohyperbolic material (**a**) and a tetra-hyperbolic material (**b**) introduced in this paper and described in detail in the discussion below. We start by considering [PERSON]'s equations for the amplitudes \(\hat{\Gamma}=\left( \begin{array}{c}E\\ H\end{array} \right)\) of the plane waves with wave vectors \(\mathbf{k}\) and frequencies \(\omega=ck_{0}\) written as  
\[k\hat{Q}\hat{\Gamma}=k_{0}\hat{M}\hat{\Gamma},\hat{Q}=\left[ \begin{array}{cc} \hat{0}&\hat{R}\\ -\hat{R}&\hat{0}\end{array} \right]\hat{R}=\left( \begin{array}{ccc}0&-k_{z}& k_{y}\\ k_{z}&0&-k_{x}\\ -k_{y}&k_{x}&0\end{array} \right) \tag{3}\]  
Consider a wave propagating in the direction \(\hat{\mathbf{k}}=\mathbf{k}/k=\left(q_{x},q_{y},q_{z}\right)\). Let us use the transformation \(\hat{T}\):  
\[\hat{T}=\left( \begin{array}{ccc}q_{x}q_{x}/q_{p}&q_{z}q_{y}/q_{\rho}&-q_{\rho }\\ -q_{y}/q_{\rho}&q_{x}/q_{\rho}&0\\ q_{x}&q_{y}&q_{z}\end{array} \right)=\left( \begin{array}{c}\mathbf{q}_{1}\\ \mathbf{q}_{2}\\ \mathbf{q}_{3}=\hat{\mathbf{k}}\end{array} \right)\]  
where \(q_{\rho}=\sqrt{q_{x}^{2}+q_{y}^{2}}\) and vectors \(q_{i}\) are the rows of matrix \(\hat{T}\). Note that \(\hat{T}=\left(\hat{T}^{T}\right)^{-1}\) and is a proper rotation \(\det\hat{T}=1\). Let us apply transformation \(\hat{T}\) to Maxwell's Equation (Equation (3))  
\[k\overline{Q\Gamma}=k_{0}\overline{M\Gamma},\overline{\Gamma}=\hat{T}_{6}\hat {\Gamma},\overline{Q}=\hat{T}_{6}\hat{Q}\hat{\Gamma}_{6}^{-1}=\left[ \begin{array} []{cc}\hat{0}&\overline{R}\\ -\overline{R}&\hat{0}\end{array} \right] \tag{4}\]  
\[\overline{R}=\left( \begin{array}{ccc}0&-1&0\\ 1&0&0\\ 0&0&0\end{array} \right)\overline{M}=\hat{T}_{6}\hat{M}\hat{\Gamma}_{6}^{-1}= \left[ \begin{array}{cc}\overline{\epsilon}&\overline{X}\\ \overline{Y}&\overline{\mu}\end{array} \right]\hat{T}_{6}=\left[ \begin{array}[ ]{cc}\hat{T}&0\\ 0&\hat{T}\end{array} \right] \tag{5}\]  
According to Equation (5) the transformed matrices \(\widetilde{a}=\overline{\epsilon},\overline{\mu},\overline{X},\overline{Y}\) can be expressed in terms of the original \(\hat{a}=\hat{e},\hat{\mu},\hat{X},\hat{Y}\) as  
\[\widetilde{a}_{ij}=q_{i}\cdot\hat{a}\cdot\mathbf{q}_{j}^{T} \tag{6}\]  
The zeroes in the bottom row of matrix \(\overline{R}\) are due to the conditions \(\mathbf{k}\cdot\mathbf{D}=0\) and \(\mathbf{k}\cdot\mathbf{B}=0\). Thus, the longitudinal components of fields \(E_{z}\) and \(H_{z}\) can be decoupled from the transverse field components in the transformed coordinates as  
\[\left( \begin{array}{c}E_{z}\\ H_{z}\end{array} \right)=-\left(