A group theoretical route to deterministic Weyl points in chiral photonic lattices

TL;DR

This paper demonstrates the existence of symmetry-induced, deterministic Weyl points with non-trivial topology in chiral cubic photonic systems, using group theory and physical realizations like nano plasmonic and photonic crystal structures.

Contribution

It introduces a group theoretical approach to identify deterministic Weyl points in chiral photonic lattices, moving beyond accidental degeneracies.

Findings

Existence of symmetry-induced Weyl points at Brillouin zone center

Realization of these points in nano plasmonic and photonic crystal systems

Topologically non-trivial dispersion with hyper-conic features

Abstract

Classical topological phases derived from point degeneracies in photonic bandstructures show intriguing and unique behaviour. Previously identified exceptional points are based on accidental degeneracies and subject to engineering on a case-by-case basis. Here we show that symmetry induced (deterministic) pseudo Weyl points with non-trivial topology and hyper-conic dispersion exist at the centre of the Brillouin zone of chiral cubic systems. We establish the physical implications by means of a $P2_13$ sphere packing, realised as a nano plasmonic system and a photonic crystal.

Tables

Table 1. Table 1: Character table for the 23 23 23 tetrahedral point group, with ω := exp ⁡ { 2 ​ ı ​ π / 3 } assign 𝜔 2 italic-ı 𝜋 3 \omega:=\exp\left\{2\imath\pi/3\right\} . The last column lists the Frobenius Schur type of the representation.

$23$	$\mathrm{𝟙}$	$3C_2$	$4C_3$	$4C_3^{\prime }$	FS
$A$	1	1	1	1	(a)
$E_{+}$	1	1	$\omega$	${\omega }^2$	(b)
$E_{-}$	1	1	${\omega }^2$	$\omega$	(b)
$T$	3	-1	0	0	(a)

Equations

$$W_{\alpha \beta }(\mathbf{k})=ıd\sum _{\gamma }{ϵ}_{\alpha \beta \gamma }{k}_{\gamma }\text{ ,}$$

$${\mathbf{p}}_i=\alpha (k_0)\sum _{j\ne i}\mathcal{G}({\mathbf{r}}_i-{\mathbf{r}}_j,k_0){\mathbf{p}}_j\text{ .}$$

$${\alpha }^{-1}(k_0){\mathbf{p}}_{\mu }=\sum _{\nu }{M}_{\mu \nu }(\mathbf{k},k_0){\mathbf{p}}_{\nu }\text{ .}$$

$${\Delta }_{\mathbf{k}=0,i}(\left\{p_{\alpha }|t_{\alpha }+T\right\}):=D_i(p_{\alpha })\text{ .}$$

$$A+E_{+}+E_{-}+3T$$

$$n_i=\frac{1}{N}\sum _R\overline{\chi }(R){\chi }_i(R)\text{ ,}$$

$$M(\lambda )\mathbf{v}=0\text{ ,}$$

$$\mathbf{v}={\mathbf{v}}_{\alpha }^{(i)}(\mathbf{k},{\mathbf{k}}_n)\text{ .}$$

$$\exp \left\{-ı{\mathbf{k}}_n\cdot \mathbf{r}\right\}{\mathbf{v}}_{\alpha }^{(i)}(\mathbf{k},{\mathbf{k}}_n)\equiv |\mathbf{k}i,n\alpha ⟩\text{ .}$$

$$M({\mathbf{k}}_n,\lambda )|\mathbf{k}i,n\alpha ⟩=0\text{ .}$$

$$\displaystyle\Big{(}M(\bm{k}_{n}^{(0)},\lambda_{0})$$

$$+\delta \mathbf{k}\cdot {\left[{\nabla }_{\mathbf{k}}M(\mathbf{k},{\lambda }_0)\right]}_{\mathbf{k}={\mathbf{k}}_n^{(0)}}$$

$$\displaystyle+\delta\lambda\left[\partial_{\lambda}M(\bm{k}_{n}^{(0)},\lambda)\right]_{\lambda=\lambda_{0}}\Big{)}$$

$$\sum _{\alpha }{c}_{\alpha }|\mathbf{k}i,n\alpha ⟩$$

$$\sum _{\beta }{\mathbf{W}}_{\alpha \beta }\cdot \delta \mathbf{k}{c}_{\beta }=-\sum _{\beta }\delta \lambda E_{\alpha \beta }{c}_{\beta }\text{ ,}$$

$${\mathbf{W}}_{\alpha \beta }:=⟨\mathbf{k}i,n\alpha |{\left[{\nabla }_{\mathbf{k}}M(\mathbf{k},{\lambda }_0)\right]}_{\mathbf{k}={\mathbf{k}}_n^{(0)}}|\mathbf{k}i,n\beta ⟩$$

$$E_{\alpha \beta }:=⟨\mathbf{k}i,n\alpha |{\left[{\partial }_{\lambda }M({\mathbf{k}}_n^{(0)},\lambda )\right]}_{\lambda ={\lambda }_0}|\mathbf{k}i,n\beta ⟩\text{ .}$$

$$P{E}_{\alpha \beta }$$

$$P{W}_{\alpha \beta }^{(n)}$$

$${\Delta }_{(\alpha \beta ),(\gamma \delta )}(P):=$$

$${\Lambda }_{(\alpha \beta n),(\gamma \delta m)}(P):=$$

$$-{\delta }_{\alpha \gamma }{\delta }_{\beta \delta }{\delta }_{nm}\text{ .}$$

$$\frac{1}{{\sum }_P}\sum _P\overline{{\chi }_D}(P){\chi }_D(P)$$

$$\frac{1}{{\sum }_P}\sum _P\overline{{\chi }_D}(P){\chi }_D(P)\text{Tr}\{R(P)\}$$

$$E=c\text{mathbbm}1$$

$$\mathbf{W}=d_1\left(\begin{array}{c}\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 1& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& 0& 1\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& 0& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right)\end{array}\right)+d_2\left(\begin{array}{c}\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 1& 0& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& 1& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)\end{array}\right)\text{ .}$$

