Cubic 3D Chern photonic insulators with orientable large Chern vectors

Abstract

Time Reversal Symmetry (TRS) broken topological phases provide gapless surface states protected by topology, regardless of additional internal symmetries, spin or valley degrees of freedom. Despite the numerous demonstrations of 2D topological phases, few examples of 3D topological systems with TRS breaking exist. In this article, we devise a general strategy to design 3D Chern insulating (3D CI) cubic photonic crystals in a weakly TRS broken environment with orientable and arbitrarily large Chern vectors. The designs display topologically protected chiral and unidirectional surface states with disjoint equifrequency loops. The resulting crystals present the following characteristics: First, by increasing the Chern number, multiple surface states channels can be supported. Second, the Chern vector can be oriented along any direction simply changing the magnetization axis, opening up larger 3D CI/3D CI interfacing possibilities as compared to 2D. Third, by lowering the TRS breaking requirements, the system is ideal for realistic photonic applications where the magnetic response is weak.

Introduction

Inspired by the discoveries of topological phenomena in solid state systems, the study of topology in the propagation of light in photonic crystals has been the subject of much recent attention^{1,2,3,4,5,6,7,8,9,10,11,12,13}. Among all topological states of matter, time-reversal symmetry (TRS) broken topological materials, such as Chern insulators (CI)^2,3 and lasers¹⁴, have been a particular focus due to their topologically protected unidirectional edge states with non-reciprocal propagation properties. In these systems, scattering processes from one boundary state into another are strongly suppressed, due to decoupling of counter-propagating 1D chiral edge channels^15,16.

Seminal works in 2D photonics have demonstrated that transverse magnetic (TM) modes in gyro-magnetic photonic crystals could mimic the Chern insulating state for light^4,5. Due to a non-zero value of the topologically invariant Chern number, these systems were shown to sustain topologically protected one-way edge states with negligible dissipation and absence of back-scattering, even in presence of impurities and lattice defects which break translational symmetry.

As originally pointed out by Ref. ¹⁷, extending these ideas to 3D is in principle possible, under some more stringent conditions. As an example, preserving the translational symmetry of the lattice, a Chern insulating phase is predicted in 3D to host chiral anomalous surface states (SS) on its boundary^{18,19,20,21,22,23}. In contrast to 2D, a 3D Chern insulator (3D CI) is a topological phase that can be characterized by three first Chern invariants–or a Chern vector C = (C_x, C_y, C_z) - defined on lower dimensional surfaces^24,25,26: such a state of matter can support chiral surface states propagating on the planes with Miller indices indicated by the Chern vector.

Previous theoretical efforts in photonics²⁷ have engineered a TRS broken insulator with a single nonzero Chern number of unit value in a 3D uniaxial structure and employing a large magnetic field. However, the Chern vector was fixed to have a single component along a preferred axis selected by the fabrication, a situation similar to that of a stack of 2D Chern layers. Moreover, this design required strong TRS breaking, which is usually an arduous challenge in photonic crystals. Finally, the Chern number value was limited to unity.

At the same time, large Chern numbers have been observed in the band-gaps of 2D square photonic crystals^28,29. Topological photonic systems with large Chern numbers can sustain multiple spatially separated edge states²⁹. These edge states allow a plethora of applications, including unidirectional multimode waveguides where information can be multiplexed through the different edge states allowing for photonic on-chip communications with higher channel capacity²⁹. Nevertheless, the value of the Chern number in these 2D systems was not determined by design, but a consequence of the particular system under study. In this sense, finding an engineering strategy to create photonic crystals with any given Chern number would be highly desirable and it is still a challenge that remains open.

In this work, we propose a method to design cubic 3D topological photonic crystals where Chern vectors of any magnitude, sign or direction can be implemented at will. Our method is based on the merging and annihilation of Weyl points through multi-fold supercell modulations in three dimensions³⁰. As a result, we obtain a 3D CI phase with the following characteristics: First, since the Chern number is additive with respect to band-folding over large supercells, the system can support arbitrarily large Chern numbers and thus the coexistence of multi-channel unidirectional surface states. Second, owing to the cubic symmetry of the underlying modulated structure, the system can exhibit any combination of nonzero elements of a Chern vector, allowing for more 3D CI/CI interfacing combinations as compared to 2D. Third, through the combined use of multi-fold supercells and reduced manipulation of Weyl points, the 3D CI phase can be realized under weak magnetization conditions, suited for realistic photonic applications. As a final step, we verify the emergence of chiral SS at interfaces between regions with differing Chern vectors, confirming the existence and spatial separation of unidirectional chiral partners. The outline of the paper is as follows: in the Results section, we provide full topological characterization of the bulk and the boundary of our photonic 3D CI. For the bulk, we show that the cubic system can support any nonzero element of the first Chern vector. The direction of the Chern vector can be tuned by changing the orientation of the applied external static magnetic field. We also prove a strategy for obtaining the 3D CI under a minimal magnetization and a method to design 3D CIs with specifically desired large Chern numbers. For the boundary, we demonstrate the emergence of unidirectional gapless SS and we analyze their 3D anomalous chiral properties. In the Methods section, we describe the numerical implementation of the technique by which we characterize the gap topology based on a 3D generalization of the photonic Wilson loop approach and provide a group theoretical analysis of the mechanism through which we create and manipulate the Weyl points and open up a topological gap by the use of supercell modulations. Further details on our design are given in the Supplementary Information.

Results

The starting point of our design is a photonic crystal with a unit cell containing four dielectric rods directed along the main diagonals of a cubic crystal (scalable lattice parameter ∣a∣). The rods meet at the origin of the unit cell, and the structure is invariant under the operations of the centrosymmetric and non-symmorphic space group (SG) Pn$\bar{3}$m (No. 224)³¹. Since we will later consider modulations of this structure, it is convenient to simulate the dielectric rods by assembling dielectric spheres with radius r = r₀ along (x, y, z)₀/∣a∣ = (t, t, t), (t, 1 − t, 1 − t), (1 − t, t, 1 − t), (1 − t, 1 − t, t) with 0 < t < 1/2, i.e. employing a spherical covering approximation³². The resulting design is shown in Fig. 1a. To obtain a TRS-preserving system, the dielectric material is described by a diagonal (isotropic) permittivity tensor, ${\varepsilon }_{TRS}=\varepsilon {{\mathbb{1}}}_{3}$, where ${{\mathbb{1}}}_{3}=\hat{{{{{{{{\bf{x}}}}}}}}}\otimes \hat{{{{{{{{\bf{x}}}}}}}}}+\hat{{{{{{{{\bf{y}}}}}}}}}\otimes \hat{{{{{{{{\bf{y}}}}}}}}}+\hat{{{{{{{{\bf{z}}}}}}}}}\otimes \hat{{{{{{{{\bf{z}}}}}}}}}$ and ε ∈ ℝ (no losses), and by unit magnetic permeability $\mu ={{\mathbb{1}}}_{3}$. To simulate the optical response of the system, we employ the MIT Photonic Bands (MPB) software package³³. As shown in Fig. 1a, with TRS, the photonic band-structure presents a three-fold degeneracy between the three lowest energy bands at the high symmetry point ${{{{{{{\bf{R}}}}}}}}=\frac{2\pi }{| a| }(1/2,1/2,1/2)$; note that, in the displayed energy window, the two lowest bands are fully degenerate. Everywhere else in the Brillouin zone (BZ), there is a gap between the second and third band. The dispersion reflects the three-fold rotational symmetry of the crystal, and so is invariant under cyclic permutations of the three k_i (i = x, y, z) directions. In order to keep the notation consistent throughout the paper and to capture all the variety of symmetry designs, we label the high symmetry points in the BZ according to the convention described in the Supplementary Note (SN) 5, i.e. in a Cartesian orthorhombic convention.

