"Cartesian light": unconventional propagation of light in a 3D superlattice of coupled cavities within a 3D photonic band gap

We explore the unconventional propagation of light in a three-dimensional (3D) superlattice of coupled resonant cavities in a 3D photonic band gap crystal. Such a 3D cavity superlattice is the photonic analogue of the Anderson model for spins and electrons in the limit of zero disorder. Using the plane-wave expansion method, we calculate the dispersion relations of the 3D cavity superlattice with the cubic inverse woodpile structure that reveal five coupled-cavity bands, typical of quadrupole-like resonances. For three out of five bands, we observe that the dispersion bandwidth is significantly larger in the $(k_x, k_z)$-diagonal directions than in other directions. To explain the directionality of the dispersion bandwidth, we employ the tight-binding method from which we derive coupling coefficients in 3D. For all converged coupled-cavity bands, we find that light hops predominantly in a few high-symmetry directions including the Cartesian $(x, y, z)$ directions, therefore we propose the name"Cartesian light". Such 3D Cartesian hopping of light in a band gap yields propagation as superlattice Bloch modes that differ fundamentally from the conventional 3D spatially-extended Bloch wave propagation in crystals, from light tunneling through a band gap, from coupled-resonator optical waveguiding, and also from light diffusing at the edge of a gap.


I. INTRODUCTION
Ever since the emergence of the field of nanophotonics, it is well-known that fruitful analogies can be drawn between the behavior of photons at the nanoscale on the one hand, and the physics of electrons, spins, and phonons in condensed matter on the other hand [1][2][3][4][5][6] . The seminal phenomenon considered in this respect was the threedimensional (3D) Anderson localization of light 7,8 -in analogy to Anderson localization of spins 9 -that continues to receive attention to date 10,11 . Other well-known examples are the analogy between weak localization or enhanced backscattering of light 12,13 and of electrons 14 , and the analogy between a complete 3D photonic band gap in a 3D photonic crystal [15][16][17] and the electronic band gap in a semiconductor crystal such as silicon or germanium 18 .
In this work, we explore the propagation of light in a 3D superlattice of coupled resonant cavities inside a 3D photonic band gap. The propagation of light in such a 3D cavity superlattice is analogous to electronic transport in an impurity band in a semiconductor 9,18,19 . Light hops from cavity to cavity throughout the 3D superlattice, which differs fundamentally from the conventional spatially-extended Bloch wave propagation outside the gap. Since the light hops predominantly in a few highsymmetry directions including the Cartesian (x, y, z) directions, we propose the name "Cartesian light" for the unusual propagation of light in the 3D superlattice of coupled cavities in a 3D photonic band gap.
In one dimension (1D), a chain of coupled resonant cavities is a well-known system that is known as a coupled resonator optical waveguide (CROW) 21 . The weak coupling between cavities in a CROW has been demonstrated at optical frequencies 22 . CROWs are widely studied for efficient nonlinear optical frequency conversion and for perfect transmission through bends, and for the one-dimensional (1D) localization of light 21,23 . Two dimensional (2D) arrays of coupled cavities have been studied, notably for unusual discrete diffraction effects 25 , for intricate coupled nanolasers 24 , for topologically-protected propagation 26 , and for transverse localization 27,28 . The coupling between the cavities at optical frequencies has been demonstrated to be significantly larger than the fabrication-induced disorder in the cavity frequencies 29 . In 3D resonator arrays without band gap, topologically-protected propagation was studied 30,31 , as well as the percolation of light through 3D lattices of coupled resonant microspheres 32 , and dynamic localization of light 33 . In 3D photonic band gap crystals in the microwave regime, slow heavy-photon propagation was reported in a 1D array of weakly-coupled cavities [34][35][36] . Numerical calculations of a 2D array of cavities embedded in a 3D woodpile photonic crystal revealed ultraslow and negative group velocities 37 . To the best of our knowledge, 3D superlattices of coupled cavities with resonances in a 3D photonic band gap have not yet been studied before.