$$ıd\sum _{\beta \gamma }{\epsilon }_{\alpha \beta \gamma }\delta k_{\gamma }{c}_{\beta }=\delta \lambda c_{\alpha }\text{ ,}$$

$$W_{\alpha \beta }(\mathbf{k})=ıd\sum _{\gamma }{ϵ}_{\alpha \beta \gamma }{k}_{\gamma }\text{ .}$$

$$ıdk\hat{\mathbf{k}}\times \mathbf{v}=\lambda \mathbf{v}$$

$$\left[\lambda ,\left(\begin{array}{c}{v}_k\\ v_{\theta }\\ v_{\varphi }\end{array}\right)\right]=\left\{\left[0,\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right)\right],\left[\pm dk,\frac{1}{\sqrt{2}}\left(\begin{array}{c}0\\ 1\\ \pm ı\end{array}\right)\right]\right\}\text{ .}$$

$$K=d{t}_1\wedge d{t}_2\left[({\partial }_2{\mathbf{e}}_1,{\partial }_1{\mathbf{e}}_2)-({\partial }_1{\mathbf{e}}_1,{\partial }_2{\mathbf{e}}_2)\right]$$

$$\mathcal{B}=ı\text{Tr}\{PdP\wedge dP\}\text{ ,}$$

$$\mathcal{B}=2d\theta \wedge d\varphi \text{Im}\left\{⟨{\partial }_{\varphi }v|{\partial }_{\theta }v⟩\right\}\text{ .}$$

$$C=\frac{1}{2\pi }{\oint }_{\mathcal{M}}\mathcal{B}$$

Full text

Current broad interest in topological phases, triggered by the discovery of the quantum Hall effect Klitzing et al. (1980) and its theoretical investigation Thouless et al. (1982); Hatsugai (1993); Kohmoto (1985), can mainly be attributed to the fact that topological features are, due to their discrete nature, insensitive to system perturbations, and can, for example, give rise to the existence of topologically induced edge states for bulk systems Qi and Zhang (2011); Hasan and Kane (2010). Plasmonic Hofmann and Das Sarma (2016) and single electron Siroki et al. (2016) surface states of Weyl semi-metals, with an isolated exceptional point of non-trivial topology, are stable against perturbations and give rise to peculiar dynamics. Recently, it has been demonstrated that topological quantization occurs in entirely classical systems such as two-dimensional (2D) photonic crystals Onoda et al. (2004); Raghu and Haldane (2008), sparking a new wave of research on photonic topology Lu et al. (2014a). In particular, topologically protected Weyl points with hyperconic dispersion have been identified in double gyroid photonic crystals with broken parity-time symmetry Lu et al. (2013). Concurrently, group theory provides a tool to predict whether a given spatio-temporal symmetry permits topologically non-trivial exceptional points, or induces them deterministically. This idea has successfully found its way and been applied to classical Lu et al. (2014b) and quantum mechanical Wieder et al. (2016); Bradlyn et al. (2016) systems. Indeed, the existence of deterministic two and three-fold degeneracies at the center of the Brillouin zone (BZ), aka the $\Gamma$ point, for cubic symmetries is well known and documented in the literature Bradley and Cracknell (2010). Recently, it has been shown that some of these degeneracies are topologically non-trivial in electronic systems Kane (1957); Bradlyn et al. (2016).

Here we show on the basis of group and perturbation theory that symmetry induced three-fold degenerate pseudo Weyl points (PWPs) exist at the $\Gamma$ point in classical (photonic) systems. They split isotropically in first order in the Bloch wave vector $\bm{k}$ for any chiral cubic space group with time reversal symmetry. The PWPs studied here constitute a deterministic 3D analog to previously studied accidental Dirac points Huang et al. (2011). We show and demonstrate that they are of non-trivial topology, leading to protected surface states. In this letter, we first derive a 3D perturbation model that leads to hyperconic dispersion with non-trivial topology, and an intermediate flat band. We then construct a minimalistic geometry, a $P2_{1}3$ sphere packing (Fig. 1), which satisfies the symmetry requirements, and apply it to a quasistatic coupled-dipole model, before discussing topologically protected surface states that emanate from a PWP in a photonic crystal analog. This underscores that the existence of PWPs, including the peculiar transport properties of associated bulk and surface states, only depends on the underlying symmetry irrespective of the particular physical realization.

The theory applies to all linear and self-consistent physical systems with time reversal invariance and chiral cubic symmetry, with dynamics described by a Fourier integral over Hilbert states $|v(\omega)\rangle\in\mathcal{H}$ which individually solve a homogeneous (generally non-linear) frequency domain eigenproblem $M(\omega)\,|v(\omega)\rangle=0$. Symmetry requires that the $\mathcal{H}$ operator $M$ commutes with all elements of the underlying space group $\mathcal{G}$, represented by $\mathcal{H}$ operators $g$. As a consequence, a set of $N$ degenerate eigenvectors $|v_{n}(\omega)\rangle$ form a representation of $\mathcal{G}$, i.e. they span an $N$-dimensional vector space with an associated algebra that is homomorphic to $\mathcal{G}$. Irreducible representations impose a lower limit on the dimensionality of the respective vector space, resulting in deterministic degeneracies Bradley and Cracknell (2010); Dresselhaus (2013); Saba et al. (2013).

First order degenerate perturbation theory and representation theory, the latter of which provides the selection rules for the matrix elements within the former, allows to derive the slopes of the bandstructure at determinstic points of degeneracy (see supplementary material for details). For deterministic three-fold degeneracies at the $\Gamma$ point, this procedure yields a perturbation matrix $W(\bm{k})$ which is valid for small $k\ll 2\pi/a$ (with lattice constant $a$) and time reversal invariance

[TABLE]

with a free parameter $d\in\mathbb{R}$ and $\alpha,\beta,\gamma$ iterating over the three partners of the irreducible representation with $\bm{k}=0$ that span the degenerate eigenspace. Note the similarity of $W_{\alpha\beta}$ to the Weyl Hamiltonian $\mathcal{W}_{ij}$: the Pauli matrices $\sigma_{ij}^{(\gamma)}$ ($i,j\in\left\{1,2\right\}$) that occur in the latter are here replaced by the 3D Levi-Civita tensor $\epsilon_{\alpha\beta\gamma}$. The first order perturbation eigenvalues corresponding to $W(\bm{k})$ are $k_{0}^{(1)}:=\omega^{(1)}/c=\{0,\pm dk\}$: they only depend on the absolute value of $\bm{k}$ and describe isotropic hyperconic dispersion.