Fig. 1: Photonic 3D CI by cubic supercell modulation at N = N_W = 2.

Each (a, b, d, e) panel shows the crystal unit cell, the irreducible Brillouin zone (IBZ) and the band-structure (BS). Frequencies f are given in reduced units, ∣a∣ being the scale invariant lattice parameter and c the speed of light. Sectors (c) and (f) contain the topological characterization via photonic Wilson loops (WL), in the Weyl semimetallic (WS) phase and in the 3D CI phase, respectively. a Photonic crystal constructed from cylinders of radius r₀ = 0.15 and dielectric constant ε = 16 at TRS. The three lowest photonic modes display a three-fold degeneracy at R and the two lowest bands are fully degenerate in the displayed energy window. b TRS breaking implemented via a gyro-electric response with ${\eta }_{z}^{{N}_{W} = 2}=16$: the bias field is adjusted in order to split the Weyl points of approximately half the BZ, i.e. at ${k}_{z}^{\pm }=\pm \! \frac{\pi }{2| a| }$, along the ${{{{{{{\bf{S}}}}}}}}{{{{{{{\bf{R}}}}}}}}{{{{{{{{\bf{S}}}}}}}}}^{\prime}$ line where ${{{{{{{{\bf{S}}}}}}}}}^{\prime}={{{{{{{\bf{S}}}}}}}}-{{{{{{{{\bf{b}}}}}}}}}_{z}$. c Electromagnetic section Chern number calculated on 2D k_z planes normal to the magnetization (upper panel) and transverse flow of the θ_y eigenvalue of the WL matrix summed over the entire subset of ν_occ bands lying below the local gap at k_z (lower panel). The discontinuity of the section Chern number at the wavevector of each Weyl point ${k}_{z}^{\pm }$ is used a measure its topological charge (q_± = 1). d Artificial folding of the bands on a N = 2 cubic supercell: Weyl points superimpose at Z in the new BZ. e Coupling and annihilation of Weyl points through a N = 2 supercell modulation with parameter r_m = r₀/20. The amplified modulation is graphically visualized via a colorbar with r_max = r₀ + r_m and r_min = r₀ − r_m. A topological direct gap (f_g) at Z is opened, with gap-to-midgap (f_g/f_m) ratio of 1.86%. The size of gap can be appropriately tuned by choosing the value of the modulation, as in the inset. f The section Chern C_z number displays constant unit value everywhere in the BZ, establishing the system to be in the 3D CI phase. Source data are provided as a Source Data file.

Following a strategy introduced in^31,34,35 we break TRS to split the three-fold degeneracy at R into two Weyl points³¹. This can be achieved by either applying an external magnetic field bias to a gyro-electric crystal^2,3 as we proceed here, or by exploiting the internal remnant magnetization of ferri-magnetic materials⁵. TRS breaking is implemented by introducing off-diagonal imaginary elements in the permittivity tensor: for the specific case of an applied magnetic field ${{{{{{{\bf{B}}}}}}}}={B}_{z}\hat{{{{{{{{\bf{z}}}}}}}}}$, the induced gyro-electric tensor is:

$${\varepsilon }_{{\eta }_{z}}=\left(\begin{array}{ccc}{\varepsilon }_{\perp }&i{\eta }_{z}&0\\ -i{\eta }_{z}&{\varepsilon }_{\perp }&0\\ 0&0&\varepsilon \end{array}\right),$$

(1)

where η_z = η_z(B_z) the bias-dependent gyro-electric parameter and ${\varepsilon }_{\perp }=\sqrt{{\varepsilon }^{2}+{\eta }_{z}^{2}}$. The gyro-electric tensor corresponding to magnetic fields in other directions can be obtained by orthogonal rotations of Eq. (1).

Under these conditions, the three-fold degeneracy splits into a pair of Weyl points (or a Weyl dipole). In order to analyze the formation and splitting of the Weyl points under TRS breaking and to predict the direction of their displacement in the BZ, we develop a k ⋅ p model around the three-fold degeneracy at R. The model, based on the group theoretical method of invariants (see Methods), allows us to conclude that Weyl points appear at inversion-symmetric positions with respect to R along the k_z direction, with a separation that can be adjusted by choosing the bias field B_z appropriately. Our MPB simulations, presented in Fig. 1b, confirm this predictions accurately. More generally, for a magnetization applied along any of the main coordinate axes x_i, the Weyl dipole is oriented along the line joining R to ${{{{{{{{\bf{R}}}}}}}}}^{\prime}={{{{{{{\bf{R}}}}}}}}-{{{{{{{{\bf{b}}}}}}}}}_{i}$, where b_i is the corresponding primitive reciprocal lattice vector. For a detailed comparison of the analytic model and the numerical simulations see Methods.

Next, we calculate the chiral topological charge q_± of the Weyl points in the k ⋅ p model using the Z2Pack numerical tool^36,37, concluding that the Weyl points have opposite valued unit charges (q_± = ±1).

We confirm these predictions by computing the topological charges directly from the MPB eigenstate solutions. To do so, we implement a numerical approach based on the analysis of the winding properties of photonic hybrid Wannier energy centers (WEC). WEC are accurately defined in the Methods section. There, we establish a mapping from electronic Wannier charge centers (WCC)^38,39 to photonics and we perform a generalization of the photonic Wilson loop approach of Ref. ⁴⁰ (initially implemented for 2D scalar waves) applicable to fully 3D electromagnetic (EM) vector fields. The results of this analysis are summarized in Fig. 1c: the top panel shows the electromagnetic Chern number C_z of the two lowest bands calculated on 2D planes orthogonal to the magnetization axis. We observe a sharp discontinuity ΔC_z = ±1 at the wavevector of each Weyl point. This is similarly reflected in the winding of the Wilson loop eigenvalues on two selected planes, as shown in the bottom panel. From the discontinuity in the section Chern number at each Weyl point, we deduce the associated topological charge⁴¹, confirming that q_± = ±1, as predicted by the k ⋅ p model.

In order to identify a geometrical perturbation able to open a band-gap, we analyze the coupling of the Weyl points from a group theoretical perspective, performing a generalization of the method of invariants now suited to capture translation breaking perturbations (see Methods). We conclude that the only deformation of the geometrical structure leading to their annihilation and to the opening of a topological gap are lattice commensurate supercell modulations. In particular, we find that it is possible to independently activate the supercell modulation along the x_i Cartesian directions and couple Weyl points generated by the corresponding magnetic field B_i. Note that, from our WEC analysis, we see that each BZ plane between the Weyl nodes carries Chern number ∣C_z∣ = 1, while each plane outside the Weyl nodes carries Chern number ∣C_z∣ = 0. When we couple the Weyl nodes by a lattice commensurate modulation, we backfold the BZ into a region commensurate with the Weyl node separation vector, which is a reciprocal lattice vector in the folded BZ. The additivity of the Chern number then ensures that every plane in the reduced BZ carries a nonzero Chern number, resulting in a 3D CI when the gap is opened³⁰. This expresses the fact that the Chern number density of our 3D system does not change as a function of the (TR-even) supercell modulation; it simply goes from being unquantized in the original system (necessitating the existence of Weyl points), to being a quantized multiple of a reciprocal lattice vector in the modulated system.

With this starting setup, in order to obtain 3D Chern insulating phases, we will follow a general three-step strategy:

• 1. First, using the external magnetic field we move the Weyl points at fractional distances of the Brillouin zone (BZ), i.e. at positions ${{{{{{{{\bf{K}}}}}}}}}_{1,2}={{{{{{{\bf{R}}}}}}}}\pm \frac{{{{{{{{{\bf{X}}}}}}}}}_{i}}{{N}_{W}}$ where ${N}_{W}\in {\mathbb{N}}$ and N_W > 1. In this way, in a further step, we will be able to couple and gap the Weyl points with a commensurate modulation of a supercell structure. Notice that larger N_W are associated to smaller splittings.