II. METHODS
In this paper, we study a 3D cavity superlattice that is embedded in a 3D photonic band gap crystal that has the inverse woodpile structure. This structure has nearly the same symmetry as a diamond crystal of carbon atoms, 39 yet thousandfold magnified, as illustrated in a youtube animation. 40 The inverse woodpile crystal structure consists of two perpendicular 2D arrays of nanopores with radius r in a high-index medium such as silicon, 39 as illustrated in Figure 1(a). Each 2D pore array corresponds to a diamond 110 crystal face. In view of the arrangement of the nanopores, it appears to be convenient to employ a tetragonal unit cell 41,42 instead of the conventional cubic unit cell 18 . The tetragonal unit cell has lattice parameters c (in the x and z directions), and a (in the y direction) in a ratio a/c = √ 2 to ensure a cubic crystal structure. More details, notably on the Brillouin zone, are presented in Appendix A.
Inverse woodpile photonic crystals possess a broad 3D photonic band gap, whose width strongly depends on the radius r of the pores 39,41,42 . For a normalized pore radius r/a = 0.24 -as considered here -a maximum relative bandwidth ∆ω pbg /ω c = 25.3% occurs for ǫ = 12.1 typical of silicon, 41,42 with ∆ω pbg the frequency width of the band gap, and ω c the band gap's center frequency. 3D inverse woodpile crystal nanostructures have been fabricated from a number of different high-refractive index backbones [43][44][45][46][47][48] . In nanophotonic experiments, the potential of silicon inverse woodpiles was demonstrated by the observation of a broad 3D photonic band gap for many angles 49 , as well as a strong spontaneous emission inhibition of embedded quantum dots 50 .
To create a resonant cavity in a inverse woodpile photonic crystal, Ref. 51 proposed a design whereby two proximal perpendicular pores have a smaller radius (r ′ < r) than all other pores, as shown in Figure 1(a). Near the intersection region of the two smaller pores, the light is confined in all three directions to within a mode volume as small as V mode = λ 3 where λ is the free-space wavelength. 51 Supercell band structures revealed up to five resonances within the band gap of the perfect crystal, depending on the defect pore radius r ′51 , that have quadrupolar symmetry. 52 The best confinement occurs for a defect radius r ′ /r = 0.5 that is also considered here. Figure 1(b) shows a 3D superlattice of cavities as is studied here, where each sphere indicates one cavity, as shown in Fig. 1(a). The cavity superlattice has lattice parameters (c x s , a s , c z s ) in the (x, y, z) directions that are integer multiples of the underlying inverse woodpile lattice parameters: c x s = M x c, a s = M y a, c z s = M z c. Here, we study the M x M y M z = 3×3×3 superlattice such that the cavities are repeated every three unit cells with lattice parameters c x s = 3c, a s = 3a, c z s = 3c. Thus, the cavity superlattice is also cubic, similar to the underlying inverse woodpile structure (see section IV for additional discussion).
We have calculated the band structure of the Design of a single cavity in an inverse woodpile photonic band gap crystal shown in a cut-out of a MxMyMz = 3 × 3 × 3 supercell that is surrounded by boxed lines. The high-index backbone is shown in gray. Two proximal smaller defect pores are indicated in green, and the cavity region is highlighted as the bright region at the center. The tetragonal lattice parameters a and c are shown, as well as the x, y, z coordinate system. (b) (x, z) and (x, y) cross sections through a 3D superlattice of resonant cavities, with red circles indicating cavities and dashed rectangles representing unit cells of the underlying inverse woodpile crystal structure (see (a)). The lattice parameters (c x s , as, c z s ) of the superlattice are shown, as well as the x, y, z coordinate system. 3D cavity superlattice using the plane-wave expansion method. 17,18,54 Using the Richardson extrapolation method allows us to estimate the frequencies in the limit of infinite grid resolution 55,56 . Details on the calculations and the convergence are given in Appendix B. All calculations were performed on the "Serendipity" cluster in the MACS group at the MESA + Institute. 57 Even on this powerful computer cluster the calculations took 210 hours. A closer inspection of the five dispersionless cavity bands in Fig. 2(b) reveals that these bands have nonzero bandwidths, indicating that cavity resonances in the M x M y M z = 3 × 3 × 3 superlattice are coupled, as is investigated in this paper. Our results agree well with a simultaneous investigation of a single cavity in an inverse woodpile crystal with finite support, studied by other numerical methods. 52 Notably, Ref. 52 also reports that the first two bands are nearly degenerate. Based on the occurrence of five cavity bands, on degeneracies between bands, and on the field distribution (reported in Ref. 51 ), it has been concluded that the resonances of the inverse woodpile cavity have quadrupolar symmetry and are the optical analogues of d-orbitals in solid-state physics 52 . Therefore, it is naively expected that neighboring cavities couple in diagonal directions.