In the following, we shall define a PWP as the exceptional point $(0,k_{0}^{(0)})$ at which the two Weyl hypercones $(\bm{k},k_{0}^{(0)}\pm dk)$ in the four dimensional $(\bm{k},k_{0})$ parameter space meet. Although the bandstructure does not support a frequency with vanishing density of states due to the flat band, this definition is justified from a topological perspective: the correlated Chern numbers can be analytically calculated when integrating the Berry curvature over a small sphere in $\bm{k}$ space for each of the three bands. They evaluate to $C=0$ for the flat band and $C=\pm 2$ for the two hyperconic bands, showing a non-trivial topological signature, similar to a genuine Weyl point with Chern numbers $C=\pm 1$.

Analysing all 3D space groups Hahn (2002), it is straightforward to show that deterministic PWPs at the center of the BZ require chiral cubic symmetry. Interestingly, the trigonal groups $P312$ (149) and $P321$ (150) have two-dimensional representations which split into an anisotropic hypercone if time inversion is present, albeit not at the $\Gamma$ point Chen et al. (2015). A closely related matter is the non-existence of deterministic Dirac points at the $\Gamma$ point of two-dimensional crystals, including the famous honeycomb lattice Lu et al. (2014b).

To demonstrate the predicted behaviour we construct a chiral cubic sphere packing which is minimalistic in the sense that it generates the lowest dimensional vector space possible in models based on for example tight binding or pair interaction. A periodic sphere packing can be constructed by placing spheres on the Wyckoff point of multiplicity $N$ within a given space group Hahn (2002). Following this procedure for any non-symmorphic chiral cubic space group $\mathcal{G}$ and Wyckoff multiplicities smaller than $12$ yields a sphere packing that has the symmetry of an achiral supergroup $\mathcal{G}_{S}$. This counter-intuitive behaviour is related to the isotropy of the sphere as seed object that allows to introduce additional spurious (irregular) rotations. These unwanted rotations can be suppressed by the finite translation part in non-symmorphic symmetries. A seed sphere on the $4a$ Wyckoff point $(x,x,x)$ of $\mathcal{G}=P2_{1}3$ thus induces a chiral cubic sphere packing with $\mathcal{G}_{S}=\mathcal{G}$ and only $4$ spheres in the unit cell. Only if $x=n/4$ ($n\in\mathbb{Z}$), the sphere packing acquires the symmetry of the achiral supergroup $Pa\bar{3}$ (205); the chiral supergroup $P4_{3}32$ (212) is induced for $x=1/8+n/2$, and $P4_{1}32$ (213) for $x=-1/8+n/2$. Note that the introduction of an octahedral isogonal point symmetry instead of the tetrahedral symmetry of $P2_{1}3$ in these cases does not impose a change of bandstructure behaviour close to $\Gamma$. Fig. 1 illustrates the sphere packing for $x=3/8\,a$ ($P4_{1}32$ symmetry).

To elucidate the physics (ahead of a concrete experimental realization) let us consider an effective plasmonic model consisting of metallic nano-spheres of radius $\rho$ in vacuum (as in Han et al. (2009); Weick et al. (2013)). The position $\bm{r}_{i}$ of sphere $i$ shall be such that the distance $d_{ij}=|\bm{r}_{i}-\bm{r}_{j}|\gg\rho$ for any pair of spheres $(i,j)$. In the quasistatic approximation, Maxwell’s equations take the self-consistent form (acting on the dipole moments ${\bm{p}_{i}}$) Draine and Flatau (1994):

[TABLE]

Here, $\alpha(k_{0})=\rho^{3}(1-3k_{0}^{2}/k_{p}^{2})^{-1}$ is the polarizability of a metallic sphere in vacuum, that is modelled by a non-dissipative Drude response with plasma wave number $k_{p}$; $\mathcal{G}(\bm{r},k_{0})$ is the dyadic Green function for the monochromatic Maxwell wave operator.

If the spheres are arranged periodically as introduced above (cf. Fig. 1), the index $i$ is conveniently substituted by a multi-index $(\bm{n},\mu)\in\mathbb{Z}^{3}\times\{1,2,3,4\}$, with $\bm{r}_{\bm{n},\mu}=\bm{T}_{\bm{n}}+\bm{r}_{\mu}$ given by the sum of the lattice vector $\bm{T}_{\bm{n}}=a\,\bm{n}$ and the position within the unit cell $\bm{r}_{\mu}$. Bloch’s theorem then implies $\bm{p}_{\bm{n},\mu}=\bm{p}_{\mu}\exp\{\imath\bm{k}\cdot\bm{T}_{\bm{n}}\}$, so that Eq. 13 reduces to a family of $12$-dimensional non-linear Hermitian eigenproblems:

[TABLE]

Numerical challenges related to the convergence of the lattice sum $M_{\mu\nu}(\bm{k},k_{0})$ and their solution are solved in the supplementary material. Since the matrix $M$ generally imposes a small perturbation to the single sphere resonance solution $K_{0}^{(n)}:=\sqrt{3}\,k_{0}^{(n)}/k_{p}=1$ (due to $\rho^{-3}\gg 1$ in Eq. 19), the eigenvalue problem is linearized by approximating $M_{\mu\nu}(\bm{k},k_{0})\approx M_{\mu\nu}(\bm{k},k_{p}/\sqrt{3})=:M_{\mu\nu}(\bm{k})$. This assumption is inadequate close to the Ewald sphere $k_{0}=|\bm{k}|$, caused by poles in the diagonal entries of $M(\bm{k}=k_{0}\bm{\hat{k}})$, however, only affecting the two modes at the top and the bottom of the bandstructure on either side of the pole, cf. dashed red line in Fig. 2(a). The eigenvalues $\lambda_{n}(\bm{k})=\alpha^{-1}(k_{0})$ ($n=1,2,\dots,N$) of $M(\bm{k})$ can be obtained numerically with low computational cost. They produce the respective dispersion relation $K_{0}^{(n)}(\bm{k})=[1-\rho^{3}\lambda_{n}(\bm{k})]^{1/2}$, as shown in Fig. 2 for parameters $x/a=0.175$, $k_{0}a/(2\pi)=0.1$ and $\rho/a=0.1$. Fig. 2(b) shows an example within our model, where the first order perturbation outweighs higher orders even for relatively large Bloch wave number $k\approx\pi/(5a)$, so that an almost perfect hypercone can be observed. We find that, in contrast to this 3D representation, all first order perturbation matrix elements vanish for two-fold degeneracies at the $\Gamma$ point for any space group, cubic or non-cubic ($K_{0}=0.995$ in Fig. 2(a) constitutes an example). These exceptional points are henceforth not lifted in first order and are topologically trivial, with Chern number $C=0$ in both bands.