• 2. Secondly, we fold the BZ by creating multi-fold (N > 1) supercells; this is achieved by replicating the original unit cell either in a cubic supercell of dimensions (N, N, N) or in a uniaxial supercell of size N directed along the magnetic field direction. This step of the procedure will merge the Weyl points, originally at K_1,2 in the natural BZ, to the same k point in the new reduced BZ, forming a four-fold degeneracy. In the SN 2, we show that fine tuning and perfect band folding are not strictly necessary for opening a gap at the Weyl points. This endows our system with a robustness and tolerance against reciprocal lattice vector mismatches.

• 3. As a third and last step, we couple and gap the opposite-charge Weyl points by spatially modulating the crystal geometry with a periodicity commensurate to the designed supercell. More specifically, we vary the radius of the cylinders through the entire supercell: numerically, this is achieved by locally changing the radius of the spheres in the covering approximation, from their original r₀ radius to the new local one r(x, y, z). Coherently to the choice made in the previous point 2, this is either done with a cubic modulation of the type: ${{\Delta }}r(x,y,z)=r(x,y,z)-{r}_{0}={r}_{m}[\cos (2\pi x/N| a| )+\cos (2\pi y/N| a| )+\cos (2\pi z/N| a| )]$ when all the Cartesian components of the modulation are turned on or with a uniaxial modulation, where only the component oriented along the magnetic field is activated, e.g. ${{\Delta }}r(x,y,z)={r}_{m}\cos (2\pi {x}_{i}/N| a| )$ for a field with B_i ≠ 0 field. More details are given in the SN 6.

Depending on the values of the parameters N_W and N, it is possible to design different tailored 3D CI phases, in particular: a cubic 3D CI with orientable Chern vectors, a 3D CI in a reduced magnetization environment and a 3D CI with tunable larger Chern numbers.

We stress that the argument of gap opening by folding and supercell modulation is very general, and can be applied as long as the constraints of commensurability between the Weyl diplacement and the supercell size are satisfied. Therefore, any other crystal structure exhibiting a pair of Weyl points could be perfectly suited to their annihilation via the mechanism proposed.

Cubic 3D CI

Our first objective is to design a cubic 3D CI with orientable Chern vectors. As we will show, this can be achieved using a cubic supercell modulation with N = N_W > 1. In order to keep the MPB simulations computationally affordable we consider the simplest case of N = N_W = 2, which requires to separate the Weyl points to half the BZ as in Fig. 1b. The effect of band folding in such a system is visualized in Fig. 1d: for a field oriented as B_z, the two Weyl points superimpose to form an artificial four-fold degeneracy at X₃ ≡ Z. More generally, on a N = N_W cubic supercell and from simple folding considerations, we expect the opposite-charge Weyl points to merge at X_i when N = N_W is even, at R_i − X_i when N = N_W is odd, where ${{{{{{{{\bf{X}}}}}}}}}_{i}=\frac{{{{{{{{{\bf{b}}}}}}}}}_{i}}{2}\equiv {{{{{{{\bf{X}}}}}}}},{{{{{{{\bf{Y}}}}}}}},{{{{{{{\bf{Z}}}}}}}}$. From this starting point, in order to realize a cubic 3D CI phase with orientable Chern vectors, all the three Cartesian components of the cubic commensurate modulation need to be simultaneously turned on. The resulting photonic band structure of the N_W = N = 2 supercell modulated structure is shown in Fig. 1e: as it can be seen, the Weyl points annihilate and open up a gap. To numerically verify the topological properties of this bulk gap in our design we compute photonic Wilson loops and analyze their winding in the BZ (see Methods). Our results, summarized in Fig. 1f, determine that the obtained insulating phase acquires a nonzero Chern number along every plane perpendicular to the magnetization axis as predicted by the k ⋅ p model. Therefore, by simply changing the orientation of the magnetization axis, it is possible to select each Cartesian component in a first Chern class vector (C_x, C_y, C_z), due to the the cubic nature of the underlying system and modulation. Note that the existence of three weak indices in 3D allows for more interfacing possibilities as compared to the 2D case, where only the trivial/TI and the opposite (or different) Chern number combinations are realizable, as discussed later.

3D CI at reduced magnetization

In our previous example we required the Weyl points to be displaced to the half of the BZ. Achieving such a condition requires large TRS breaking parameters. In our simulations, fulfilling this requirement implied using a magnetization bias corresponding to ${\eta }^{{N}_{W} = 2}=16$. Note that, to date, large gyrotropic parameters have been experimentally achieved in photonic crystals only in the microwave frequency regime via ferri-magnetic materials^5,42 and that the gyrotropic response of most currently known dielectric materials is weak. Therefore, in this section, we suggest a way to hugely reduce the magnetization requirements for obtaining CI phase by employing multi-fold supercells and by increasing the intensity of the modulation. Instead of displacing the Weyl points to half the BZ and applying a supercell modulation over two original unit cells, we now move the Weyl points to a smaller fractional distance of the BZ and apply a supercell modulation over a larger number of original supercells to merge and gap the Weyl points appropriately. The resulting 3D CI phase still displays the same Chern number as in the maximally TRS broken system with N = N_W = 2, but it occurs in a largely reduced magnetic field environment due to the smaller k-space displacement of the Weyl points.

For example, making N = N_W = 3, which corresponds to a dipole separation of one third of the BZ and spatially modulating the structure over 3 original unit cells, one can get the same topological phase as in the N = N_W = 2 example. In order to keep the calculations computationally feasible, we simulate a uniaxial system of size (1, 1, N). Nevertheless, the concept is readily generalizable to cubic supercells. Under this construction, the CI phase is achieved at ${\eta }_{z}^{{N}_{W} = 3}=7.8$. As it can be naturally expected, the resulting 3D CI suffers a moderate reduction in the band-gap. However, as shown in Table 1 and Fig. 2 where we compare the designs at N = N_W = 2 and N = N_W = 3, the later design presents a good compromise between gap size and TRS-breaking amplitude, considering the large advantage coming from a large drop in the required magnetization bias (from ${\eta }^{{N}_{W} = 2}=16$ down to ${\eta }^{{N}_{W} = 3}=7.8$). This could be of particular interest for photonic applications where the magnetic response is weak. The gyrotropic parameter value can be further decreased, as shown in the SN 3, by modulating over even larger supercells N_W = 5, 6, 7 and by optimizing the modulation intensity in order to partially compensate for the band-gap decrease. However we cannot indefinitely iterate this procedure down to zero bias since a compromise on the band-gap is always unavoidable: indeed, in the limit of very large N, there is no splitting of Weyl points and thus no TRS broken gap can clearly be opened, ${\lim }_{N\to \infty }{f}_{g}=0$.

Table 1 Reduction of the magnetic bias η and of the topological gap-to-midgap ratio f_g/f_m(%), still achieving the same Chern number C. Source data are provided as a Source Data file.

Fig. 2: 3D CI in a reduced magnetization environment. Uniaxial supercells (1, 1, N) with modulation parameter r_m = r₀/20 and WLs on selected planes for k_z/2π = 0.3.

a, b 3D CI with C_z = 1 at large magnetization ${\eta }^{{N}_{W} = 2}=16$, corresponding to the maximum Weyl dipole separation and band-gap. c, d 3D CI in a reduced magnetic environment ${\eta }^{{N}_{W} = 3}=7.8$. The band-gap only suffers a moderate contraction, yet the Chern vector and the topological properties are preserved. Source data are provided as a Source Data file.

Figure 2 and Table 1 show the drop in the magnetization as compared to the decrease in the topological gap computed for uniaxial supercells at N = 2 and N = 3 with the same supercell modulation parameter r_m.

3D CI with larger Chern numbers