B. Dispersion bandwidths
For a 1D coupled-resonator optical waveguide (CROW), it is well-known that the coupling coefficient along the waveguide is proportional to the dispersion bandwidth. 21,36,53 A straightforward extension of this notion to a 3D cavity superlattice is to consider the dispersion bandwidth in various crystal directions, since this is straightforward to derive from photonic band structures as in Figure 2. For a given crystal direction characterized by wave vector k, the dispersion bandwidth is defined as in other words, the absolute value of the difference between the maximum and minimum frequencies on a trajectory in reciprocal space between the origin Γ and the edge of the Brillouin zone k BZ in the direction of k. As an example, for the m = 3 coupled-cavity band in Figure 2(b), between Γ and Z the minimum and maximum frequencies are nearly the same (ω = 0.521) hence the bandwidth ∆ω is nearly zero. Between Γ and U the minimum and maximum frequencies differ much more (ω = 0.520 to 0.521) hence the dispersion bandwidth is much greater in the diagonal direction. A polar plot of the dispersion bandwidth ∆ω versus wave vector k in the (k X , k Z ) plane is shown in Figure 3(a). From the band frequencies mentioned above, the dispersion bandwidth is very small in the real-space x and z-directions (corresponding to ΓX, ΓZ, respectively, in reciprocal space). In the diagonal directions that correspond to the ΓU high-symmetry trajectory the dispersion bandwidth is much greater. As seen from a given central cavity in real space, the wave vector k is then directed towards a second nearest neighboring cavity in the diagonal 1/ √ 2.(1, 0, 1) direction (see Appendix A). The polar plot of the dispersion bandwidth for m = 3 therefore looks like a quadrupolar radiation pattern. Based on the 1D CROW-reasoning given above, one tentatively infers that light is transported through the 3D cavity superlattice preferentially in the xz-diagonal (corresponding to ΓU ) directions. Figure 3(b) shows the dispersion bandwidth ∆ω in the (k Y , k U ) plane. The largest bandwidth occurs at about 45 o off the (XU Z) plane, which corresponds to the ΓR high-symmetry direction. The bandwidth in the ΓYdirection is small, from which one tentatively infers that there is little light transport in the y-direction in real space. The dispersion bandwidth for the m = 4, 5 bands is shown in Figures 4(a,b). For the m = 4 band, the dispersion bandwidth is large in the diagonal directions that correspond to the ΓU high-symmetry directions, and it is smaller in the x and z-directions (ΓX,ΓZ, respectively). Compared to the m = 3 band, the dispersion bandwidth for the m = 4 band appears to be less strongly direc-tional. To quantify the directionality, we consider a directionality D ratio between the maximum and the minimum bandwidths D = ∆ω max /∆ω min in the (k X , k Z ) plane, which yields a directionality of about D = 4 that is much lower than D = 15 for the m = 3 band.
For the m = 5 band in Figure 4(c), the polar plot of the dispersion bandwidth looks very much like a quadrupolar emission pattern. The bandwidth is small in the diagonal directions that correspond to the ΓU high-symmetry directions, about 2.5× smaller than for the m = 4 band. The dispersion bandwidth is much smaller in the x and z directions, corresponding to a large directionality D = 26.