The universality of our results is vividly demonstrated if we replace the small metallic spheres by larger spheres of radius $\rho/a=0.25$ (fill fraction of $\pi/12\approx 26\%$), made of a high refractive index material with $n=4$, thus constructing a photonic crystal analog. The associated bandstructure (calculated with MPB Johnson and Joannopoulos (2001)) close to $k_{0}a/(2\pi)=0.5$ (supplementary figure $1$) resembles Fig. 2. A partial band gap opens in the projected bulk bandstructure with respect to a $[001]$ inclination in Fig. 3(a): this is the blue area of all $(\bm{k}_{\parallel},k_{0})$ for which at least one bulk mode exists for arbitrary $k_{z}\in\mathbb{R}$ Joannopoulos et al. (2008). Since the PWP degeneracy as well as the four-fold degeneracy at $R$ (projected onto $A$) is protected by cubic symmetry, this gap can be opened completely by e.g. perturbing the sphere positions (supplementary figure $2$). Topological surface states exist in the band gap at the interface between two enantiomorphic structures (with same bulk bands, but opposite chirality and Chern characteristics): Fig. 3(a) shows the surface mode dispersion of $12$ unit cells of a right handed crystal ($x/a=0.175$) and $12$ unit cells of a left handed crystal ($x/a=-0.175$) stacked in $[001]$ direction in a supercell geometry. The space group of the supercell is monoclinic with $P2_{1}/c$ (14) symmetry (note, however, that the Bravais lattice is tetragonal).

The supercell symmetry including time reversal requires all modes along $Z-A-X$ (Fig. 3(a), left inset) to be two-fold degenerate. However, the space group representations along $\Gamma-Z$ and $X-\Gamma$ (including $\Gamma$ itself) are 1D. The surface states still stick together. This can be understood as follows: consider a single interface surface mode along $\Gamma-Z$. Along this path, $\bm{k}$ is invariant under the two-fold (screw) rotation $(2)$ in $P2_{1}/c$ Hahn (2002) (note that our $x$ axis corresponds to their $y$ axis). This screw rotation further maps the field profile from one interface to the other, so that a 1D representation has to have intensities of equal magnitude on both interfaces. Since the two interfaces are well separated by a zero field bulk region, Maxwell’s equations are also satisfied for the same frequency by a field that is non-zero at one of the two interfaces only. The mode must thus be two-fold degenerate. Close to the $\Gamma$ point, the decay length becomes larger than $6$ unit cells, so that the argument breaks down and the degeneracy is lifted. For $X-\Gamma$, the same line of thought applies to the glide plane $(4)$.

Fig. 3(b) demonstrates that the modes within the bulk gap are indeed localized at the surface, in contrast to the bulk modes within the blue region. There are two degenerate mode pairs at $\bm{k}=0.2\pi/a\times(1,0,0)^{T}$ (brown points $2$ and $3$ in Fig. 3(a), only one states is shown, respectively). The degeneracy splits for $\bm{k}=0.2\pi/a\times(\cos(\phi),\sin(\phi),0)^{T}$ with arbitrary $\phi$ (green points for $\phi=0.28\pi$).

We have thus shown that surface states exist. But are these also of topological nature? The conventional path $\Gamma{-}Z{-}A{-}X{-}\Gamma$ does not reveal the topological nature of the four surface states emanating from the PWP. To show that these are, indeed, protected, we follow Wan et al. (2011) and consider the cylinder $\bm{k}(\varphi,k_{z})=(k\cos(\varphi),k\sin(\varphi),k_{z})^{T}$ (with constant $k$ and $-\pi/a<k_{z}\leq\pi/a$, $0<\varphi\leq 2\pi$). This cylinder constitutes a closed surface (a torus) in $\bm{k}$-space within which the bandstructure exhibits a band gap, so that a gap Chern number (sum over all bands below the gap) is well defined. The change in gap Chern number $|\Delta C|$ across an interface equals the number of topologically protected surface states that connect the bulk bands below the gap with those above Raghu and Haldane (2008); Lu et al. (2014a). The gap Chern number for the above torus and a hyperconic band at a PWP is given by $|C|=2$, as shown above (note that only the Chern number of the PWP at the gap frequency needs to be considered). This results in $8$ surface bands for the supercell geometry with two $|\Delta C|=4$ interfaces, four of which are observed along the half cylinder in Fig. 3(a): each of these bands touches and connects the projected bulk bands above and below the gap and thus is, veritably, protected.

In this letter, we have shown that isotropic hyperconic dispersion can be found at the $\Gamma$-point of chiral cubic photonic lattices. The associated exceptional PWP points exhibit the topological characteristics of a double Weyl point. A natural application for the unique dispersion behaviour of these PWPs are zero refractive index materials that have been suggested previously in the context of accidental Dirac points in two-dimensional photonic crystals Huang et al. (2011). The localization of associated protected surface states and their flatness yield a gigantic density of states, making PWP systems an ideal starting point to explore topological lasing applications Harari et al. (2016) in three dimensions.

I The canonical representation of the plasmonic sphere packing at $\Gamma$