Lastly, we show that our design strategy can also be used to design photonic TIs wih larger Chern numbers. This can be achieved by modulating over even multi-fold supercells with N = 2n > 2, $n\in {\mathbb{N}}$, while keeping N_W = 2. The use of larger supercells permits folding of the BZ multiple times. In the band-folding process the Chern number contribution in each folded region of the BZ adds up. We thus expect the gap resulting for such a modulated system to achieve larger Chern numbers C_i according to the following relation: C_i = n. To prove this, we build uniaxial supercells of size (1, 1, 2n) with n = 1, 2, 3, 4. These crystalline supercells are magnetized along the ${{{\hat{{{{{\bf{z}}}}}}}}}$ direction, creating Weyl points at half of the original BZ (N_W = 2). After folding, we find that the Weyl points are superimposed at R − Z ≡ S if n is even and at R if n is odd. As a final step we activate the modulation along the z direction. In Fig. 3 we then calculate the Chern number of the band gaps in these systems using photonic Wilson loops (WL) by analyzing their winding in the BZ. The modulated supercells with n = 1, 2, 3, 4 under the appropriate TRS breaking acquire, as predicted, Chern numbers C_z = 1, 2, 3, 4 respectively. We restrict our calculations of uniaxial systems due to computational limitations, nevertheless, the argument for the Chern number growth holds similarly in the cubic case, when all the three components of the modulation are turned on.

Fig. 3: Generating larger Chern numbers.

Annihilation of Weyl points at N_W = 2 over multi-fold supercells of even N. a Unit Chern number 3D CI. b, c, d Increasing Chern numbers C according to the relation C = N/2. Calculations performed on 1D supercells (1, 1, N) with modulation parameter r_m = r₀/20 and WL on selected planes for k_z/2π = 0.3. Source data are provided as a Source Data file.

Chiral surface states

Finally, to characterize the bulk-boundary correspondence in the designed systems, we now analyze the emergence of SS at the interface between the cubic 3D CI and a trivially gapped photonic crystal. Other CI/CI interfacing possibilities are discussed in the SN 5. Finding a proper insulating interface is an important requirement to prevent propagation of edge modes in free space due to modes living in the light cone. Furthermore, we also must avoid the formation of dangling defect states due to lattice mismatches. This is usually a quite difficult task in 3D, due to limited available band-gap geometries as compared 2D (details on the trivial interface in the SN 8). To keep the simulations numerically affordable, we stick to the case of a cubic supercell with N = N_W = 2 and analyze a topological slab with normal vector oriented along $\hat{{{{{{{{\bf{x}}}}}}}}}$, in presence of a B_z field. From the bulk-boundary considerations, we expect unidirectional chiral SS to appear on the planes parallel to the magnetization (i.e., with normal vectors perpendidcular to the magnetization direction). Surface states are considered unidirectional in the following sense: The component of the group velocity (or Poynting vector) normal to the magnetization direction has a well defined sign i.e. surface states cannot back-scatter along this specific direction. This component will be later denoted as conserved component.

With this setup, we characterize the hallmarks of chiral SS propagation using a combined real-reciprocal space analysis. Fig. 4a shows the band-structure for the (100)-surface, confirming the emergence of chiral SS connecting the lower and upper bands and fully crossing the band-gap. To better visualize the SS energy dispersion, in Fig. 4b we consider a 3D surface plot, out of which we take the midgap equifrequency cut shown in Fig. 4c. We observe the emergence of 3D Chern Fermi loops which are the natural evolution of the Fermi arcs of the photonic Weyl semimetallic phase. In the SN 11, we show how the SS of the WS phase evolve into the 3D Chern SS as an consequence of Weyl points annihilation. As we will show now, these Fermi loops can be separated in real space, i.e. we can associate those with positive group velocity component normal to the magnetization to a surface of the slab and those with negative one to the other surface. We will therefore consider them as disjoint. In order to establish the relation between counter-propagating modes with respect to the direction orthogonal to the magnetization axis $\hat{{{{{{{{\bf{z}}}}}}}}}$ and the interface normal $\hat{{{{{{{{\bf{x}}}}}}}}}$, i.e. $\hat{{{{{{{{\bf{y}}}}}}}}}=\hat{{{{{{{{\bf{z}}}}}}}}}\times \hat{{{{{{{{\bf{x}}}}}}}}}$, we analyze the propagation of individual edge channels at fixed k_z. As indicated by red/blue colors in Fig. 4(b, c), modes propagating with positive transverse group velocity v_y > 0 appear on one side of the topological slab, their flow being compensated by counter-propagating v_y < 0 partners located on the other surface of the slab. We define as chiral partners, the pair of surface states living on the opposite sides of the slab, moving with opposite component of the group velocity which is normal to magnetization axis. This feature is visualized in Fig. 4(d, f) where we select a pair of chiral partners for explanatory purposes and display their electric field profile in real space on cross sectional view of the crystal slab. We conclude that each disjoint piece of the SS energy sheet in Fig. 4c corresponds to v_y > 0 and v_y < 0: this spatial separation of chiral partners, provided by the bulk, is the protection mechanism which prevents the back-scattering of one state into the other. Because of this, the presence of touching points in the SS dispersion between different chiral partners (red and blue lines in Fig. 4c), is purely accidental. As so, these crossings occur between states that reside on opposite sides of the slab and are physically separated in real space by the bulk. Therefore they cannot gap out, up to exponentially small finite size effects, and are protected by the spatial separation separation of chiral partners on opposite surfaces.

Fig. 4: 3D CI surface states.

Projected band-structure on the surface BZ with Miller index (100), for the interface between N_triv = 8 cells of trivial insulator and N_TI = 8 cells of TI, under a ${{{{{{{\bf{B}}}}}}}}={B}_{z}\hat{{{{{{{{\bf{z}}}}}}}}}$ field and with unit Chern number. Edge states cross the topological gap, highlighted in the shaded region. The extended BZ is fully displayed from −Z to Z, since it lacks of TRS due to application of a $\hat{{{{{{{{\bf{z}}}}}}}}}$ directed bias field. b 2D dispersion of topological surface states on the projected BZ, in the energy range of the topological band-gap. c Equifrequency loops for a cut taken at midgap: the arrows on the plot indicate the direction of the group velocity v_g = ∇_kf∣a∣/c = (v_y, v_z). Blue and red colors correspond to chiral partners with opposite v_y and same optical chirality $\bar{c}$ which are located on opposite left/right sides of the interface (L, R). d Edge states dispersion along a direction normal to both the interface and the external magnetic field; highlighted in circles is a pair of counterpropagating m − th edge channels with ±k_y: the spatial profile of their total electric field is shown in two following panels. e Counterpropagating chiral partners located on the two opposite left/right surfaces of the sample (L, R): wavefront propagation in the k direction indicated by White arrows on a xy cross section. f Spatial profile of the total electric field on a yz planar cut: green and blue arrows indicate respectively the group velocity v_g and the total Poynting vector S, the last being obtained from a spatial integral of the quantity Re(E × H^*), calculated directly from the EM field profiles. Differently from what expected from a conventional QHE edge channel, the Poynting vector displays a component parallel to the external magnetic field (S_z ≠ 0), suggesting both QHE and CME features and the possibility of energy flow along the magnetization axis: we excluded such a possibility in the SN 9, by summing up the contributions of all the edge channels constituting the surface state. Source data are provided as a Source Data file.

Fig. 5: Topological characterization via Wilson loops.

Wilson loops for the system in the Weyl semimetal phase (in blue) and in the cubic 3D C_zI phase (in red). Each θ eigenvalue of the WL is related to the WEC in terms of the i-th lattice parameter ∣a_i∣ by 2πWEC_n,i ≡ ∣a_i∣θ_n,i, with ∣a_i∣ = ∣a∣ and ∣a_i∣ = N∣a∣ for the WS and the 3D CI respectively. a, b Individual photonic WECs: flow in the k_x direction on a selected fixed k_z plane, for each of the bands lying below the local gap. c, d Net photonic WEC flow, obtained by summing up the contribution of the individual WECs shown in the above panels: the net winding appears now clearly. Source data are provided as a Source Data file.