C. Coupling coefficients
To understand the coupling between the cavities in a 3D cavity superlattice more fundamentally, we derive the coupling coefficients of light from the dispersion relations using the tight-binding method, see Appendix C for details. Figure 2(b) shows that the m = 3 band is accurately described by the tight-binding model. It appears that only 7 independent coupling coefficients κ are needed in the tight-binding model, namely for the realspace directions x, y, z, xz − diagonal (corresponding to ΓU in reciprocal space), xy − diagonal (corresponding to ΓS), yz − diagonal (corresponding to ΓT ), and xyz − diagonal (corresponding to ΓR). The reasons are as follows: since the inverse woodpile cavity has mirror symmetry with respect to the (y, z) and (x, y) planes, the coupling coefficients in the +x and +z directions are symmetry related to those in the −x and −z directions, respectively, and the coefficients in the xz − diagonal directions are symmetry related to each other. The coupling coefficients in the +y and −y directions are equal by reciprocity, see Appendix D.
The coupling coefficients are given in Table I. For the m = 3 band the coupling coefficients of light are overlaid on the cavity superlattice structure in Figure 5. In the x and z directions, the coupling coefficients are relatively large and positive, in the xz − diagonal directions the coupling coefficients are large and negative, in the y direction the coupling coefficient is about 10× smaller, and in all other directions the coupling coefficients are vanishingly small (typically 100× less).
Remarkably, the simultaneous occurrence of large coupling coefficients with near-vanishing dispersion bandwidths means that 1D CROW-like arguments do not hold for 3D cavity superlattices. In other words, the bandwidth in a particular crystal direction for a 3D cavity superlattice is not necessarily proportional to the coupling coefficient in the same direction in real space. The small difference between the x and z coefficients confirms that the x and z directions are not symmetry related for the inverse-woodpile cavity, as opposed to the perfect inverse-woodpile structure. 59 According to an 1D CROW-like argument, the large The calculated coupling coefficients have in addition to the real part also an imaginary part, that is at least 100× smaller than the real part, that does not have physical significance, and that is not reported here. In addition to the 7 coefficients κ per band (see text), we also provide the β, defined in Eq. (C18) in Appendix C.
Coupling coefficients of light from a central cavity (black circle) to neighboring cavities (other circles) for the m = 3 coupled-cavity band, indicated with arrows. Nonzero coefficients only occur in the (x, y) and (x, z) planes. Blue and red indicate negative and positive coupling coefficients, respectively, as shown by the color bar. The x, y, z coordinate system is shown. This figure has been made using Par-aView 58 . coupling coefficient in the xz −diagonal directions agrees with the observation of a large dispersion bandwidth in the diagonal U direction, see Fig. 3. However, the negative sign disagrees with the fact that at the U point the band frequency is lower than at the Γ point, since the reverse is true for an 1D CROW, see Ref. 21 .
The κ y coupling coefficient in the y direction is small, in agreement with the band frequencies that are almost the same at Γ and Y . It is remarkable that the κ y coupling coefficient is 10× smaller than the coefficients for the x or z directions, while the nearest neighbor distance is only a ( √ 2) greater than in the x or z directions, which would correspond to only a (exp( √ 2) ≈ 4×) smaller coefficient for cavities coupled by evanescent Bloch modes. We surmise that the quadrupolar field pattern of each cavity in the xz plane couples poorly to a neighboring cavity in a neighboring xz plane. Therefore, light mostly hops in 2D (x, z)-layers, which is analogous to 2D electron transport in graphite or graphene layers. 60,61 Since the light propagates very unusually by hopping only in a few discrete directions, we propose the name "Cartesian light" for the propagation of light in a 3D cavity superlattice.
The coupling coefficients of light for the m = 4, 5 bands are shown in Fig. 6. For the m = 4 band, the nonzero coefficients are κ x , and κ z . The coupling coefficients to all other neighboring cavities vanish, including the coefficient κ y in the y direction. In the hopping of the m = 4 band we find the ultimate Cartesian light: light hops only in the x-z directions. For the m = 5 band, the nonzero coefficients are κ x , κ y , κ z , κ xz , κ xy , κ yz , and κ xyz , that is, nonzero coupling coefficients to all neighboring cavities. In the (x, z) plane, there is a mix of positive x and z coupling coefficients, and negative xz −diagonal coupling coefficients, which is similar to the m = 3 band.
We have not analyzed the m = 1, 2 bands of coupledcavity modes, since the band structures do not converge monotonically with increasing spatial resolution, as is elaborated in Appendix B.