The existence of three-fold degeneracies at the $\Gamma$ point in chiral tetrahedral and octahedral space groups is enumerated in the literature Bradley and Cracknell (2010). This general results stems from the fact that space group representations at $\Gamma$ can be written as a direct product between the representation of the isogonal point group and the trivial representation of the translation group, i.e.

[TABLE]

It is straightforward to verify that $\Delta_{\bm{k}=0,i}$ is an irreducible representation of the space group if $D_{i}$ is an irreducible representation of the isogonal point group. There are therefore 3D irreducible representations ($T$) for any chiral cubic space group associated with $\bm{k}=0$ (rigorous analysis Bradley and Cracknell (2010) shows that there are no others apart from the above): we tabulate the characters for tetrahedral point group in table 1, and refer to Saba et al. (2013) for the octahedral group. This generic finding is very useful for the application of the perturbation theory below. Additionally, we here derive the canonical representation of the plasmonic sphere packing and thereby the exact split of the underlying $12$ dimensional vector space into irreducible representations.

The representation matrix of the vector space $V$ corresponding to different spheres (within the equivalence class of translations) is given by permutation matrices. The following table lists the spheres $\nu$ onto which a sphere $\mu$ is mapped under the symmetry operation $\alpha$ in $P2_{1}3$. The symmetry indices are adapted from Hahn (2002).

[TABLE]

The space within which the polarisation vector is rotated is the three-dimensional Euclidean space $W$. We only list the relevant trace (or character) of the rotation matrix $\chi(\phi)=1+2\cos(\phi)$ here, where $\phi$ denotes the angle of rotation. The characters of the symmetry operations in $V$ and $W$, and of the canonical representation for the system vector space $V\otimes W$ are summarised below.

[TABLE]

The characters of the irreducible representations of the associated tetrahedral point group are shown in Tab. 1. Note that we use the Schoenflies notation $23$ in this section in order to avoid confusion with the symbol $T$ that is conventionally used for three-dimensional irreducible representations. The canonical representation at the $\Gamma$ point reduces into

[TABLE]

(where $(+)$ denotes a direct sum operation in this context) according to the reduction formula Dresselhaus (2013)

[TABLE]

where $n_{i}$ quantifies how often an irreducible representation $i$ is contained in the reducible representation with characters $\chi(R)$ ($\,\overline{\cdot\vphantom{\bar{\cdot}}}\,$ denotes complex conjugation), $N$ is the number of elements $R$ in the group, and $\chi_{i}(R)$ is the character of the irreducible representation.

Note that the representation $E_{\pm}$ form a two-fold degenerate pair $E$ because the system is also invariant under time inversion, as evident from the Schur type $(b)$ in Tab. 1 here. For details on time inversion symmetry in addition to space group symmetries see for example Herring Herring (1937) and Frobenius Frobenius (1968), page 354 ff. The predicted split is thus reproduced in Fig. 2 in the main manuscript: the bandstructure shows one trivial $A$ representation, one 2D $E$ representation, and three 3D $T$ representations, two of which produce genuine hyperconic dispersion, while the third is degenerate in first order and corresponds to the special case $d=0$ below.

II Perturbation theory and selection rules

Reduction into the irreducible representations of the respective subgroup of the space group $\mathcal{G}$ as described in Saba et al. (2013) predicts how and whether a degenerate mode splits when going away from that point. If not perturbed along a high symmetry direction in reciprocal space it always splits into 1D representations as the Bloch representations, characterized by $\bm{k}$, with respect to the invariant subgroup of translations generate a full star, with the same dimension as the underlying point group.

We are here however concerned with the order of splitting: for linear Weyl dispersion, the first order cannot vanish. We here derive an approach that is analogous to, but more general than the well established $\bm{k}\cdot\bm{p}$ perturbation theory Dresselhaus (2013). Consider a generic non-linear eigenproblem

[TABLE]

characterized by an operator $M(\lambda)$ that commutes with all operators of $\mathcal{G}$. The eigenvectors $\bm{v}$ are partners $(\bm{k}_{n},\alpha)$ of the irreducible representations of the space group enumerated by $(\bm{k},i)$. In this context, $\bm{k}$ is a Bloch wave vector in the irreducible Brillouin zone and $i$ denotes different induced representations from the small representations of the associated little group, whereas $\bm{k}_{n}$ iterates over different representatives in the star of $\bm{k}$ and $\alpha$ over partners of the small representations of the little group. We hence write:

[TABLE]

Working in the spatial Fourier domain, we can suppress the explicit $\bm{k}_{n}$ dependence in the eigenvector, dividing by the phase factor $\exp\left\{\imath\bm{k}_{n}\cdot\bm{r}\right\}$ corresponding to a translation group representation. At an arbitrary point within the BZ we define

[TABLE]

Importantly, these are partners $\alpha$ of the irreducible representations of the little group at $\bm{k}_{n}$. The eigenvalue equation becomes

[TABLE]

For $\bm{k}_{n}^{(1)}=\bm{k}_{n}^{(0)}+\delta\bm{k}$, with small $|\delta\bm{k}|\ll\pi/a$, and hence $\lambda_{1}=\lambda_{0}+\delta\lambda$, a first order Taylor expansion of $M$ around $M(\bm{k}_{n}^{(0)},\lambda_{0})$ and an expansion of the eigenfunctions into the aforementioned partners $\alpha$ (cf. representation theorem (iii) in Saba et al. (2013)) yields:

[TABLE]

where we have assumed that $i$ labels the irreducible representation corresponding to the eigenvalue $\lambda_{0}$ and used the fact that $M(\bm{k}_{n}^{(0)},\lambda_{0})|\bm{k}j,n\beta\rangle\neq 0$ for any representation $j\neq i$ (i.e. we exclude accidental degeneracies), so that a finite coefficient would contradict Eq. 6 for arbitrary but small $\delta\bm{k}$. In the weak form, Eq. 6 reads:

[TABLE]

with

[TABLE]

and

[TABLE]

