Designing Photonic Topological Insulators with Quantum-Spin-Hall Edge States using Topology Optimization

Designing photonic topological insulators is highly non-trivial because it requires inversion of band symmetries around the band gap, which was so far done using intuition combined with meticulous trial and error. Here we take a completely different approach: we consider the design of photonic topological insulators as an inverse design problem and use topology optimization to maximize the transmission through an edge mode with a sharp bend. Two design domains composed of two different, but initially identical, C$_\text{6v}$-symmetric unit cells define the geometrical design problem. Remarkably, the optimization results in a photonic topological insulator reminiscent of the shrink-and-grow approach to quantum-spin-Hall photonic topological insulators but with notable differences in the topology of the crystal as well as qualitatively different band structures and with significantly improved performance as gauged by the band-gap sizes, which are at least 50 \% larger than previous designs. Furthermore, we find a directional beta factor exceeding 99 \%, and very low losses for sharp bends. Our approach allows for the introduction of fabrication limitations by design and opens an avenue towards designing PTIs with hitherto unexplored symmetry constraints.


I. 1. INTRODUCTION
Defects in photonic crystals (PCs) such as cavities or waveguides allow for confining light to subwavelength dimensions [1] and has enabled a wealth of studies of light-matter interaction at the nanoscale [2]. However, photonic crystals are also very sensitive to structural disorder, which results in substantial variations in the resonance frequencies and quality of PC cavities [3] and backscattering-induced Anderson localization of PC waveguides [4], which are serious impediments for applications of PCs. Photonic topological insulators (PTIs) provide a fundamentally different approach to confinement where light is confined as symmetry-protected edge states between PCs with different band topologies and such edge modes may be robust against backscattering [5,6].
While PTIs have been realized in several photonic platforms featuring broken time-reversal symmetry [7][8][9], which allows unidirectional modes similar to the quantum Hall effect [10], these approaches require materials and geometries incompatible with inorganic semiconductors and planar technology, e.g., III-V materials for active communication devices or solid-state quantum optics and silicon photonics for chip-scale routing and optical interconnects. More recently, time-reversal-symmetric all-dielectric PC-based PTI designs for the photonic quantum-spin-Hall effect have been proposed [11], which was later extended to planar PTIs [12] in which a number of chiral quantum optical effects [13] have been demonstrated experimentally [14,15]. As an alternative approach, breaking parity symmetry leads to valley-Hall PTIs [16,17].
A key ingredient in PTI design is band inversion, i.e., the symmetry of the photonic bands must be swapped around the photonic band gap in order to induce the topological transition forming the edge state. This is in general a non-trivial task and thus remarkably few designs are available in the literature. The starting point for all planar time-reversal-symmetry-protected quantum-spin-Hall PTI designs so far is photonic graphene, which features a doubly degenerate Dirac point that upon suitable lattice perturbations may be gapped in topologically non-trivial ways. Notably, all dielectric photonic-crystal-based quantumspin-Hall PTI designs to date employ the shrink-and-grow strategy first used by Wu and Xu [11]. Here we take an entirely different approach: we consider the design of a PTI as an inverse numerical optimization problem seeking the maximum transmission of edge states between two crystals with C 6v symmetry by manipulating the crystal geometry. We solve this problem using density-based topology optimization (TO) [18], an optimization-based design method utilizing adjoint sensitivity analysis [19] and gradient-based optimization algorithms [20] to efficiently solve design problems with potentially billions of degrees of freedom. Density-based topology optimization has previously been successfully applied across numerous areas, including recent optimization and design of complex structures and materials in mechanics [21][22][23][24], optics and photonics [25][26][27][28], thermo-fluidics [29], and acoustics [30,31]. Very recently, TO was used to devise new designs of acoustic topological insulators [32] and here we apply this approach for the first time to the design of PTIs. Remarkably, our approach generates topologically non-trivial bands that resemble those obtained from the use of the shrink-and-grow strategy although these appears spontaneously through the optimization process and not as a imposed design constraint. While our approach thus qualitatively reproduces concepts and ideas already available in the literature, the resulting unit cells are entirely different with non-trivial geometric features, which turn out to significantly improve figures of merit: The relative band gaps of the two crystal phases are 19% and 6%, we find a directional beta-factor > 99%, and the bending losses are extremely small.