Furthermore, in order to investigate the possibility of energy propagation along the magnetization axis, we analyze the Poynting vector associated to our SS. We find that, even if individual edge channels display nonzero propagation along the magnetization direction, e.g. as in Fig. 4f, integrating the total contribution of entire SS yields no net energy transport along the bias field, as expected due to equilibrium conditions (see SN 9). Interestingly, analyzing the polarization state of each edge channel, we observe a well defined sign of the spatially-averaged optical chirality $\bar{c} \, > \, 0$^{43,44,45,46,47,48} through the entire SS (see SN 10).

Discussions

In this work, we developed a strategy to induce annihilation of Weyl points through cubic and multi-fold supercell modulations, allowing us to achieve a photonic 3D CI phase with the following characteristics:

First, arbitrarily large Chern numbers can be achieved by design, allowing for multi-modal propagation of topological surface states^28[,29. On the one hand, the system with Chern number N supports N equifrequency loops. These N equifrequency loops are compressed into a folded BZ that is 1/N the size of the original BZ. In this sense, if we are interested in quantities integrated over the BZ, we cannot expect an increase in extensive quantities such us the total field intensity. However, if we are interested in addressing states at a particular wavevector, which is a reasonable constraint in photonic systems, then the modulation has allowed us to address N chiral surface modes with equivalent reciprocal lattice vectors, i.e. achieve unidirectional multiple surface mode operation. The capability of designing photonic systems with large Chern numbers in 3D could find interesting applications in the development of the emergent field of topological lasers^14,49 with a larger number of unidirectional SS.

Second, we showed that any element of the first Chern class vector can be selected by simply changing the magnetization direction, allowing for unique 3D CI/3D CI interfacing combinations as compared to 2D. For example, considering that every planar cut parallel to the magnetization direction is capable of supporting anomalous surface states, it could be worth investigating what occurs at a 3D C_yI/3D C_zI planar interface. Similarly, owing to cubic symmetry of the 3D CI system, it would be possible to tune the relative angular phase in the supercell modulation in order to either break spatial inversion or not, which could lead to interesting surface states at the boundary of an obstructed atomic insulator (OAI) 3D CI and a 3D CI³⁰. This may allow to develop of interesting photonic analogues of axionic responses^26,30 that we are investigating as part of a future work. Another possibility allowed just in 3D, could be to arrange different 3D CIs around an inert core, with the 3D CI composing each panel having Chern vector (C_x, C_y, C_z) oriented to point inwards (e.g. fixing a 3D + C_xI on a left $\hat{{{{{{{{\bf{x}}}}}}}}}$ panel). Such a 3D interfacing arrangement, originally proposed in Ref. ²⁶ as a possible realization of a magneto-electrical (ME) coupler in the field of electronics, has not yet a realization or equivalent in photonics. Analysis of all these challenging designs is left for further investigation.

Third, we showed the TRS breaking parameters required to induce this 3D CI phase can be substantially diminished by the use of larger supercells, which can enable the realization of a 3D CI phase also in photonic systems where the magnetic response is weak or it is not possible to manipulate largely the Weyl points in the BZ. We also intend to emphasize that the strategy we devised in our paper is material agnostic, and can be easily adapted to any to-be-discovered experimental platform. In that sense, our work provides a roadmap to future experimental exploration of topological photonic crystals by showing how to reduce the needed magnetic response.

Finally, we showed that 3D CI photonic phase obtained displays chiral surface states on the planes orthogonal to the magnetization. As a remarkable signature of this, we observed the formation of disjoint equifrequency loops structures associated to the spatial separation of optically-chiral and counter-propagating partners. In conclusion, our system provides a realization of a photonic 3D CI state of matter in a fully cubic platform, with large Chern vectors engineered by design and in a weakly magnetic environment.

Methods

EM section Chern number from Wilson loop approach in 3D

In order to provide a topological characterization of the bulk of the 3D CI photonic phase, we employ a mapping between electronic hybrid Wannier charge centers (WCC)^38,39 and photonic hybrid Wannier energy centers (WEC) in 3D. Hybrid WEC of the type θ_n,y, which are localized in the y-direction and flowing in the k_x transverse direction³⁸, are computed for each fixed k_z plane from the subset of n = 1, . . , N_b bands below the local gap at k_z, as follows³⁹:

$${\theta }_{n,y}({k}_{x})=-{{{{{{{\rm{Im}}}}}}}}{{{{{{\mathrm{log}}}}}}}\,({w}_{n}({k}_{x})),$$

(2)

where w_n are the eigenvalues of the Wilson loop (WL) operator. The real space correspondence in terms of the i-th lattice parameter ∣a_i∣ is given by 2πWEC_n,i ≡ ∣a_i∣θ_n,i. As shown in^50,51, the WL operator can be numerically implemented as a path-ordered product of overlap matrices

$$\hat{W}({k}_{x})= \bar{\mathop{\prod}\limits_{{k}_{{y}_{i}}\in l}} {\hat{M}}^{{k}_{{y}_{i}},{k}_{{y}_{i+1}}},$$

(3)

evaluated on a k_y discretized closed loop l crossing the BZ torus with

$${\hat{M}}_{m,n}^{{k}_{{y}_{i}},{k}_{{y}_{i+1}}}=\langle {u}_{m,{k}_{{y}_{i}}}| {u}_{n,{k}_{{y}_{i+1}}}\rangle .$$

(4)

Here u_m,k(r) represents the periodic part of the electromagnetic Bloch wavefunctions ${{{\Psi }}}_{m,{{{{{{{\bf{k}}}}}}}}}({{{{{{{\bf{r}}}}}}}})=\left(\begin{array}{c}{{{{{{{{\bf{E}}}}}}}}}_{m,{{{{{{{\bf{k}}}}}}}}}({{{{{{{\bf{r}}}}}}}})\\ {{{{{{{{\bf{H}}}}}}}}}_{m,{{{{{{{\bf{k}}}}}}}}}({{{{{{{\bf{r}}}}}}}})\end{array}\right)$, for each mode of momentum k and band m, according to the relation Ψ_m,k(r) = e^{−ik⋅r}u_m,k(r). In order to perform the correct equivalence between WCC and WEC, the scalar product needs to be weighted over the medium permittivity ε, permeability μ and bianisotropy χ tensors^52,53:

$$\langle u| v\rangle =\int {d}^{3}{{{{{{{\bf{r}}}}}}}}{u}^{{{{\dagger}}} }\left(\begin{array}{cc}\varepsilon &\chi \\ {\chi }^{{{{\dagger}}} }&\mu \\ \end{array}\right)v$$

(5)

As detailed in Ref. ⁴⁰, in order to cure the effect of a possible arbitrary phase gained by the fields from numerical evaluation at different points on the l path, we enforce periodic boundary conditions at the endpoints ${k}_{{y}_{i}}$ and ${k}_{{y}_{i}+{b}_{y}}$, where b_y is the reciprocal y lattice vector. Differently from 2D, we can avoid the singularity at ω = 0, k = 0^54,55 by looking only at planar cuts which do not include the Γ point; this is sufficient for determining the Chern vector of our model. From the winding of the hybrid WEC in the BZ, the electromagnetic Chern number can be directly calculated as follows:

$$2\pi {C}_{z}^{EM}=\int\nolimits_{-\pi }^{\pi }\mathop{\sum }\limits_{n=1}^{{\nu }_{occ}}d{\theta }_{n,y}({k}_{x})$$

(6)

i.e. integrating along the k_x closed path across the BZ and summing up the contribution of transverse y WEC for the entire set of ν_occ bands lying below the local gap. Notice that all the definitions remain valid under cyclic permutations (k_i, θ_j, C_k) of the Cartesian indexes (ijk) = (xyz). Fig. 5 displays hybrid WEC for the 3D CI phase of Fig. 1f, comparing their individual to their total net contribution. We note that, for our particular non-bianisotropic (χ = 0) system, it is possible to perform a simplification in the computation, decoupling electric and magnetic fields: under these conditions⁵⁶, the Chern number calculated in terms of just the electric (C^E) or just the magnetic field (C^M) is related to the total electromagnetic Chern number (C^EM) as:

$${C}^{EM}=\frac{{C}^{E}+{C}^{M}}{2}={C}^{E}={C}^{M}.$$

(7)

Analytical models

In this section, we set up two symmetry-adapted models that describe: first, the threefold degeneracy at R and its splitting in k space when an external magnetic field B is applied, to give rise to a pair of Weyl points; second, the merging of the Weyl points by a lattice-commensurate modulation into a gapped topological phase.

Both analytical models are based on the standard group theoretical method of invariants, which can be found in the literature⁵⁷. This allows us to find an expansion in powers of the wave vector k of the photonic energy bands ω, able to replicate the photonic modes dispersion. In order to do so, we construct an effective energy dispersion operator H(k), expressed in terms of the space group irreducible representations bases, able to capture the photonic modes symmetry properties. In the photonic context, H(k) can be viewed as a perturbative expansion of the Maxwell-Bloch operator acting on the electromagnetic fields in the first order formulation of Maxwell’s equations⁵⁸. For other applications of this approach to photonic systems, see e.g. Refs. ^3,59.

In what follows, we adopt the notation convention taken in the Bilbao Crystallographic Server (BCS)⁶⁰, unless otherwise stated, and express reciprocal lattice vectors in reduced units 2π/∣a∣ = 1.

Threefold degenerate model at R

The following model describes the local behavior of the modes which are threefold degenerate at point R = (1/2, 1/2, 1/2) (see Fig. 6(a)). From the numerical computations, we know that this degeneracy is related to the three dimensional small representation ${R}_{4}^{-}$ of the little group of R, given the transformation properties of these modes under the elements of the space group $Pn\bar{3}m$ (No. 224).

Fig. 6: Comparison of the numerical band-dispersion with the analytical model.

Extracting the coefficients of the analytical k ⋅ p model in the vicinity of R. Here, we are considering an extrapolation of the model at a distance of δk = 0.008 from R. The empty circles are for numerical computed bands v_1,2,3 while the lines for the analytical dispersion w_1,2,3. a In presence of TRS, we conclude that b₀ = −2.9a₀ > 0. b In a weak field η = 0.5, we obtain that α₀ ~ − β₀ since the third band does not move in energy and α₀ ~ 0 since the vertical displacement of the two lowest degenerate modes is equal and opposite. Weyl points appear at positions: ${k}_{z}^{\pm }=\pm \sqrt{\frac{| {\delta }_{0}{B}_{z}| }{{b}_{0}}}$, which here is at ${k}_{z}^{\pm }=\pm \! 0.004$. The δ₀B_z > 0 and δ₀B_z < 0 cases are equivalent upon reversing the two lowest photonic modes. Source data are provided as a Source Data file.

Following the method of invariants, we first note that the product ${R}_{4}^{-* }\times {R}_{4}^{-}$ is decomposed into small irreducible representations (irreps) at Γ as:

$${R}_{4}^{-* }\times {R}_{4}^{-}={{{\Gamma }}}_{1}^{+}+{{{\Gamma }}}_{3}^{+}+{{{\Gamma }}}_{4}^{+}+{{{\Gamma }}}_{5}^{+}$$

(8)

using the character orthogonality relations and where '*' denotes complex conjugation.

A general state in this three-band space can be expanded in the basis $\{\left|{\phi }_{i}\right\rangle \}$ of the ${{{{{{{{\rm{R}}}}}}}}}_{{{{{{{{\rm{{4}}}}}}}^{-}}}}$ representation as:

$$\left|{\psi }_{{{{{{{{\bf{k}}}}}}}}}\right\rangle ={c}_{i}({{{{{{{\bf{k}}}}}}}})\left|{\phi }_{{{{{{{{\bf{k}}}}}}}}}^{i}\right\rangle ,$$

(9)

adopting the Einstein summation convention. The energy expectation value is a scalar invariant, which is computed as:

$$\langle H\rangle =\left\langle {\psi }_{{{{{{{{\bf{k}}}}}}}}}\right|H\left|{\psi }_{{{{{{{{\bf{k}}}}}}}}}\right\rangle ={c}_{i}^{* }({{{{{{{\bf{k}}}}}}}}){c}_{j}({{{{{{{\bf{k}}}}}}}})H{({{{{{{{\bf{k}}}}}}}})}_{ij}.$$

(10)

We seek combination of bilinears ${c}_{i}^{* }{c}_{j}$ transforming as the irreps above and take the Hermitian scalar product with functions of k and B with the same symmetry properties. From each term with ${c}_{i}^{* }{c}_{j}$ in this scalar product, it is easy to obtain the matrix elements of H(k, B). The energy scalar is written in this scalar product form:

$$\langle H\rangle =\mathop{\sum}\limits_{\alpha ,i}{C}_{i}^{\alpha }{({q}_{i}^{\alpha })}^{* }\cdot {p}_{i}^{\alpha }$$

(11)

where the sum runs over the irreps in the decomposition and {q_i} and {p_i} are the symmetry-adapted bases of the state coefficients and k and B, respectively. The coupling constants ${C}_{i}^{\alpha }$ are parameters of the model.

To find the bases of bilinears in the wave coefficients transforming as the irreps above, we use the representation ${\rho }_{{R}_{4}^{-}}$ of the generators (omitting inversion for brevity, as it is represented by the negative identity matrix):

$$\begin{array}{ll}{\rho }_{{R}_{4}^{-}}({2}_{001})&=\left(\begin{array}{rcl}-1&0&0\\ 0&-1&0\\ 0&0&1\end{array}\right)\,{{\mbox{,}}}\,\ {\rho }_{{R}_{4}^{-}}({3}_{111}^{+})=\left(\begin{array}{rcl}0&1&0\\ 0&0&1\\ 1&0&0\end{array}\right)\,{{\mbox{,}}}\,\\ &{\rho }_{{R}_{4}^{-}}({2}_{110})=\left(\begin{array}{rcl}0&1&0\\ 1&0&0\\ 0&0&-1\end{array}\right),\end{array}$$

(12)

where, for convenience, we have labeled the matrices by the rotation part of the symmetry element only. Note that these differ from the BCS data by a permutation of the basis, chosen to rearrange H in a more convenient form. Since it is a non-symmorphic space group, some operations have fractional translations. In Seitz notation, these are:

$$\left \{{2}_{001}| \frac{1}{2},\frac{1}{2},0\right\},\left \{{2}_{110}| \frac{1}{2},\frac{1}{2},0 \right \}.$$

(13)

We then find functions of k and B with the same transformation properties, up to second order in the wave vector. In principle, the magnetic field could be strong. Therefore, the criterium to choose the maximum power of B that is included for each order in k is to exhaust all the possibilities in the irrep decomposition. This way, we ensure that all the couplings allowed by symmetry are included for a given order in the wave vector. Finally, we require H(k, B) to be Hermitian.

Following this procedure, we find that the most general expression for the energy operator is:

$$H=\, ({a}_{0}{k}^{2}+{\alpha }_{0}{B}^{2}){{\mathbb{1}}}_{3}+i{\delta }_{0}{\varepsilon }_{jkl}{B}_{l}+{b}_{0}\left(\begin{array}{rcl}{k}_{x}^{2}&0&0\\ 0&{k}_{y}^{2}&0\\ 0&0&{k}_{z}^{2}\end{array}\right)\\ + \,{\beta }_{0}\left(\begin{array}{rcl}{B}_{x}^{2}&0&0\\ 0&{B}_{y}^{2}&0\\ 0&0&{B}_{z}^{2}\end{array}\right)+{c}_{0}\left(\begin{array}{rcl}0&{k}_{x}{k}_{y}&{k}_{x}{k}_{z}\\ {k}_{x}{k}_{y}&0&{k}_{y}{k}_{z}\\ {k}_{x}{k}_{z}&{k}_{y}{k}_{z}&0\end{array}\right)\\ +\, {\gamma }_{0}\left(\begin{array}{rcl}0&{B}_{x}{B}_{y}&{B}_{x}{B}_{z}\\ {B}_{x}{B}_{y}&0&{B}_{y}{B}_{z}\\ {B}_{x}{B}_{z}&{B}_{y}{B}_{z}&0\end{array}\right),$$

(14)

where k = (k_x, k_y, k_z) is the wave vector measured from the point R and we employ real coefficients (Latin when referring to k and Greek to B).

One can check that the energy operator is invariant under the little-group symmetries as it verifies, for every operation g = {R∣t}:

$${\rho }_{{R}_{4}^{-}}(g)H{\rho }_{{R}_{4}^{-}}{(g)}^{-1}=H(R{{{{{{{\bf{k}}}}}}}},R{{{{{{{\bf{B}}}}}}}}).$$

(15)

We also have imposed that the model be invariant when both the system and the external magnetic field B are transformed by time reversal Θ. We can express the TR operation as Θ = Uκ, where U is a unitary matrix and κ is the complex conjugation operator. Then, the TRS condition reads:

$${{\Theta }}H({{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{B}}}}}}}}){{{\Theta }}}^{-1}=U{H}^{* }({{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{B}}}}}}}}){U}^{-1}=H(-{{{{{{{\bf{k}}}}}}}},-{{{{{{{\bf{B}}}}}}}})$$