D. Propagation in 3D on the superlattice
We now discuss the propagation in real space, and why the dispersion bandwidth of the m = 3 band is much larger in the U direction than in the other directions in the (k X , k Z ) plane in reciprocal space. We first discuss the dispersion bandwidth in the xz − diagonal directions that are symmetry related to each other, and that correspond to the ΓU direction in reciprocal space.
Let us consider the relative phase of the resonating cavities for a wave vector at the U -point (k = k BZ = U ), as shown in Figure 7. Since this coupled resonance eigenmode is a Bloch wave of the 3D superlattice (with a phase front as indicated in Fig. 7) it is clear that this collective oscillation differs fundamentally from waveguiding behavior in a 1D CROW; in other words, the superlattice Bloch modes differ from the ones in a CROW.
In Figure 7, neighboring cavities in the xz − diagonal direction resonate in-phase with each other. Neighboring cavities in the x, or z-direction resonate out-of-phase with each other. Hence, there is a checkerboard pattern of two sublattices of cavities that resonate out-of-phase with each other. The in-phase resonance of neighboring cavities in the xz-diagonal direction is energetically favorable for the negative coupling coefficient in the xzdiagonal directions, as can be understood from Eq. (C19) in Appendix C. The out-of-phase resonance of neighboring cavities in the x, or z-direction is energetically favorable for the positive coupling coefficient in the x, or z-direction. Hence, the checkerboard pattern of out-ofphase resonating cavities is energetically favorable for all coupling coefficients in the (x, z)-plane, and the band frequency is relatively low at the U high-symmetry point.
We now consider the superlattice Bloch mode at the k = Γ point (the origin in reciprocal space). All cavities resonate in-phase with each other, which is energetically favorable for the negative coupling coefficient in the xzdiagonal direction. However, it is not energetically favorable for the positive coupling coefficients in the x, or z-direction. Hence, the band frequency is relatively high at the Γ point, and higher than at the U -point, in agree-z x k u FIG. 7. Relative phase of the resonating cavities for a coupledcavity mode with a wave vector at the U high-symmetry point. As seen from a given cavity, the wave vector k is directed towards an xz diagonally neighboring cavity. The black dashed line indicates a wavefront of the Bloch wave. The blue and red cavities resonate out-of-phase with each other. The couplings are indicated with arrows, which are green if the corresponding cavities resonate with an energetically favorable relative phase, as is the case for all couplings. ment with the observation of large dispersion bandwidth in Figure 2(b).
We now discuss the dispersion bandwidth in any direction in the (k X , k Z ) plane other than the xz-diagonal direction. For example, the relative phase of the resonating cavities at the X high-symmetry point is shown in Fig. 8. Only the checkerboard pattern of out-of-phase resonating cavities at the U high-symmetry point is energetically favorable for all coupling coefficients in the (x, z) plane. For any other point in the Brillouin zone, the relative phase of the resonating cavities is not energetically favorable for all coupling coefficients in the (x, z) plane. Hence, the dispersion bandwidth is small for directions other than the xz-diagonal direction, in agreement with the superlattice band structures in Figure 2  Relative phase of the resonating cavities for a coupled-cavity mode with a wave vector at the X highsymmetry point. As seen from a given cavity, the wave vector k is directed towards a nearest-neighboring cavity. The black dashed line indicates a wavefront of the Bloch wave. The blue and red cavities resonate out-of-phase with each other. The couplings are indicated with arrows, green arrows indicate an energetically favorable relative phase, and red arrows indicate an unfavorable phase. band gap, from coupled-resonator optical waveguiding, and also from light diffusing at the edge of a gap.
(1) A characteristic feature of the 3D superlattice Bloch modes is that they are constructed from modes where the light field is hopping from lattice site to lattice site. In other words, the field pattern has its maxima on the lattice sites (i.e. the cavities) and decays exponentially in between lattice sites, since isolated cavities are tuned inside the photonic band gap where the wave vector is complex. In contrast, the Bloch modes outside the photonic band gap are constructed from purely real modes of propagation; there is no reason for field maxima to be located on preferred positions in the crystal.
(2) A second characteristic feature of the 3D superlattice Bloch modes is that they are genuine modes of propagation centred within the 3D photonic band gap. In this sense, they are distinguished from the modes in photonic crystals with finite support that were recently described in Ref. 63 . In that study, it was found that the finite extent of a photonic band gap crystal leads to the filling of the density of states (DOS) in the band gap by states that are centred outside the band gap, while extending into the band gap due to their substantial band width.
(3) The 3D Cartesian superlattice propagation differs fundamentally from the propagation in lower-dimensional 1D (a CROW) and 2D arrays of cavities. Firstly, in section III D we have already discussed that the superlattice modes differ fundamentally from those of a CROW. In other words, a 3D superlattice does not seem to be a "3D CROW". Secondly, if we perturb the frequency of one of the cavities in a superlattice, a bound state appears instantaneously in 1D and 2D, whereas a threshold frequency difference is required in 3D, see Ref. 20 .
(4) The propagation of light in a 3D cavity superlattice in a photonic band gap differs fundamentally from directional diffusion that was identified for 3D photonic band gap crystals with a certain degree of disorder 38 . In the latter case, the modes of propagation are not waves but diffusive. Moreover, the typical frequencies are at the edge of the band gap, hence outside the gap, as opposed to the cavity superlattice modes that reside within the band gap, see Fig. 2.

B. Crystal structures of the cavity superlattice
If the magnification factors of the superlattice's lattice parameters, compared to the underlying crystal structure's lattice parameters, fulfill M i = (M j , M k ) (i, j, k = x, y, z), the cavity superlattice is not cubic anymorein contrast to the underlying inverse woodpile structure -but has become tetragonal. In the most general case with M x = M y = M z the superlattice has different cavity spacings in each direction (x, y, z); the superlattice has then become orthorhombic. Given that cavities in an inverse woodpile structure are necessarily located along the smaller-pore line defects, we currently doubt whether it is feasible to realize other 3D Bravais superlattices.
We have seen that for the m = 3, 4, 5 bands, the coupling coefficient in the y-direction is smaller than in the (x, z) plane, by typically 10×. To make the hopping of light more 3D, it is necessary to increase the coupling coefficient in the y-direction compared to the coefficients in the (x, z)-directions, which can be achieved by a closer cavity spacing in the y-direction, for example, in a M x M y M z = 3 × 2 × 3 supercell in case of M x = M z = 3 (as studied here) or in general for supercells with M y < (M x , M z ).
Conversely, if it is desired to realize a superlattice with effectively 2D transport of light in (x, z) planes, the superlattice parameter in the y-direction should be made greater than the ones in the (x, z)-directions (M y >> (M x , M z )). In this situation, the 2D transport of light may hold analogies to that of charge carriers in graphite layers 60 or in high-T c superconductors.