The partial derivative $\partial\lambda$ is trivially invariant under all operations $P$ of the little group at $\bm{k}_{n}^{(0)}$. $M$ is invariant under all space group operations, and $P\in\mathcal{G}$, while $P|\bm{k}i,n\alpha\rangle=\sum_{\beta}D_{\alpha\beta}(P)|\bm{k}i,n\beta\rangle$. Hence, the action on $W$ and $E$ is generally given by

[TABLE]

where $R(P)$ denotes the Rodrigues’ matrix representation that acts in the standard 3D Euclidean vector space. Since $\delta\bm{k}$ is arbitrary, according to Eq. 7, $\bm{W}$ and $E$ are also required to be invariant under all $P$. Introducing the respective direct product space, $E$ and $W$ must hence for all $P$ be in the kernel of $\Delta(P)$ and $\Lambda(P)$, respectively, where

[TABLE]

This reveals the group theoretically allowed form of $E$ and $W$, substantially reducing their degrees of freedom and providing selection rules for many of their entries. $E$ and $W$ are respectively non-zero only if $\Delta+{\mathbbm{1}}$ and $\Lambda+{\mathbbm{1}}$ contain the trivial representation if decomposed into the irreducible representations of $i$ via Eq. 4 (representation theorem (ii) in Saba et al. (2013)), i.e.

[TABLE]

More precisely, the integer number on the left side of Eq. 10 and Eq. 11 is equal to the dimension of the intersection of the nullspaces of $\Delta$ and $\Lambda$, respectively, and hence specifies the degrees of freedom in $E$ and $W$ (but not their form).

III Perturbation theory at $\Gamma$

Consider now a perturbation around the centre of the BZ using the equations derived in the previous section. The isogonal little group is here the isogonal point group of the lattice, i.e. for chiral cubic symmetry either tetrahedral ($23$) or octahedral ($432$). We here discuss the tetrahedral symmetry, but note that the main results regarding PWPs are exactly the same for the octahedral point group. Consider first a non-degenerate mode with trivial irreducible representation $A$ (see table 1). While Eq. 10 is satisfied, Eq. 11 is not. $W$ is therefore trivial, with all its elements vanishing and Eq. 7 yielding $\delta\lambda=0$. A non-degenerate mode is henceforth flat at the centre of the BZ in a cubic symmetry, even in the absence of inversion or time reversal symmetry. Note, that the $E_{\pm}$ modes show the same behaviour, so that a deterministic two-fold degenerate mode, as present in the case of time reversal only in tetrahedral symmetries is also predicted to be flat at the centre of the BZ. Further note that the two-fold $E$ representation in the octahedral point group (table 3.33 in Dresselhaus (2013)) also yields a flat band.

The representation $T$ in table 1 (the Rodrigues’ representation itself) that corresponds to a three-fold degeneracy, however, satisfies both Eq. 10 and Eq. 11, with $\text{dim}(\text{ker}(\Delta))=1$ and $\text{dim}(\text{ker}(\Lambda))=2$. We compute the intersection of the kernels of $\Delta(P)$ and $\Lambda(P)$ for all $P$ numerically, and find explicitly for the $T$ degeneracy that $E$ is a scalar matrix

[TABLE]

and that $\bm{W}$ is given by

[TABLE]

In the presence of time reversal symmetry we further require that $\,\overline{M\vphantom{\bar{M}}}\,(-\bm{k})=M(\bm{k})$, and hence $c\in\mathbb{R}$ and $W_{\gamma}$ Hermitian, so that Eq. 7 corresponds to the perturbation matrix given in eq. $1$ of the main manuscript:

[TABLE]

where we set w.l.o.g. $d_{1}=-\imath d$, $d_{2}=\imath d$ with $d\in\mathbb{R}$ and $c=-1$.

IV Chern numbers for the hyperconic bands

We here calculate the Chern numbers for the three bands that meet at the PWP. We start from the perturbation matrix (Eq. $1$ in the main manuscript, derived above)

[TABLE]

To study the topological nature of the PWP, it is convenient to calculate the Chern number over a small sphere in reciprocal space, with the PWP at its centre. We first introduce spherical coordinates with $\bm{k}:=k(\cos{\phi}\sin{\theta},\sin{\phi}\sin{\theta},\cos{\theta})^{T}$ and associated basis vectors $\{\bm{\hat{k}},\bm{\hat{\theta}},\bm{\hat{\phi}}\}$. With the eigenvectors expanded as $\bm{v}=v_{k}\bm{\hat{k}}+v_{\theta}\bm{\hat{\theta}}+v_{\phi}\bm{\hat{\phi}}$, the perturbation equation

[TABLE]

has the solution set

[TABLE]

We note that $\bm{v}$ above is the vector over the coefficients $c_{\alpha}$ of partners $|\alpha\rangle$ of the irreducible $T$ representation. These are by definition orthonormalized and $\bm{k}$ independent. Topologically, the discussion of any physical Hilbert space therefore collapses to the three-dimensional vector space of the coefficients. The Gaussian curvature

[TABLE]

is a measure of intrinsic curvature of the tangent space of a surface spanned by local basis vectors $\bm{e}_{1}$ and $\bm{e}_{2}$. The Berry curvature is a similar measure concerning a generic complex vector field $|v(\bm{x})\rangle$ that is attached to a manifold $\mathcal{M}$ parametrized by $\bm{x}$. It can be conveniently written as the two-form

[TABLE]

with the projector $P:=|v\rangle\langle v|$ and its exterior derivative $dP=\sum_{i}dx_{i}\partial_{x_{i}}P$ in the tangent space of $\mathcal{M}$. Note that the wedge product does not vanish as $P$ is an operator, but the Berry curvature is easily shown to be real. In our case, the manifold is 2D, with $x_{1}=\theta$ and $x_{2}=\phi$, so that the Berry curvature reduces to

[TABLE]

Similar to the Euler number $\chi:=(2\pi)^{-1}\int_{\mathcal{M}}K$ that approaches integer values for any closed surface and is directly related to its genus, the Chern characteristics

[TABLE]

is also quantized and constitutes a topological invariant of the vector field. The flat band eigenvector $\bm{v}=\bm{\hat{k}}$ is real and trivially leads to a vanishing Berry curvature when substituted into Eq. 12. The eigenvectors of the hyperconic bands $\bm{v}_{\pm}=(\bm{\hat{\theta}}\pm\imath\bm{\hat{\phi}})/\sqrt{2}$ on the other hand yield a Berry curvature of