(16)

where the unitary 3 × 3 matrix part has the simple form $U={{\mathbb{1}}}_{3}$. Evaluating the model at the point R, in the presence of TRS, allows us to fix some of the coefficients by comparing with the numerical simulations, as described in Fig. 6. This yields b₀ ~ 2.9a₀.

When a magnetic field is applied along one of the coordinate axes ${{{{{{{\bf{B}}}}}}}}={B}_{i}{\hat{{{{{{{{\bf{x}}}}}}}}}}_{i}$, the energy dispersion of the three photonic modes along the line parallel to the field is:

$$\left\{\begin{array}{ll}{\omega }_{1}=({a}_{0}+{b}_{0}){k}_{i}^{2}+({\alpha }_{0}+{\beta }_{0}){B}_{i}^{2}\quad &\\ \kern-2.1pc {\omega }_{2}={a}_{0}{k}_{i}^{2}+{\alpha }_{0}{B}_{i}^{2}-{\delta }_{0}{B}_{i}\quad &\\ \kern-2.1pc {\omega }_{3}={a}_{0}{k}_{i}^{2}+{\alpha }_{0}{B}_{i}^{2}+{\delta }_{0}{B}_{i}.\quad &\end{array}\right.$$

(17)

This further fixes α₀ ~ −β₀ ~ 0 and shows that the magnetic field fully lifts the threefold degeneracy. We see in Fig. 6(b) that the band curved upwards in energy will cross with one of the remaining two, giving rise to a Weyl point. The strength of the magnetic field tunes where this crossing happens along this line, according to the expression:

$${k}_{i}^{\pm }=\pm \sqrt{\frac{| {\delta }_{0}\ {B}_{i}| }{{b}_{0}}}.$$

(18)

The same happens in the opposite direction along the same line, hence the ± sign. This shows that a Weyl dipole appears along the line parametrized by k_i and that the position of the nodes can be tuned by the magnetic field strength B_i.

Coupling of Weyl points by supercell modulation

Because the validity of the previous analysis is limited to the neighborhood of the point R, we construct another model that expands directly around the Weyl points. In particular, we can fix the magnetic field to ${{{{{{{\bf{B}}}}}}}}={B}_{z}{{{\hat{{{{{\boldsymbol{z}}}}}}}}}$ and tune the strength B_z to create a pair of Weyl points at K_1,2 = (1/2, 1/2, ±1/4). The Weyl nodes can then be coupled by a supercell modulation that doubles the real-space unit cell in the ${{{\hat{{{{{\boldsymbol{z}}}}}}}}}$ direction. This corresponds to the uniaxial N_W = N = 2 case in the main text, i.e. a (1, 1, 2) supercell. A generalization to cubic (N, N, N) supercells is straightforward, since each Cartesian component of the modulation can be turned on independently.

The compatibility relations from R into the T = (1/2, 1/2, u) line yield:

$${R}_{4}^{-}(3)\to {T}_{1}(1)+{T}_{5}(2),$$

(19)

where the dimensions of the small irreps are in parentheses. The magnetic field splits the states of T₅ and one of them is degenerate with the T₁ state at K_1,2. We set up a model that describes these six photonic modes and its coupling by a cell modulation commensurate with the lattice.

Both K₁ and K₂ belong to the same star and are related by inversion. The representation that acts on the six photonic states is obtained from the space group representation induced from the direct sum T₁ + T₅. We then restrict this representation only to the K_1,2 arms and consider the elements which either leave K_i invariant or relate one to the other. The subspace of these two arms is invariant under all these elements, which form a group that we denote by G_W.

Let us call D the representation of G_W so obtained. D is divided into blocks arising from the T₁ and T₅ irreps, hence we write D = D₁ + D₅. The direct product of the full space group irrep can be used to find the reduction of the product D^* × D into small irreps of the little groups at Γ and ${{{{{{{{\bf{X}}}}}}}}}_{3}=\frac{{{{{{{{{\bf{b}}}}}}}}}_{z}}{2}=(0,0,1/2)$. The result is shown in Table 2. Note that the label for the X point differs from that in the main text (where it is called X₂) and was chosen for consistency with the BCS notation.

Table 2 Decomposition of the product of representations D^* × D. Each row labels the decomposition of the term-wise product, considering that D = D₁ + D₅ is derived from a space group irrep induced from the sum of small irreps T₁ + T₅, which is then restricted to the arms K_1,2. The terms in the reduction are small irreps of the little groups at Γ and X₃ = (0, 0, 1/2). The labels for the X irreps are those from the point ${{{{{{{\bf{X}}}}}}}}\equiv {{{{{{{{\bf{X}}}}}}}}}_{2}=\frac{{{{{{{{{\bf{b}}}}}}}}}_{y}}{2}=(0,1/2,0)$, since X₃ is in the star of X.

We may divide the 6 × 6 energy operator matrix into 3 × 3 blocks. Then, the diagonal blocks are identified with the Weyl points and the off-diagonal ones with the modulation that couples the states at both nodes. In view of Table 2, we can also anticipate that the Γ irreps will yield the elements of the diagonal blocks, and the X_m (m = 1, 2, 3, 4) irreps the off-diagonal ones. This is consistent with lattice translations having non-trivial representation for the X_m irreps, which means that they couple non-equivalent points in the BZ.

The matrices of D that are needed to find the symmetry-adapted bases are the following (again, we omit the representation matrix of inversion for brevity):

$$\rho _{D_1}(2_{001}) \, = \, \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right),\rho _{D_1}(m_{010}) = \left( \begin{array}{cc} {e^{ - i\pi /4}} & 0 \\ 0 & {e^{i\pi /4}} \end{array} \right),\\ \rho _{D_1}\left( {4_{001}^ + } \right) \, = \,\left( \begin{array}{cc} {e^{i\frac{\pi }{4}}} & 0 \\ 0 & {e^{ - i\pi /4}} \end{array} \right)$$