II. 2. MODEL
Designing a PTI can be viewed as a problem of determining the material configuration in two PC phases, which, due to their different Bloch-function symmetries, form a topological edge state at their interface. The starting point of our model is a geometric domain split between two phases in a checkerboard topology as shown in Figure 1A. In a PTI, light impinging from P1 can only scatter to P2 and P4 but not to P3 and thus our hypothesis is that by maximizing the transmission to P2 and P4 while minimizing the transmission to P3, the TO will be forced to generate a PTI. As we will show, this indeed turned out to be the case. We would like to note that we have also investigated other design domain topologies and find that the inclusion of the sharp bend between the two phases is essential to generate a PTI. Here we consider two PC phases constructed from hexagonal unit cells exhibiting C 6v -symmetry [33], leaving a twelfth of each hexagonal cell designable for each phase as shown in Figure 1B. The PC phases are constrained to consist of two materials, silicon and air. The numerical model used in the design procedure detailed in the following considers a spatial domain, Ω ∈ R 2 , cf. Figure 1C. For the design procedure, the two phases are distributed in Ω as shown using black and orange outlines. In this configuration, the PC phases exhibit 180 o rotational symmetry around the centre of the hexagonal material slab. Light impinging from P1 reaches a crossing between the different crystal phases, and at the crossing, light cannot propagate forwards, it can only scatter to the left or right. (B) Symmetry constraints for a single unit cell, which in both phases is restricted to C6v symmetry. The designable part of the unit cell is denoted ΩD,i (dark gray) and mirrored across the bold black line. The triangular region (light gray + dark gray) is replicated by 60 o · n, n ∈ {1, 2, ..., 5} rotations to fill the hexagonal unit cell. (C) The full optimization domain showing the PC phase distribution. Light is impinging on port P1 from an excitation source whose position is indicated with a black dot highlighted by a red circle at port P1 and the intensity at ports P2-P4 is recorded. The target function maximizes the intensity at P2 and P4 while minizing the intensity at P3. The regions containing the first and second crystal phases are outlined using orange and black denoted Phase 1 and Phase 2, respectively. A number of subdomains, highlighted using various colors, We model the light propagation in the frequency domain, where E denotes the electric field, k 0 = 2πν c the free-space wave number with ν being the frequency and c the speed of light in vacuum, r is the material-dependent relative permittivity and S denotes the excitation source. The model domain Ω is truncated using a perfectly matched layer (PML) [34]. In the design process, artificial attenuation is introduced in the subdomain surrounding the hexagonal slab to prohibit the E-field from propagating around the outside of the slab from P1 to P2-P4. The attenuation is introduced in Eq. (1) through a complex value of the relative permittivity, r = 1 + 10i, in the region surrounding the hexagonal slab.
The problem of designing the geometries for the two PC phases is recast as a continuous optimization problem and solved using topology optimization as detailed in [32] and references therein. In short, the transmission from P1 to P2 and P4 is sought maximized for a number of target frequencies while simultaneously minimizing the transmission from P1 to P3 and ensuring The final PTI material configuration. Both the shape and topology of the second phase has changed from the initial guess while only the shape has changed for the first phase. The interface between the two phases is highlighted in red.
that both PC phases have bulk bandgaps at the targeted frequencies. The power flow to P2-P4 as well as inside the bulk of the PC phases for a given excitation S(r) and material configuration, r (r), is estimated by evaluating the spatially normalized integral of the Poynting vector over Ω defined in Eq. (4). The target function, Φ Total , as defined in Eq. (2) is maximized by (re)distributing silicon and air in the periodic unit cells of each crystal phase, Ω D,i , i ∈ {1, 2} (see Figure 1B), while ensuring that constraints enforcing bandgaps in the bulk of the PC phases, given in Eq. (3), are satisfied, with, Here P denotes the Poynting vector and τ a set of scaling constants. The subdomains Ω Φ are all shown in Figure 1C.
The target function consists of three contributions: Φ Max,1 measuring the power flow to P2 and P4 (the region Ω Max,1 ), Φ Max,2 measuring the power flow at the centre of the hexagonal slab (the region Ω Max,2 ), and Φ Min , which measures the power flow to P3 (the region Ω Min ). The four constraint functions: Φ BG,j (ξ), j ∈ {1, 2, 3, 4} are computed by evaluating the power flow in the regions Ω BG,j . The choice of the target function in Eq. (2) and the configuration of Ω Max,i and Ω Min assures that by solving the optimization problem the vast majority of the power coupled into the hexagonal slab at P1 will flow to P2 and/or P4 with only minimal power flowing to P3. The inclusion of Φ Max,2 in the target function ensures that light propagates to the center of the domain in the initial steps of the optimization but has negligible influence on the final result. The choice of Eq. (3) and Ω BG,j assures that only a negligible amount of power flows into the bulk of the PC phase, effectively ensuring bulk bandgaps in both PC phases at the targeted frequencies.