[TABLE]

and hence a Chern characteristics of $C=\pm 2$, in direct analogy to the Euler number $\chi=2$ of the sphere.

V Additional figures

Fig. 4 shows the plain bulk bandstructure of the photonic crystal sphere packing close to the PWP of interest around $k_{0}a/(2\pi)=0.5$. In Fig. 5 the sphere positions in the unit cell are distorted, thereby destroying cubic symmetry and both degeneracies at the $\Gamma$ and at the $R$ point. We note that, despite the fact that the Bravais lattice remains simple cubic by construction, the path followed in Fig. 5 is thus in principle inappropriate, but kept for the sake of comparability with Fig. 4. A full bandgap is, however, observed for the particular realization in the whole (triclinic) Brillouin zone. Opening the gap does not change the topological characteristics of the bands above and below the gap. Frequency isolated protected surface states can hence be found for the perturbed structure.

VI The eigenvalue equation for the sphere packing

Since a proper derivation of Eq. 1 of the main manuscript, reproduced here for clarity

[TABLE]

seems to be absent in the literature, we here rigorously show its validity in the assumed limits. We model the dielectric matrix intrinsically via an effective displacement field with the corresponding permittivity $\varepsilon_{d}$, and the polarisation $\bm{P}_{m}$ of the metallic spheres explicitly. The macroscopic Maxwell equations for monochromatic fields can then be transformed into a wave equation for the electric field $\bm{E}$ and the polarisation field $\bm{P}_{m}$, that vanishes in the dielectric matrix:

[TABLE]

We have here introduced the wave number $\kappa:=\sqrt{\varepsilon_{d}}k_{0}$ and the Maxwell operator $\mathcal{M}$ for convenience.

Eq. 14 can be inverted using the dyadic Green function

[TABLE]

that solves the equation $\mathcal{M}\mathcal{G}_{\mathcal{M}}(r)={\mathbbm{1}}\delta^{(3)}(\bm{r})$, with the (retarded) Helmholtz Green function $\mathcal{G}_{\mathcal{H}}(r)=-(4\pi r)^{-1}\exp\left\{\imath\kappa r\right\}$ that solves

[TABLE]

We assume that the metallic material response $\varepsilon_{m}(k_{0})$ is known and frequency dependent. Therefore, with $\bm{E}=4\pi/(\varepsilon_{m}-\varepsilon_{d})\,\bm{P}_{m}$, Eq. 14 transforms into a relation on the metallic subdomain only:

[TABLE]

The polarisation field inside a single sphere can be considered constant for sufficiently small frequencies, in zero order in $k_{0}\rho\ll 1$.111Note that this condition does not impose a strong constraint on $k_{0}$ itself, since we assume small spheres $\rho\ll\min(d_{\mu\nu})<a$. We hence make the ansatz

[TABLE]

and $V_{S}=(4\pi\rho^{3})/3$ the sphere volume. Substitution into Eq. 16 and testing with the basis functions $\chi(\bm{r}-\bm{r}_{i})$ yields an equation for the dipole moments $\bm{p}_{i}$ only:

[TABLE]

Here,

[TABLE]

is a double integral over balls $\mathcal{B}_{\rho}(\bm{r})$ with radius $\rho$ at position $\bm{r}$. For any two spheres $i\neq j$ and in the limit of small spheres $\rho\ll d_{ij}$, $I_{\rho}(\bm{r}_{i},\bm{r}_{j})$ reduces to $V_{S}\,\mathcal{G}_{\mathcal{M}}(\bm{r}_{i}-\bm{r}_{j})$. For $i=j$, we use the following theorem.

Theorem I: Let

[TABLE]

be a double integral over balls (at the origin) $\mathcal{B}_{\rho}$ of a function $f(\bm{x})$ that only depends on the difference in coordinates $\bm{x}=\bm{r}-\bm{r}^{\prime}$. The integral can then be expressed as a single integral over a Ball $\mathcal{B}_{2\rho}$

[TABLE]

Proof: A change of coordinates $(\bm{r},\bm{r}^{\prime})\mapsto(\bm{x},\bm{R}:=\bm{r}+\bm{r}^{\prime})$ yields:

[TABLE]

where, with the binary operator $(+)$ on spatial domains denoting the Minkowski sum,

[TABLE]

Using the translation invariance of the volume and a mixed volume expansion Burago and Zalgaller (1988), we arrive at:

[TABLE]

In the last line we used that the volume is twice that of a spherical cap with sphere radius $\rho$ and height $\rho-x/2$. The proposition follows from resubstitution. $\square$

Applying theorem I with $f(\bm{x})=\mathcal{G}_{\mathcal{M}}(\bm{x})$ and the low frequency condition on the Green function yields the simplified expression

[TABLE]

with only two independent entries $\mathcal{G}_{\mathcal{M}}^{(zz)}$ and $\mathcal{G}_{\mathcal{M}}^{(xz)}$. The integral over the latter entry vanishes as the integrand is antisymmetric in both the azimuth $\phi$ and the polar angle $\theta$. The former integral evaluates to the simple expression

[TABLE]

Substitution of these results into Eq. 17 yields

[TABLE]

For spheres made of a lossless Drude material $\varepsilon_{m}(k_{0})=1-k_{p}^{2}/k_{0}^{2}$, embedded in vacuum $\varepsilon_{d}=1$, this equation is identical to Eq. 13. We note that this describes an assembly of point dipoles with polarisability $\alpha$, that are interacting via the Maxwell Greens function without self-interaction. This simple picture apparently describes the correct physics in the assumed limits, although it is not possible to derive Eq. 13 within its scope.

VII The lattice sum

Substituting the Bloch form $\bm{p}_{\bm{n},\mu}=\bm{p}_{\mu}\exp\left\{\imath\bm{k}\cdot\bm{T}_{\bm{n}}\right\}$ into Eq. 13 yields

[TABLE]

where we define $\mathcal{G}(0):=0$ to simplify the notation.