(20)

and

$$\begin{array}{ll}& \kern-7pc {\rho }_{{D}_{5}}({2}_{001})=\left(\begin{array}{cc}-{1}_{2}&0\\ 0&-{1}_{2}\\ \end{array}\right)\,{{\mbox{,}}}\,\\ & \kern-1pc {\rho }_{{D}_{5}}({m}_{010})=\left(\begin{array}{cccc}0&{e}^{-i3\pi /4}&0&0\\ {e}^{i\pi /4}&0&0&0\\ 0&0&0&{e}^{-i\pi /4}\\ 0&0&{e}^{i3\pi /4}&0\\ \end{array}\right)\,{{\mbox{,}}}\,\\ &{\rho }_{{D}_{5}}({4}_{001}^{+})=\left(\begin{array}{cccc}0&{e}^{-i3\pi /4}&0&0\\ {e}^{-i3\pi /4}&0&0&0\\ 0&0&0&{e}^{-i\pi /4}\\ 0&0&{e}^{-i\pi /4}&0\\ \end{array}\right).\end{array}$$

(21)

Following the same procedure as in the Methods section, we obtain the matrix representation of the energy operator. As a last step, we impose TRS. The unitary part of the TR operation Θ = Uκ is in this case

$$U=\left(\begin{array}{cc}0&1\\ 1&0\end{array}\right)\otimes \left(\begin{array}{ccc}1&0&0\\ 0&-1&0\\ 0&0&1\end{array}\right).$$