C. Disorder
We have studied the dispersion and hopping for superlattices without disorder. Let us briefly comment on the sensitivity of the results to a small degree of disorder, since we performed calculations for several grid resolutions in Appendix B, and a change of the grid resolution implies a slight shift in the geometry. On the one hand, we observed that for all five bands of coupledcavity modes, the center frequency is highly sensitive to the grid resolution and therefore to disorder. This is likely the result of the lightning rod effect of the cavity mode field pattern identified in Ref. 51 , wherein the inverse woodpile cavity resonances have regions of high intensity at sharp corners in the dielectric material. If such sharp corners are slightly distorted, it is quite conceivable that the overlap with the field pattern changes, leading to a change in resonance frequency. On the other hand, we observe for bands m = 3, 4, 5 that the features of the dispersion bands remain the same while the grid resolution is increased. Therefore, we expect the coupling coefficients for the m = 3, 4, 5 bands to be robust to small degrees of disorder.
A 3D cavity superlattice is the photonic analogue of the Anderson model for spins and electrons 9 , albeit in the limit of zero disorder. The 3D cavity superlattice also corresponds to the Hubbard model without interactions 64,65 . We anticipate that the present study may form the basis for further exploration of the physics of the 3D Anderson model for nanophotonic cavity superlattices that will proceed by introducing controlled degrees of disorder in the cavity resonance frequencies.

D. Outlook
A possible application would be a scalable, coherently linked network of NV-based registers, see Ref. 66 .

V. SUMMARY
We have studied for the first time ever the propagation of light in a 3D cavity superlattice within a 3D photonic band gap. Such a 3D cavity superlattice is the photonic analogue of the Anderson model in the limit of zero disorder. The light hops only in a few high-symmetry directions including the Cartesian (x, y, z) directions, therefore we propose the name "Cartesian light ". 3D Cartesian hopping of light in a 3D band gap yields propagation as superlattice Bloch modes that differ fundamentally from the conventional 3D spatially-extended Bloch wave propagation in crystals, from light tunneling through a band gap, from coupled-resonator optical waveguiding, and also from light diffusing at the edge of a gap. The large coupling coefficients in the Cartesian directions occur simultaneously with a near vanishing dispersion bandwidth in this direction. This means that 1D CROW-like do not arguments hold for 3D cavity superlattices. The unusually small dispersion bandwidth in the Cartesian directions is a result of interplay between positive and negative coupling coefficients in the Cartesian and diagonal directions.

VI. ACKNOWLEDGMENTS
We thank Bill Barnes (Exeter, Twente) for coining the term "Cartesian light" and for discussions, and we thank Geert Brocks, Ad Lagendijk, Jan Klärs and Bart van Tiggelen for useful discussions. This research is supported by the 4TU federation, by the FOM/NWO programme "Stirring of light!," the STW/NWO-Perspectief program "Free-form scattering optics", the "Descartes-Huygens" prize of the French Academy of Sciences, and the MESA + Institute for Nanotechnology section Applied Nanophotonics (ANP).
Appendix A: Brillouin zone and tetragonal unit cell First Brillouin zone of the inverse woodpile crystal structure showing the high symmetry points X, Y, Z, S, R, T, U and the origin called Γ. Figure 1(a) shows a tetragonal representation of the M x M y M z = 3 × 3 × 3 supercell of the cubic inverse woodpile crystal structure including a pair of proximal defect nanopores that form a resonant cavity. In terms of the unit vectors of the conventional cubic diamond structure, 18 the tetragonal unit supercell has unit vec- Figure 9 shows the first Brillouin zone of the inverse woodpile cavity superlattice with characteristic highsymmetry points. The main axes are given by X = [π/c x s , 0, 0], y = [0, π/a s , 0], and Z = [0, 0, π/c z s ]. The ΓX and ΓZ directions are notable as they correspond to waves propagating along each set of nanopores, while the diagonal ΓU direction lies in between. Due to the geometry of the cavity (composed of two proximal defect pores), the ΓX and ΓZ directions are not symmetry related, in contrast to the underlying inverse woodpile crystal structure. 59