Eq. 2 from the main manuscript, namely

[TABLE]

results from a lattice vector shift $\bm{T}_{m}\mapsto\bm{T}_{m}+\bm{T}_{\bm{n}}$, with the matrix

[TABLE]

This numerical lattice sum that is truncated by for example $|\bm{n}|<N$ converges very slowly to the infinite series Eq. 20.222Precisely speaking, it does not converge at all in the vicinity of Ewald’s sphere $|\bm{k}|=k_{0}$ due to a ringing like phenomenon (spatial and temporal dispersion counteract and prevent the summands from alternating in one direction). Convergence is generally guaranteed by the alternating phase factor, but the rate of convergence is governed by the long range $r^{-3}$ interaction from the Green function, and can only be algebraic. This problem is well known and dates back to Ewald Ewald (1921), and also appears in the study of the electronic structure in crystals Hove and Krumhansl (1953).

In order to circumvent the numerically expensive and inaccurate real space representation in Eq. 20, we follow an approach introduced by Silveirinha and Fernandes Silveirinha and Fernandes (2005) for the Helmholtz Green function $\mathcal{G}_{\mathcal{H}}(\bm{r})$. Their main idea is to split the lattice sum into real space contribution and a spectral contribution, using a radial weight function $f(r)$ that approaches $1$ at exponential (or even Gaussian) rate for $r\rightarrow\infty$ and vanishes at the origin $r=0$, thereby splitting $M=Mf(r)+M[1-f(r)]$. The first term is then effectively calculated in real space, whereas the second term converges rapidly in the spectral representation, for which a plane wave expansion yields a summation over reciprocal vectors $\bm{G}$ instead of the real space vectors $\bm{T}$.

This technique can be trivially generalised to the dyadic case exploiting Eq. 15. To highlight the importance of this method, we note that the procedure enables us to obtain converged solutions, with relative errors in the eigenvalues at machine precision, for cut-offs as small as $N=5$. For comparison, the naive representation Eq. 20 yields relative errors as large as $10^{-2}$ for $N=100$. Under these circumstances, numerical summation over almost $10^{7}$ elements per Bloch wave vector makes band structure calculations slow, while still inaccurate and suffering from ringing artefacts. On the contrary, the eigenvalue problem itself is of dimension $12$ in our example and comes at virtually no computational cost.

VIII Hamilton formulation and energy flow

Consider the coupled harmonic oscillator Hamiltonian

[TABLE]

with the dipole moments $\bar{\bm{p}}_{n}$ taking the role of canonical coordinates, with conjugated momenta $\bar{\bm{\Pi}}_{n}=\frac{d\bar{\bm{p}}_{n}}{dt}$.333Note that the dipole moments $\bm{p}_{n}$ in Eq. 13 denote the phasors of monochromatic physical dipole moments $\bar{\bm{p}}_{n}$. Further note that the above Hamiltonian has the odd dimension $(\text{charge}\times\text{length}/\text{time})^{2}$. This has no direct implication on the following, particularly on the dimensions of the mode velocity. However, scaling the coordinates with $\sqrt{m}/e$, with the effective mass $m$ and the charge $e$ of the microscopic charges, yields correct physical dimensions. The canonical coordinates and momenta are proportional to the centre of charge with respect to the corresponding sphere centre and the average charge flow, respectively (see also Weick et al. (2013)).

With the single particle resonance $\omega_{r}:=\omega_{p}/\sqrt{3}$ and the interaction matrix given by $I_{nm}:=-\rho^{3}\mathcal{G}(\bm{r}_{n}-\bm{r}_{m})$ for $n\neq m$ and $I_{nn}=0$, it is straightforward to show that the equations of motion of this Hamiltonian reduce to Eq. 13 in the monochromatic case.

The Hamilton formulation has the advantage, that we can easily derive expressions for the average energy flow of the plasmonic modes. The single particle energy $E_{i}=\frac{1}{2}\left(\bar{\bm{\Pi}}_{i}^{2}+\omega_{r}^{2}\bar{\bm{p}}_{i}^{2}\right)$ has the time derivative

[TABLE]

That means that the amount of energy that flows from point $i$ to point $j$ per time is given by $-\omega_{r}^{2}\bar{\bm{\Pi}}_{i}^{T}I_{ij}\,\bar{\bm{p}}_{j}$. At the same time, $-\omega_{r}^{2}\bar{\bm{\Pi}}_{j}^{T}I_{ji}\,\bar{\bm{p}}_{i}=-\omega_{r}^{2}\bar{\bm{p}}_{i}^{T}I_{ij}\,\bar{\bm{\Pi}}_{j}$ obviously flows from $j$ to $i$. The total power that is transferred from $i$ to $j$ is hence given by the difference of these two contributions (as physically required, it is antisymmetric by construction):

[TABLE]

The symmetric energy flow between $i$ and $j$ can now be defined as

[TABLE]

The total energy flow of a system of spheres, averaged over time for the monochromatic field is then given by

[TABLE]

whereas the time average of the total energy is

[TABLE]

The mode velocity is then given by

[TABLE]

For the periodic system, we can replace $\sum_{n}\mapsto N\sum_{\mu}$ and $\sum_{nm}\mapsto N\sum_{\mu\nu}$ in Eq. 23, if we also replace $I_{nm}\mapsto I_{\mu\nu}:=-\rho^{3}M_{\mu\nu}$ and $\bm{p}_{n}\mapsto\bm{p}_{\mu}$. The procedure is similar to the one that led to Eq. 20 above. For the flow, we have

[TABLE]

so that, in the lossless Drude case, the average energy velocity reduces to the simple expression for the group velocity, as is well known for photonic crystals Joannopoulos et al. (2008)

[TABLE]

where we have used Eq. 19 and the Hellmann-Feynman theorem. Note that for any more sophisticated material model $\alpha(k_{0})$, the substitution of Eq. 22 and Eq. 23 into Eq. 24 yields the general form

[TABLE]

for which we implicitly assume $\alpha^{-1}=\alpha^{-1}(K_{0}(\bm{k}))$. This expression does not reduce to the group velocity Eq. 25, even for a more general lossless dispersion.