III. 3. RESULTS
We model a silicon slab perforated with air holes, which is implemented in a 2D using a reduced effective permittivity of rSi = 9.61 for silicon while rair = 1 is assumed for the air. A lattice constant of a = 735 nm is chosen for the hexagonal unit cells, cf. Figure 1A. The excitation, S(x), is chosen to be a transverse electric (TE) polarized point source. The design problem defined by Eqs. (1)-(3) is solved for the two frequencies ν 1 = 184 THz and ν 2 = 187 THz (approximately corresponding to the wavelengths λ 1 = 1630 nm, λ 2 = 1604 nm) using a min/max approach [18]. The initial material distribution is chosen to be a photonic crystal structure with equilateral triangular holes and a bandgap in the considered frequency range, inspired by the recent work of Barik et al. [12], see Figure 2A. Note that the starting material configuration includes a smooth transition in the dielectric function around the edges, which is a result of the filtering process used in the TO algorithm, see e.g. [30].
The material distributions for both PC phases resulting from solving the design problem are shown in Figure 2B. The interface between the phases is highlighted in red, and black/white corresponds to silicon/air, respectively. From the figure it is clearly seen that the shape of the air holes in the left PC phase has changed as a result of the design process, while for the right PC phase both a shape change and a change of the topology of the air holes has occurred.
To demonstrate the topological behaviour of the designed PTI, Figure 3 shows the max-normalized electric field intensity on a logarithmic scale, for the chosen excitation at ν = 184 THz (λ = 1630 nm), for the initial and final material distributions. From Figure 3A it is clearly observed that the initial PC configuration exhibits a bandgap at 184 THz with the energy dropping several orders of magnitude as the field extends into the crystal. In contrast, for the topological edge-state generated by our TO design shown in Figure 3B, the vast majority of the electric field entering at P1 propagates to either P2 or P4 while very little energy propagates to P3 (≈ 1%), indicating the successful design of a back-scattering-suppressing structure. Chiral quantum optical effects are pronounced in PTI edge states due to the spin-momentum locking of the quantum-spin-Hall states. This effect is very prominent in our geometry, as shown in Figure 4A, which shows the emission intensity from emitters  with opposite chirality (spin), placed at (x s , y s ) = (0nm, 70nm) relative to the centre of the unit cell of Phase 1 closest to the interface. By extracting the emitted intensity going left or right for a given emitter chirality, we extract directional β-factors for left-and right-hand circularly polarized emitters as with Ω L and Ω R highlighted in the bottom panel of Figure 4A using green and blue respectively. We exctract β L = β R > 0.995, corresponding to a perfectly deterministic directional light-matter interface [13] within the numerical error of our calculations. As a demonstration of the back-scattering-suppressed field propagation past sharp bends offered by the PTI, two additional slab models are constructed. The first model consists of a slab with four 120 o bends along the interface between the two PC phases, the geometry is shown in the bottom panel of Figure 4B. The second model consists of a straight line interface between the two PC phases identical to the one shown in the bottom panel of Figure 4A. For both models, a point source emitting a TE polarized field at 184 THz is placed outside the slab on the left hand side and aligned with the PC phase interface position. The resulting electric-field intensities in the slabs are shown in the top and middle panels of Figure 4B. The top panel shows the slab with bends while the middle shows the slab with the straight interface. An inset showing a magnified image of the output is included at the top right of both panels, showing qualitatively identical intensity. To obtain a quantitative measure of the loss from traversing the four bends, the integral of |E| 2 in the regions highlighted in blue in the bottom panel of Figure 4A and purple in the bottom panel of Figure 4B, is computed. Hereby, we find that the loss is less than 0.04 dB per bend.
Finally, we have calculated the band structures of both PC phases as well as the dispersion relation of the edge state. The result is presented in Figure 5 where Figure 5A shows a unit cell from each of the two PC phases while Figure 5B shows the associated band structures for the lowest six bands. The left column shows the full band structures and the right column shows a zoom around the target frequencies with an overlay highlighting the bulk bandgaps. From the band structure for phase 1 it is seen that bands 2 and 3 as well as bands 4 and 5 experience a degeneracy at the Γ-point. For phase 2, a degeneracy is only observed at the Γ-point for bands 2 and 3 whereas band 4 has a bandgap on both sides. The degeneracies at the Γ-point form the foundation for the formation of the PTI edge state as the two PC phases are brought into contact. Figure 5C shows the band diagram for the supercell, obtained by imposing periodic boundary conditions on the left/right of the supercell and Neumann boundary conditions on the top/bottom. The frequencies targeted in the topology optimization process are included as horizontal magenta lines. The introduction of the Neumann boundary conditions leads to a number of artificial standing wave Bloch-modes becoming solutions to the eigenvalue problem solved to obtain the band structure. To verify that these modes are indeed artificial Bloch-modes, we have calculated band diagrams for each of the two PTI phases as shown in Figure 5D and observe that these modes clearly are numerical artefacts stemming from the resonances in each of the two phases. By filtering these artificial modes away, one arrives at a pair of counter-propagating edge modes inside the bulk-bandgap region highlighted in Figure 5C using blue and red according to the pseudo-spin associated with each mode. The horizontal purple lines have been added in the diagram to denote the frequencies at which the PTI was designed. The band structure in Figure 5C clearly shows that interface edge-states exist in the bulk bandgap for the two PC phases. A small bandgap is found in the edge-states, which may be attributed to the interface itself breaking the C 6v symmetry [15]. The PTI phases have a central frequency of 188.6 THz (203.3 THz) and bandgap of 11.3 THz (39.5 THz) for Phase 1 (Phase 2), which corresponds to relative bandgap widths of 0.06 and 0.19, respectively. This is 50 % or more larger than the relative bandgaps found by Barik et al. [12], which directly shows the power of employing TO to design PTIs.

IV. 4. CONCLUSION
We have used inverse design by TO to generate a PTI structure without imposing any geometrical constraints beyond the C 6v symmetry. It is interesting to note that although our algorithm is set up to maximize the transmission of light through a crystal with two different phases including a sharp bend, a number of interesting features and phenomena emerge as indirect benefits, such as larger band gaps than previous designs, extremely low bending losses, as well as (within the numerical error) unity directionality of embedded emitters. Neither of these were design constraints but simply emerge as a result of the TO, which is a direct testimony to the power of TO when applied to complex inverse problems. Notably, it would be rather trivial to extend our formulation to, e.g., also target a widening of the bandwidth by introducing more target frequencies which are spread further apart but perhaps even more importantly, our method may be employed to generate entirely new PTI designs based on different crystal symmetries and constraints, which is an exciting prospect for further work.