Appendix B: Calculations and convergence
We have calculated the band structure of 3D cavity superlattices using the plane-wave expansion method, 17,18 where we employed the well-known MIT PhotonicBands (MPB) package 54 . We performed calculations with increasing spatial grid resolutions of 12 × 17 × 12, 24 × 34 × 24, 48 × 68 × 48, and 96 × 136 × 96 per unit cell of the underlying inverse woodpile crystal. Although the 12 × 17 × 12 calculation is a replication of Ref. 51 , the results do not agree perfectly since the subpixel averaging has been updated in Version 1.5 of the MPB code that we use here. We verified that with an older version of MPB our calculations agree exactly with Ref. 51 . To study the convergence of the defect bands with increasing spatial resolution, we show in Figure 10(a,b) the m = 1 and m = 3 defect bands, respectively, for all resolutions considered. The general trend observed in the data -also seen with the other three defect bands -is that the average frequency of the bands decreases with increasing resolution, where the initial decrease is fast, whereas for increasing resolution the decrease slows down, and the Richardson-extrapolated (or converged) band frequencies are reached at a resolution of 96 × 136 × 96.
The initial refinement from 12×17×12 to 24×34×24 is remarkably large in comparison to the bandwidth of each defect band. We speculate that the shifts are related to the "lighting rod"-behavior of the field patterns of a cavity mode, whereby the high fields are concentrated at sharp corners 51 . A change in the grid resolution effectively corresponds to a change in the geometry. At low resolution, the field pattern likely "misses" the sharp fea-tures of the high-index material, whereas with increased resolution the field pattern increasingly fits within the high-index material, causing a decrease of the field energy -in view of the variational principle 17 -as observed in Figure 10.
To quantify the rate of convergence, we calculate the convergence order following Ref. 55 . We assume that the grid spacing h is sufficiently refined for the error E to asymptotically approach zero as with ω(h) the frequency calculated for grid spacing h at a given wave vector, ω exact the exact frequency at the same wave vector, C a constant. The convergence order p is obtained from a sequence of three resolutions as 55 where we used the 12 × 17 × 12, 24 × 34 × 24 and 48 × 68 × 48 resolution results. From data as shown in Figure 10(a,b), we obtain p = 1.85, 1.96, 2.24, 2.16, 2.10 for m = 1, 2, 3, 4, 5, respectively, in close agreement with the convergence order p = 2 for the plane-wave expansion method. 54 We therefore conclude that the band frequencies are converging accurately, as expected. Moreover, the convergence of the band frequencies allows us to use Richardson extrapolation to obtain the band frequency ω h=0 in the limit of infinite grid resolution (h = 0). 55,56 The frequency ω h=0 is estimated as The Richardson-extrapolated frequencies are shown in Fig. 10(a,b) as a function of wave vector for the m = 1 and m = 3 defect bands, respectively. The Richardsonextrapolated frequencies are slightly below the frequencies for the finest grid resolution (96 × 136 × 96) as expected in case of convergence 55 . Nevertheless, unexpected features were found in the dispersion relations as a function of spatial grid resolution. Figures 10(c,d) show the dispersion relations of the m = 1 and m = 3 defect bands that are overlaid for all grid resolutions by subtracting the average band frequencies from the data in Figures 10(a,b). For the m = 1 defect band, we observe that for resolutions 12 × 17 × 12 and 24 × 34 × 24, the band between the T and Y highsymmetry points has a maximum at the 4 th symbol from the left, midway in between T and Y . For 48×68×48 and 96 × 136 × 96, the maximum has moved to the 1 st symbol from the left, or one sixth of the way from T to Y . Thus, there is a qualitative change in the dispersion relations as the grid is refined. Similar observations were made on the m = 2 defect band. Therefore, we do not trust the m = 1, 2 defect bands sufficiently to derive coupling coefficients. In contrast, for the m = 3 defect band, the maximum in frequency is always at the T point, thus the shape is preserved, and the band readily converges. Similar observations as for m = 3 were made on the m = 4, 5 defect bands. Therefore, we trust the convergence of the m = 3, 4, 5 bands sufficiently to warrant the extraction of coupling coefficients.