(22)

We are interested in modulations implemented by physically altering the dielectric structure of the crystal. Therefore, the modulation itself is considered to transform trivially under Θ.

TRS forbids one of the couplings for each one of the X₂, X₃ and X₄ modulations. Furthermore, analyzing the effect of these three by numerically diagonalizing the H matrix shows that only modulations transforming in the X₁ representation can gap out the Weyl points. For examples of modulations that do not open a gap, see SN 12.

The representation matrices used to obtain the X₁ modulation terms are the following:

$$\begin{array}{ll}{\rho }_{{X}_{1}}({2}_{001})&=\left(\begin{array}{cc}1&0\\ 0&1\\ \end{array}\right)\,{{\mbox{,}}}\,\ {\rho }_{{X}_{1}}({m}_{010})=\left(\begin{array}{cc}0&-1\\ 1&0\\ \end{array}\right)\,{{\mbox{,}}}\,\\ &{\rho }_{{X}_{1}}({4}_{001}^{+})=\left(\begin{array}{cc}0&-1\\ 1&0\\ \end{array}\right).\end{array}$$

(23)

We present the expression along the T line of H(k_z, B_z) with only the X₁ couplings included:

$$H({k}_{z},{B}_{z})=\left(\begin{array}{cc}{W}_{+}&{X}_{1}\\ {X}_{1}^{{{{\dagger}}} }&{W}_{-}\end{array}\right),$$

(24)

$${W}_{\pm }=\left(\begin{array}{rcl}a{k}_{z}^{2}\pm b{k}_{z}+\alpha {B}_{z}^{2}&0&0\\ 0&c{k}_{z}^{2}\pm d{k}_{z}+\beta {B}_{z}^{2}&{B}_{z}\delta \pm \gamma {k}_{z}{B}_{z}\\ 0&{B}_{z}\delta \pm \gamma {k}_{z}{B}_{z}&c{k}_{z}^{2}\pm d{k}_{z}+\beta {B}_{z}^{2}\end{array}\right),$$

(25)

$${X}_{1}({p}_{1},{q}_{1})=\left(\begin{array}{rcl}{C}_{1}({p}_{1}-i{q}_{1})&0&0\\ 0&{C}_{2}{p}_{1}-i{C}_{3}{q}_{1}&0\\ 0&0&-i{C}_{2}{q}_{1}+{C}_{3}{p}_{1}\end{array}\right),$$

(26)

where the coordinates (p₁, q₁) transform as X₁ and parametrize the modulation strength, and C_1,2,3 are real coupling constants, while the rest of parameters are also real. The k_z component is taken from the point where the Weyl points merge after the cell folding.

The effect of this modulation is visualized in Fig. 7. The band-gap opened via supercell modulation in the TRS broken system is shown to be a Chern gap in the main text, by numerical means.

Fig. 7: Effect of the X₁ modulation and Weyl points coupling.

From the analytical model, a four-fold degeneracy point is achieved at the folding condition ${B}_{z}=\frac{\delta }{\alpha -\beta }$, as shown in panel a. When a X₁ modulation is introduced, the Weyl dipole couples and a band-gap is opened, as showed in panel b. Parameters: a = −7, b = −30, α = 1, c = 1, d = 30, γ = 1, β = 0.45, δ = 7.5, C_1,2,3 = 1, p₁ = 7, q₁ = 0.

From the analytical model, we also observe that, to exactly superimpose the Weyl points, we need to tune the magnetic field to the folding condition:

$${B}_{z}=\frac{\delta }{\alpha -\beta },$$

(27)

as can be seen by diagonalizing the matrix H(0, B_z).

As stated before, this model addresses the case where N = N_W = 2. When N_W = 2 is fixed but N = 2n with n integer (see Fig. 3), the modulation belongs to point (0, 0, 1/N) and must enter at order n in H. Therefore, unless N_W = N = 2, the modulation will belong to the high-symmetry line Δ. The expression for the modulation for every N = N_W > 2 is given in the SN 13.

Example of cell modulation

We use the projectors onto the i-th basis element in the space of the irrep X₁:

$${P}_{ii}\propto \mathop{\sum}\limits_{g\in G}{X}_{1}^{* }{(g)}_{ii}g$$

(28)

where g runs over the little co-group at X₃ = (0, 0, 1/2) and we disregard any normalization factors. Applying these to an arbitrary function f(z), we find:

$$\begin{array}{ll}&{X}_{1} :\ \left\{\begin{array}{l}{P}_{11}f\propto f(z)+f(-z)\quad \\ {P}_{22}f\propto f(z)-f(-z)\quad \end{array}\right.\\ {X}_{2},{X}_{3},{X}_{4} \kern-2pc &:\ {P}_{ii}f=0,\quad i=1,2.\end{array}$$

(29)

Therefore, given functions of z that under lattice translations obey ${{{{{{{\bf{T}}}}}}}}f={e}^{i{{{{{{{{\bf{X}}}}}}}}}_{3}\cdot {{{{{{{\bf{T}}}}}}}}}f$, those that provide a basis for this irrep are one even and one odd, respectively. This proves that a modulation of the radius of the rods ${{\Delta }}r=r(z)-{r}_{0}={r}_{m}\cos (2\pi z/N| a| )$ with N = 2 belongs to X₁. In particular, it is parametrized by (p, q) = (p, 0) in the model given in Methods section. We also note that, when the modulation is cubic, i.e. ${{\Delta }}r=r(x,y,z)-{r}_{0}={r}_{m}[\cos (\pi x/| a| )+\cos (\pi y/| a| )+\cos (\pi z/| a| )]$, it is only the z dependent part that is responsible for gapping the Weyl points generated in a B_z field.

We also note that our derivation was performed employing Hermitian perturbations. However, via the introduction of non-Hermitian perturbations in the model, it could be possible to incorporate the effect of losses (or gain) in the system. Non-Hermitian terms usually lead to a spread of the Chern bands along the imaginary axis, which therefore transform into band regions in the complex plane. As long as the effect is limited enough not to lead to merging of the two bands, it is generally possible to separate the two bands by a line-gap. Indeed, as shown in Ref. ⁶¹, in presence of a line-gap, Chern insulators are stable with respect to non-Hermitian lossy effects.