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

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.

Current broad interest in topological phases, triggered by the discovery of the quantum Hall effect [1] and its theoretical investigation [2][3][4], 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 [5,6]. Plasmonic [7] and single electron [8] 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 [9,10], sparking a new wave of research on photonic topology [11]. In particular, topologically protected Weyl points with hyperconic dispersion have been identified in double gyroid photonic crystals with broken parity-time symmetry [12]. 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 [13] and quantum mechanical [14,15] systems. Indeed, the existence of deterministic two and three-fold degeneracies at the center of the Brillouin zone (BZ), aka the Γ point, for cubic symmetries is well known and documented in the literature [16]. Recently, it has been shown that some of these degeneracies are topologically non-trivial in electronic systems [15,18].
Here we show on the basis of group and perturbation theory that symmetry induced three-fold degenerate pseudo Weyl points (PWPs) exist at the Γ point in classical (photonic) systems. They split isotropically in first order in the Bloch wave vector 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 [19]. 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 * Email: m.saba@imperial.ac.uk † Email: o.hess@imperial.ac.uk dispersion with non-trivial topology, and an intermediate flat band. We then construct a minimalistic geometry, a P 2 1 3 sphere packing (Fig. 1), which satisfies the symmetry requirements, and apply it to a quasistatic coupleddipole 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 in-arXiv:1706.06030v1 [physics.optics] 19 Jun 2017 tegral over Hilbert states |v(ω) ∈ H which individually solve a homogeneous (generally non-linear) frequency domain eigenproblem M (ω) |v(ω) = 0. Symmetry requires that the H operator M commutes with all elements of the underlying space group G, represented by H operators g. As a consequence, a set of N degenerate eigenvectors |v n (ω) form a representation of G, i.e. they span an N -dimensional vector space with an associated algebra that is homomorphic to G. Irreducible representations impose a lower limit on the dimensionality of the respective vector space, resulting in deterministic degeneracies [16,20,21].
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 Γ point, this procedure yields a perturbation matrix W (k) which is valid for small k 2π/a (with lattice constant a) and time reversal invariance with a free parameter d ∈ R and α, β, γ iterating over the three partners of the irreducible representation with k = 0 that span the degenerate eigenspace. Note the similarity of W αβ to the Weyl Hamiltonian W ij : the Pauli matrices σ ) that occur in the latter are here replaced by the 3D Levi-Civita tensor αβγ . The first order perturbation eigenvalues corresponding to W (k) are k (1) 0 := ω (1) /c = {0, ±dk}: they only depend on the absolute value of 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 (k, k (0) 0 ± dk) in the four dimensional (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 k space for each of the three bands. They evaluate to C = 0 for the flat band and C = ±2 for the two hyperconic bands, showing a non-trivial topological signature, similar to a genuine Weyl point with Chern numbers C = ±1.
Analysing all 3D space groups [17], it is straightforward to show that deterministic PWPs at the center of the BZ require chiral cubic symmetry. Interestingly, the trigonal groups P 312 (149) and P 321 (150) have two-dimensional representations which split into an anisotropic hypercone if time inversion is present, albeit not at the Γ point [23]. A closely related matter is the non-existence of deterministic Dirac points at the Γ point of two-dimensional crystals, including the famous honeycomb lattice [13].
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 [17]. Following this procedure for any non-symmorphic chiral cubic space group G and Wyckoff multiplicities smaller than 12 yields a sphere packing that has the symmetry of an achiral supergroup 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 G = P 2 1 3 thus induces a chiral cubic sphere packing with G S = G and only 4 spheres in the unit cell. Only if x = n/4 (n ∈ Z), the sphere packing acquires the symmetry of the achiral supergroup P a3 (205); the chiral supergroup P 4 3 32 (212) is induced for x = 1/8 + n/2, and P 4 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 P 2 1 3 in these cases does not impose a change of bandstructure behaviour close to Γ. Fig. 1 illustrates the sphere packing for x = 3/8 a (P 4 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 ρ in vacuum (as in [24,25]). The position r i of sphere i shall be such that the distance d ij = |r i − r j | ρ for any pair of spheres (i, j). In the quasistatic approximation, Maxwell's equations take the self-consistent form (acting on the dipole moments p i ) [26]: Here, α(k 0 ) = ρ 3 (1 − 3k 2 0 /k 2 p ) −1 is the polarizability of a metallic sphere in vacuum, that is modelled by a nondissipative Drude response with plasma wave number k p ; G(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 (n, µ) ∈ Z 3 × {1, 2, 3, 4}, with r n,µ = T n + r µ given by the sum of the lattice vector T n = a n and the position within the unit cell r µ . Bloch's theorem then implies p n,µ = p µ exp{ık · T n }, so that Eq. 13 reduces to a family of 12-dimensional non-linear Hermitian eigenproblems: Numerical challenges related to the convergence of the lattice sum M µν (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 (n) 0 in Eq. 19), the eigenvalue problem is linearized by ap- . This assumption is inadequate close to the Ewald sphere k 0 = |k|, caused by poles in the diagonal entries of M (k = k 0k ), 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 λ n (k) = α −1 (k 0 ) (n = 1, 2, . . . , N ) of M (k) can be obtained numerically with low computational cost. They produce the respective dispersion relation K (n) 0 (k) = [1 − ρ 3 λ n (k)] 1/2 , as shown in Fig. 2 for parameters x/a = 0.175, k 0 a/(2π) = 0.1 and ρ/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 ≈ π/(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 Γ point for any space group, cubic or noncubic (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 ρ/a = 0.25 (fill fraction of π/12 ≈ 26%), made of a high refractive index material with n = 4, thus constructing a photonic crystal analog. The associated bandstructure (calculated with MPB [27]) close to k 0 a/(2π) = 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 (k , k 0 ) for which at least one bulk mode exists for arbitrary k z ∈ R [28]. 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 P 2 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 Γ − Z and X − Γ (including Γ itself) are 1D. The surface states still stick together. This can be understood as follows: consider a single interface surface mode along Γ − Z. Along this path, k is invariant under the two-fold (screw) rotation (2) in P 2 1 /c [17] (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 Γ point, the decay length becomes larger than 6 unit cells, so that the argument breaks down and the degeneracy is lifted. For X − Γ, 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 k = 0.2π/a × (1, 0, 0) T (brown   Fig. 3(a), only one states is shown, respectively). The degeneracy splits for k = 0.2π/a × (cos(φ), sin(φ), 0) T with arbitrary φ (green points for φ = 0.28π).
We have thus shown that surface states exist. But are these also of topological nature? The conventional path Γ−Z−A−X−Γ does not reveal the topological nature of the four surface states emanating from the PWP. To show that these are, indeed, protected, we follow [30] and consider the cylinder k(ϕ, k z ) = (k cos(ϕ), k sin(ϕ), k z ) T (with constant k and −π/a < k z ≤ π/a, 0 < ϕ ≤ 2π). This cylinder constitutes a closed surface (a torus) in kspace 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 |∆C| across an interface equals the number of topologically protected surface states that connect the bulk bands below the gap with those above [10,11]. 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 |∆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 Γ-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 [19]. 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 [31] in three dimensions.

ACKNOWLEDGMENTS
This work was supported by the EPSRC through program grant EP/L027151/1. We would like to thank Dr. Paloma Arroyo Huidobro, Gleb Siroki and Prof. Sir John Pendry for helpful discussions. The existence of three-fold degeneracies at the Γ point in chiral tetrahedral and octahedral space groups is enumerated in the literature [16]. This general results stems from the fact that space group representations at Γ 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.

Supplementary Material
It is straightforward to verify that ∆ 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 k = 0 (rigorous analysis [16] shows that there are no others apart from the above): we tabulate the characters for tetrahedral point group in table I, and refer to [21] 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 ν onto which a sphere µ is mapped under the symmetry operation α in P 2 1 3. The symmetry indices are adapted from [17]. 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 χ(φ) = 1 + 2 cos(φ) here, where φ 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 ⊗W are summarised below.
The characters of the irreducible representations of the associated tetrahedral point group are shown in Tab. I. 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 Γ point reduces into (where (+) denotes a direct sum operation in this context) according to the reduction formula [20] where n i quantifies how often an irreducible representation i is contained in the reducible representation with characters χ(R) ( · denotes complex conjugation), N is the number of elements R in the group, and χ i (R) is the character of the irreducible representation. Note that the representation E ± 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. I here. For details on time inversion symmetry in addition to space group symmetries see for example Herring [32] and Frobenius [33], 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 G as described in [21] 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 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 k · p perturbation theory [20]. Consider a generic non-linear eigenproblem characterized by an operator M (λ) that commutes with all operators of G. The eigenvectors v are partners (k n , α) of the irreducible representations of the space group enumerated by (k, i). In this context, 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 k n iterates over different representatives in the star of k and α over partners of the small representations of the little group. We hence write: Working in the spatial Fourier domain, we can suppress the explicit k n dependence in the eigenvector, dividing by the phase factor exp {ık n · r} corresponding to a translation group representation. At an arbitrary point within the BZ we define Importantly, these are partners α of the irreducible representations of the little group at k n . The eigenvalue equation becomes M (k n , λ)|ki, nα = 0 .
For k n + δk, with small |δk| π/a, and hence λ 1 = λ 0 + δλ, a first order Taylor expansion of M around M (k (0) n , λ 0 ) and an expansion of the eigenfunctions into the aforementioned partners α (cf. representation theorem (iii) in [21]) yields: where we have assumed that i labels the irreducible representation corresponding to the eigenvalue λ 0 and used the fact that M (k (0) n , λ 0 )|kj, nβ = 0 for any representation j = i (i.e. we exclude accidental degeneracies), so that a finite coefficient would contradict Eq. 6 for arbitrary but small δk. In the weak form, Eq. 6 reads: with The partial derivative ∂λ is trivially invariant under all operations P of the little group at k n . M is invariant under all space group operations, and P ∈ G, while P |ki, nα = β D αβ (P )|ki, nβ . Hence, the action on W and E is generally given by where R(P ) denotes the Rodrigues' matrix representation that acts in the standard 3D Euclidean vector space. Since δk is arbitrary, according to Eq. 7, 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 ∆(P ) and Λ(P ), respectively, where ∆ 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 ∆ and Λ, respectively, and hence specifies the degrees of freedom in E and W (but not their form).

III. PERTURBATION THEORY AT Γ
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 I). While Eq. 10 is satisfied, Eq. 11 is not. W is therefore trivial, with all its elements vanishing and Eq. 7 yielding δλ = 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 ± 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 [20]) also yields a flat band.
The representation T in table I (the Rodrigues' representation itself) that corresponds to a three-fold degeneracy, however, satisfies both Eq. 10 and Eq. 11, with dim(ker(∆)) = 1 and dim(ker(Λ)) = 2. We compute the intersection of the kernels of ∆(P ) and Λ(P ) for all P numerically, and find explicitly for the T degeneracy that E is a scalar matrix and that W is given by In the presence of time reversal symmetry we further require that M (−k) = M (k), and hence c ∈ R and W γ Hermitian, so that Eq. 7 corresponds to the perturbation matrix given in eq. 1 of the main manuscript: where we set w.l.o.g. d 1 = −ıd, d 2 = ıd with d ∈ 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) 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 k := k(cos φ sin θ, sin φ sin θ, cos θ) T and associated basis vectors {k,θ,φ}. With the eigenvectors expanded as v = v kk + v θθ + v φφ , the perturbation equation We note that v above is the vector over the coefficients c α of partners |α of the irreducible T representation. These are by definition orthonormalized and k independent. Topologically, the discussion of any physical Hilbert space therefore collapses to the three-dimensional vector space of the coefficients. The Gaussian curvature is a measure of intrinsic curvature of the tangent space of a surface spanned by local basis vectors e 1 and e 2 . The Berry curvature is a similar measure concerning a generic complex vector field |v(x) that is attached to a manifold M parametrized by x. It can be conveniently written as the two-form with the projector P := |v v| and its exterior derivative dP = i dx i ∂ xi P in the tangent space of 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 = θ and x 2 = φ, so that the Berry curvature reduces to Similar to the Euler number χ := (2π) −1 M K that approaches integer values for any closed surface and is directly related to its genus, the Chern characteristics is also quantized and constitutes a topological invariant of the vector field. The flat band eigenvector v =k is real and trivially leads to a vanishing Berry curvature when substituted into Eq. 12. The eigenvectors of the hyperconic bands v ± = (θ ± ıφ)/ √ 2 on the other hand yield a Berry curvature of B = ± dθ ∧ dφ ∂ θθ , ∂ φφ = ± sin(θ) dθ ∧ dφ , and hence a Chern characteristics of C = ±2, in direct analogy to the Euler number χ = 2 of the sphere. Fig. 4 shows the plain bulk bandstructure of the photonic crystal sphere packing close to the PWP of interest around k 0 a/(2π) = 0.5. In Fig. 5 the sphere positions in the unit cell are distorted, thereby destroying cubic symmetry and both degeneracies at the Γ 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 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 ε d , and the polarisation 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 E and the polarisation field P m , that vanishes in the dielectric matrix: We have here introduced the wave number κ := √ ε d k 0 and the Maxwell operator M for convenience. Eq. 14 can be inverted using the dyadic Green function We assume that the metallic material response ε m (k 0 ) is known and frequency dependent. Therefore, with E = 4π/(ε m − ε d ) P m , Eq. 14 transforms into a relation on the metallic subdomain only: The polarisation field inside a single sphere can be considered constant for sufficiently small frequencies, in zero order in k 0 ρ 1. 1 We hence make the ansatz and V S = (4πρ 3 )/3 the sphere volume. Substitution into Eq. 16 and testing with the basis functions χ(r − r i ) yields an equation for the dipole moments p i only: Here, is a double integral over balls B ρ (r) with radius ρ at position r. For any two spheres i = j and in the limit of small spheres ρ d ij , I ρ (r i , r j ) reduces to V S G M (r i − r j ). For i = j, we use the following theorem. Theorem I: Let be a double integral over balls (at the origin) B ρ of a function f (x) that only depends on the difference in coordinates x = r − r . The integral can then be expressed as a single integral over a Ball B 2ρ Proof: A change of coordinates (r, r ) → (x, R := r + r ) yields: , where, with the binary operator (+) on spatial domains denoting the Minkowski sum, This numerical lattice sum that is truncated by for example |n| < N converges very slowly to the infinite series Eq. 20. 2 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 [35], and also appears in the study of the electronic structure in crystals [36].
In order to circumvent the numerically expensive and inaccurate real space representation in Eq. 20, we follow an approach introduced by Silveirinha and Fernandes [37] for the Helmholtz Green function G H (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 → ∞ and vanishes at the origin r = 0, thereby splitting M = M f (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 G instead of the real space vectors 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 with the dipole momentsp n taking the role of canonical coordinates, with conjugated momentaΠ n = dpn dt . 3 The I µν := −ρ 3 M µν and p n → p µ . The procedure is similar to the one that led to Eq. 20 above. For the flow, we have 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 [28] v = ω r ∇ k K 0 (k) = ∇ k ω(k) , where we have used Eq. 19 and the Hellmann-Feynman theorem. Note that for any more sophisticated material model α(k 0 ), the substitution of Eq. 22 and Eq. 23 into Eq. 24 yields the general form for which we implicitly assume α −1 = α −1 (K 0 (k)). This expression does not reduce to the group velocity Eq. 25, even for a more general lossless dispersion.