Appendix C: Photonic tight-binding method
We employ the well-known tight-binding method to model the dispersion of the defect bands 18,21 . In developing the tight-binding approximation, we assume that in the vicinity of each lattice point the full periodic superlattice dielectric function, ǫ, can be approximated by the dielectric function, ǫ Ω , of a single cavity located at the lattice point. This assumption is valid for the cavity superlattice, since the superlattice dielectric function differs only from ǫ Ω at the defect pores in the superlattice that do not form part of the cavity, all of which are at least a lattice constant away. We also assume that the modes of the cavity are well localized; i.e., if E m Ω is a mode of a cavity at the origin, with Ω m the resonance frequency of the single cavity and c the speed of light, then we require that E m Ω (r) be very small when r exceeds a distance of the order of the lattice constant, which we shall refer to as the "range" of E m Ω . Let us briefly verify this assumption: for the defect bands, the mode volume is about V ≈ λ 3 , with λ the free-space wavelength 51 . Hence, the electric field attenuates over a range of about λ/2 in a given direction. We model the cavity mode as an evanescent plane wave that attenuates by a factor of 1/e over this distance, which yields an imaginary part of the wave vector of k ′′ = 2/λ. The nearest neighbor distances are equal to the lattice parameters of the superlattice (c x s , a s , c z s ) that are multiples of the inverse woodpile lattice parameters (c x s , a s , c z s ) = (3c, 3a, 3c) (with c = a/ √ 2), and the reduced frequency ω ′ = a/λ of a cavity resonance is equal to ω ′ = 0.52...0.54 (for m = 3, 4, 5). For the directions x or z, this yields a product of the imaginary part of the wave vector k ′′ and the nearest neighbor distance ∆r of Thus the electric field intensity from one cavity has attenuated to as little as exp(−2.2) ≈ 0.11 at a nearest neighboring cavity in the directions x or z, and for the y-direction the decay is even greater, thereby readily fulfilling the requirements of the tight-binding method.
In the tight-binding method, we write the superlattice dielectric function ǫ as ǫ = ǫ Ω + ∆ǫ(r), where ∆ǫ(r) contains all corrections to the cavity dielectric function required to produce the full periodic dielectric function of the superlattice. Since the product ∆ǫ(r)E m Ω (r), though nonzero, is exceedingly small, we might expect the solution to the full superlattice Maxwell equations to be quite close to the cavity wave function E m Ω (r) or to wave functions with which E m Ω (r) is degenerate. Based on this expectation, one seeks an E(r) that can be expanded in a relatively small number of localized cavity wave functions: The wave vector k describes the relative phase of the resonating cavities and ensures that k satisfies the Bloch theorem for the Brillouin zone of the superlattice of cavities.
If we multiply the superlattice Maxwell equations ǫ Ω (r)E n * Ω (r) · E m Ω (r)dr = δ mn , we arrive at an eigenvalue equation that determines the coefficients b m (k) and the Bloch frequencies (ω(k)/c) 2 : The first term on the right of Eq. (C9) contains integrals of the form ∞ −∞ drE n * Ω (r) · ǫ Ω (r)E m Ω (r − R).
We interpret our assumption of well-localized cavity modes to mean that Eq. (C10) is small compared to unity. We assume that the integrals in the third term on the right of Eq. (C9) are small, since they also contain the product of two cavity wave functions centered at different sites. Finally, we assume that the second term on the right of Eq. (C9) is small because we expect the cavity wave functions to become small at distances large enough for the periodic dielectric function to deviate appreciably from the cavity one. Consequently, the right-hand side of (C10) (and therefore [(Ω n /c) 2 −(ω(k)/c) 2 ]b n ) is always small. This is possible if (Ω n /c) 2 − (ω(k)/c) 2 is small whenever b n is not (and vice versa). Thus (ω(k)/c) 2 must be close to a cavity mode, say (Ω 0 /c) 2 , and all the b n except those going with that mode and modes degenerate with (or close to) it in frequency must be small: The coupling coefficients are extracted from the dispersion by means of a least-squares fit of the right-hand side of Eq. (C20) to the dispersion throughout the whole first Brillouin zone. We represent the first Brillouin zone by a grid of 12×12×12 cubes. We fit over the k vectors in the middle of the cubes. The results of the tight-binding calculations are discussed in Section III. We note that for all nonzero coupling coefficients, the imaginary part is of the order of at most 10 −6 , which is 100× to 1000× less than the real part (of the order of ∼ 10 −4 for m = 3, 5 and ∼ 10 −3 for m = 4). Therefore, we are confident that the coupling coefficients are physically significant. As expected, β is negligibly small, with a real part of the order of ∼ 10 −10 .