Novel topological beam-splitting in photonic crystals

We create a passive wave splitter, created purely by geometry, to engineer three-way beam splitting in electromagnetism in transverse electric polarisation. We do so by considering arrangements of Indium Phosphide dielectric pillars in air, in particular we place several inclusions within a cell that is then extended periodically upon a square lattice. Hexagonal lattice structures more commonly used in topological valleytronics but, as we discuss, three-way splitting is only possible using a square, or rectangular, lattice. To achieve splitting and transport around a sharp bend we use accidental, and not symmetry-induced, Dirac cones. Within each cell pillars are either arranged around a triangle or square; we demonstrate the mechanism of splitting and why it does not occur for one of the cases. The theory is developed and full scattering simulations demonstrate the effectiveness of the proposed designs.


Introduction
Beam-splitters play a key role in many optical devices, in particular interferometers, with a consequently broad variety of applications ranging from astrophysics [1] to quantum computing [2] as well as to many areas of optical electronics such as optical modulators [3] for fibre optic telecommunications amongst much more. Indeed to achieve complete control over the flow of light, power division and redirection devices are required, of which beam-splitters are those most commonly utilised; a recurrent theme is the desire to have broadband lossless splitters, capable of multiple and tunable re-direction, and with dimensions that are sub-micron, or at most microns, in size.
Several different beam-splitting approaches have been successfully implemented ranging from coupled bent dielectric slab waveguides or ridge waveguides atop substrates, to photonic crystal/ grating devices using waveguide or self-collimation or more recently non-reciprocal media or using topological designs. The dielectric waveguide approach is exemplified for polarisation beam splitters, as used for Mach-Zehnder interferometers, by [3,4] and recently with more compact designs [5]. Splitters based upon dielectric waveguides traditionally have disadvantage of scale, even for high-index dielectrics requiring them to be of the order of several wavelengths in length to minimise scattering and losses at corners or shallow bends. This scale limitation, and a desire to minimise radiation losses at the bend, motivated [6] to take advantage of photonic waveguides created by removing rows or columns of the crystal array; here we design the topological analogues. These photonic waveguide devices have been successfully applied to beam-splitters and branched waveguides with T-shaped and Y-shaped branches implemented by [7][8][9], amongst many others, and the insertion of additional, or modified, array elements in the neighbourhood of the corner [10] can improve transmission or broaden frequency range; this has spawned ever more elaborate designs to optimise for various scenarios [11,12]. A related, but different, approach to creating splitters within photonic crystals is to take advantage of self-collimation [13] and self-guiding created by strong dynamic anisotropy [14]. Beam splitters, including three-way splitters [15], are created and applied to devices [16,17] using lattice array alterations to create effective mirrors or partial mirrors.
More recently, ideas from topological insulators [18] have been transposed into photonics [19], as reviewed in [20], showing promise for robust one-way edge states protected against disorder by topology. This promise is tempered by the requirement for time reversal symmetry (TRS) to be broken, and the electron spin to be mimicked by pseudospin in continuum systems. A simpler, passive and time reversal symmetric, but less robust approach is to attempt to reproduce the valley-Hall effect and utilise ideas from the field of valleytronics [21], for instance [22] create dielectric photonic topological arrangements leading to reflectionless guiding and designs for optical delay lines. Topological designs have been recently implemented for telecommunication wavelengths [23] on a CMOS-compatible chip thus bringing these concepts closer to application. These valley-Hall devices are locally topologically nontrivial however globally trivial, and therefore cannot draw upon the full power of the analogy with topological insulators, but do have advantages in terms of simplicity of construction as one need only break spatial inversion or reflectional symmetry, together with suppressing backscatter. Given the emergence of topological guiding there is now interest in developing this for photonic circuits, [24], and a natural drive to explore the potential of these new ideas. The vast majority of this valleytronics literature takes advantage of periodic hexagonal or honeycomb lattices, utilising ideas from graphene, in particular the symmetry properties of the hexagonal Brillouin zone and symmetry induced Dirac cones at the KK vertices. Perturbations of the structure, that break symmetries, then gap the Dirac points giving topologically nontrivial band-gaps and well-defined KK valleys. The valleys have opposite chirality and are related by parity and/or reflectional symmetry as well as TRS. A key point is that by engineering a large Fourier separation between the valleys one suppresses intervalley scattering and ultimately the valley is used as an information carrier [20]. One negative that emerges is geometrical: The hexagonal systems can only create two-way energy-splitters [25], and three-way splitting would require other geometries for which symmetry induced Dirac points do not occur-hence the recipe outlined above cannot be employed.
Fortunately, one can engineer accidental Dirac points [26,27] for square systems (and indeed for all two-dimensional systems, most usefully for those that do not have symmetry induced Dirac cones) and, although they are no longer at the high symmetry points one can mimic some of the effects found in hexagonal systems for the usual symmetry induced Dirac points, see for instance [28]. Recently [29] showed beam splitting for a model system of waves on elastic plates, for a fourth order elastic plate equation and masses of infinitesimal radius, here we explore to what extent three-way beam splitters can be employed in electromagnetism. We do this via a combination of group and k·p theory, [30], coupled with detailed numerical simulations to extract the interfacial edge states, zero-line modes (ZLMs), and further we take these idealised edge states and perform numerical simulations, using [31], showing transport around right-angled bends and three-way energy splitting. A typical splitter that we construct is shown in Fig. 1, it is constructed from dielectric inclusions arranged upon a square lattice, within each elementary cell there is an arrangement of inclusions and the choice of arrangement is critical for splitting; we illustrate this importance by contrasting two arrangements. A judicious choice of inclusion arrangements, and the connection of quadrants of material yield the passive splitter that, as in Fig. 1, takes an incoming wave and then splits it in three.
We operate at telecommunication wavelengths using a photonic crystal (PC) consisting of Indium Phosphide (InP) dielectric pillars in air. Drawing upon the design and fabrication of a dielectric carpet in [32], we contrast two designs of PCs within a square lattice array; the first of which has 6 InP dielectric pillars of diameter 200 nm, height 2µm, and minimum center-to-center spacing of 250nm arranged around an equilateral triangle, Fig. 2(a), and the second has 8 pillars with similar dimensions arranged around a square as in Fig. 2(b). The refractive index of InP is n = 3.16, and we consider primarily (TM) polarization whereby the electric field is parallel to the pillars' axis and discuss the analogous transverse electric (TE) polarization, the magnetic field is parallel to the pillar axis, in the concluding remarks. These two designs exemplify the potential of topological designs for energy splitting, and in particular how to achieve three-way splitting, for intrachip communication devices. Figure 1: A three-way splitter designed on a square lattice using accidental degeneracies and geometry (right panel) shows the junction region between four quadrants each containing a set of dielectric inclusions placed around a rotated square. The quadrants differ in only the rotation given to the inclusion arrangement within it. In the left panel, the wave energy incoming along the leftmost interface is split at the junction into three interface waves. In this simulation the normalized angular frequency ωd/c is 4.78 and we treat TM polarisation.
The Maxwell equations split naturally into p and s polarizations, with forcing created by an electric line source or magnetic current dipole at position r s respectively, [33], as TM and TE, with 0 , µ 0 ( r , µ r ) as the permittivity and permeability in-vacuo (and relative values), ω is frequency with time-harmonic waves, exp(−iωt), assumed, j T and I s being currents. We will often use a normalised frequency, ωd/c, hereafter where d is the pitch of the lattice and c the speed of light, c 2 = 1/µ 0 0 . The TM field is driven by a line monopole source at δ r s and the TE field by a line dipole. For the polarised fields, taking a Cartesian coordinate system , where x is the in-plane (or transverse) variable and x 3 is the out-of-plane (or longitudinal) variable, we use invariance along x 3 to split the vector Maxwell system into two polarizations: The TM polarization has E l = E(x)e 3 (TE similarly has H l = H (x)e 3 ), that is, the field is perpendicular to the pillars, and we concentrate from hereon on the TM polarized field. We take the relative permittivity to be spatially dependent, i.e. r ≡ r (x), with the remaining parameters constant, and thus for source-free TM fields we have that where ∇ 2 is Laplacian with respect to the in-plane variable x, and the material dependence is encapsulated in a(x) = µ 0 0 µ r r (x). Eq. (3) carries the assumption of material isotropy and scalar permittivity and permeability that is lossless. Moreover, for the pillar crystal a(x) is a piecewise constant function as ε r is 1 in air and 3.16 2 in the pillars, and µ r = 1 everywhere (a non-magnetic medium). This is equivalent to Helmholtz equations in each homogeneous phase, coupled through continuity conditions across the pillar interfaces of the field E and its normal derivative ∂nE = n · ∇E, with n the normal to the interface.

Engineering a non-symmetry repelled Dirac cone for a square lattice
To create three-way splitting we first need to engineer accidental degeneracies along the Brillouin zone as shown in Fig. 2(c), boundary, it is not sufficient to simply generate Dirac points, and the arrangement of inclusions in the interior of the physical elementary cell, and their symmetry, plays a key role. We utilise perturbation theory and symmetry arguments to design the structures we study; the inclusions are either in a triangular formation (and have C 3v symmetry) or are arranged along a square (and have C 4v symmetry), see Fig. 2(a,b). As we shall see the interaction of the internal inclusion symmetries with the symmetries of the lattice lead to fundamentally different behaviours; ultimately we shall see that the C 3v case is unable to create splitters for reasons uncovered by symmetry arguments. We begin by determining the criteria under which accidental Dirac points are created, then gap them using geometrical rotations that break symmetry to create band-gaps with protected edge states; these edge states form the building blocks of the splitters we design.

Perturbation theory and band interactions
We begin by considering infinite periodic media and interpret the dispersion diagrams of Fig.  3 that contrast the triangular, C 3v , and square, C 4v , cases. We take the inclusion arrangements shown in Fig. 2(a,b); when we rotate the inclusion arrangements, the rotation is taken around the centroid. The most notable difference between the non-rotated cases, Fig. 3(a) and (b), are the single Dirac point for bands 4 and 5 in (a) vis-a-vis two in (b); the most important bands are coloured red and correspond to n = 3, 4, 5, 6 in the index notation we adopt. As advertised in the introduction a symmetry breaking perturbation, in this case a rotation anti-clockwise, gaps Dirac points to create band-gaps as indicated in Fig. 4(a,b). A critical issue, of course, is where the accidental Dirac points arise, or indeed whether they arise at all, [26] and this requires an analysis of the band structure. A minor point is that, due to the precise position of the inclusions, spaced at 0.25d, the C 4v case has an additional glide and reflectional symmetry and there is a Dirac point exactly at N for bands 5 and 6, but this has no influence on the analysis here.
To understand how the bands interact we consider the eigenfunctions E nκ (x), with n the band index and κ the Bloch wavevector, i.e. the Bloch momentum vector in reciprocal space, in the first Brillouin zone such that with ω nκ the eigenfrequency. We consider Bloch waves and the eigenstates |E nκ relate to periodic eigenstates |u nκ , which form a complete basis, via with the eigenstates orthonormal, i.e. nκ |E nκ E nκ | =1, E nκ |E mκ = δ mn δ κκ . The completeness of the periodic eigenstate basis, means we can expand |u nκ in terms of the basis set {u jκ 0 (x)} where κ 0 is fixed, to deduce where ∆κ = κ − κ 0 and m running over all positive integers. Assuming the perturbation |∆κ | 1 the governing equation becomes This relationship connects the eigenstates at a given point in reciprocal space, κ 0 , to the frequency, and hence determines the local dispersion relation, and can be used to investigate band repulsion or attraction. We concentrate upon the behaviour near the point N as both sets of inclusions share σ v mirror symmetry; the C 3v case does not have σ h symmetry and this leads to the "missing" Dirac point along M X in Fig. 3(a). Looking further at the band diagrams, say that of Fig. 3(a) one observes, at the wavenumber N, close to the Dirac point of interest, that the eigenstates of bands 3 − 6 have characteristic orbitals familiar from quantum mechanics [34]. The symmetries inherent in these eigenstates motivate the use of the mathematical language of symmetry, that is, group theory: The point group symmetry of the structure is G Γ = C 2v ; this is also the point group symmetry at N, G N = C 2v ( Table 1). The C 2v point group arises from a combination of spatial (reflectional) and time-reversal symmetries; the latter relates κ → −κ. The irreducible representations (IRs) at N are one-dimensional hence there is no symmetry induced degeneracy; however, two of the IRs can be tuned such that an accidental degeneracy, that is not symmetry repelled, is created along N X (and also along M X for the C 4v case).
The four bands highlighted in red in Fig. 3(a,b) (bands numbered 3 − 6 inclusive) are associated with the eigenstates shown alongside the dispersion curves; these match the basis function symmetries of the C 2v group given in Table 1(a); hence this indicates that bands 3 − 6 are symmetry induced and the sequential ordering of them (lowest to highest) is deduced numerically, via the eigenstates, as: {B 2 , A 1 , B 1 , A 2 } and these are shown in Fig. 3. The rotated, symmetry-breaking, cases of Fig. 4(a,b) retain, at least broadly, the same eigenstate structure.
Intuitively, we expect the two bands (band 4 and 5) forming the accidental degeneracy, A 1 , B 1 , to be strongly coupled with each other; the other two symmetry induced bands, B 2 , A 2 , will have limited influence on the local curvature or slope of the A 1 , B 1 bands. The even more distant bands (bands 1 and 2 and those above 6) will have negligible effect on the A 1 , B 1 bands and we quantify this by separating the bands into the the symmetry set eigenstates (SSE, bands 3 − 6) and the remainder that lie outside the SSE. Eq. (5) becomes Multiplying (6), by the conjugated states E * pκ 0 (x) or E * qκ 0 (x) (where p ∈ SSE, q SSE) and integrating over the primitive cell in physical space we obtain the following two equations, and the H ab are explicitly, where we recall that here the function a(x) is piecewise constant. A useful point is that, using the symmetries of the eigenstates one can select which H ab are zero. Neglecting the terms coupling states outside the SSE to each other we obtain, for p ∈ SSE, that Our main interest is in the neighbourhood of the Dirac point so we set n = p ∈ SSE, and perturb with κ = κ 0 + ∆κ, then Therefore the second summation in (10), coupling states within the SSE to those outside, falls into second-order, i.e. O |∆κ | 2 , hence the first-order equation is the 4 × 4 system, where ∆ω p = ω pκ − ω pκ 0 and n, p ∈ SSE. Equation (12) contains the information that allows us to determine whether an accidental Dirac degeneracy will occur. One can take this further and notably, the higher-order corrections, that encompass the coupling between bands within the SSE to those outside, provide the band curvature details away from a locally linear point. In that instance, Eq. (12) is modified to a 4 × 4 matrix eigenvalue problem, where the Hamiltonian, with components H pm , is Hermitian [27].

Compatibility relations and creating accidental degeneracy along N X
Bands vary continuously, except possibly at accidental degeneracies where mode inversion may occur, which in turn leads to a discontinuity of the intersecting surfaces. Hence when moving along a continuous band of simple eigenvalues the eigenstates continuously transform; departing from the high symmetry point N, the associated IRs describing the transformation properties of the eigenstates smoothly transition into IRs that belong to the point groups along N Γ or N X.
In physical space both cellular structures we consider in Fig. 2 have σ v spatial symmetry, this is equivalent to σ h symmetry in Fourier space. From the definition of a point group symmetry, i.e. any symmetry operatorR ∈ G Γ that satisfies,Rκ = κ mod G, where G is a reciprocal lattice basis vector, implies that κ ∈ N X solely has the mirror symmetry operator, σ h within its point group. Similarly, for a κ ∈ ΓN, only the vertical mirror symmetry operator, σ v satisfies the point group criterion. The symmetries of the eigenstates, for a κ belonging to either of the paths, N X and N Γ, are shown within the basis functions column of Table 1 We now simplify the 4 × 4 system of Eq. (12) by assuming it is dominated by the two strongly coupled bands, bands 4 and 5 with IRs A 1 and B 1 , which will result in a 2 × 2 system. As we move away from N, the A 1 , B 1 IRs belonging to the C 2v table, from continuity of the bands, transform into the IRs of the σ v,h table and reveal compatibility relations of the IRs. For the symmetry σ h the eigenstates at N and along N X satisfy the following, whereP is the projection operator. The bands (A 1 , B 1 ) at N are therefore compatible with (A, B) along N X. Physically, this transition is also evident from the eigenstates. Similarly, at N and along N Γ, the eigenstates transform under σ v as, The bands (A 1 , B 1 ) at N are therefore compatible with (A, A) along N Γ which implies band repulsion and an inability to have a band crossing along N Γ as we observe. Importantly, note that, in deriving equation (12) we have only assumed that κ 0 belongs to a particular symmetry set band (surfaces 3 − 6) (the band at κ 0 must be continuously connected to the same band at N). Therefore, the compatibility relations allow us to choose any expansion point along the the path ΓN X where the eigenfunction basis set, (5), transforms accordingly i.e.
In order to solve the 2-band eigenvalue problem, Eq. (12), we compute the determinant of the truncated matrix, where parity simplifies the Hermitian matrix; the eigenstates are evaluated at κ 0 . Solving the eigenvalue problem yields the following result, where the ± corresponds to the A 1 , B 1 bands, respectively. This result implies that the A, B bands have an identical slope, albeit with opposite gradients; hence, if, at N an instance can be found where ω B 1 > ω A 1 then the bands will invariably cross along the path N X. We are not guaranteed that an accidental degeneracy must occur along N X as parameters (the radii of inclusions, number of inclusions, permittivity etc) could occur with ω B 1 < ω A 1 , but this inequality gives a useful criteria for their existence, or otherwise. This parametric freedom afforded by inclusion change in geometry or material tunes increases, or decreases, in the slope thereby increasing or decreasing the distance between N and the Dirac point. Note that the Dirac cone occurs along the spatial symmetry path, σ h , of the structure due to the opposite parities of the A, B bands; band repulsion occurs along the N Γ path [27] thereby resulting in a partial band gap along N Γ. If ω B 1 > ω A 1 , then the partial gap along N Γ isolates the Dirac cone along a portion of the IBZ path, ΓN X.
We have designed situations where either two, or four, (for C 3v and C 4v respectively) pairs of Dirac points are created by an accidental degeneracy and we have shown how they are created and given a prescription for their occurence. These Dirac points have been gapped by a symmetry breaking perturbation and band gaps have been created.

Edge states: Zero line modes
We aim to use our knowledge of the designed band-gaps to create situations whereby an interface will support edge states. We begin by taking a half-space of one medium and place it above another; the only difference between the two media being that the symmetry of the cells in the upper and lower media is broken by clockwise and counter-clockwise rotations respectively (see Fig. 5(a,b) for C 3v ). This has the effect, in Fourier space, of interchanging N and M, the Berry curvatures [35] are opposite in sign at N and M, and hence at the interface between such media and create valley Hall edge states; these are named zero line modes (ZLMs) due to their origin as arising from opposite Berry curvature in the adjoining media.
We compute ZLMs by considering finite tall ribbons one cell thick with Bloch conditions  the polarisation chosen, the even modes are the physically relevant ones. The eigenstates and dispersion curves are shown in in Fig. 6 and the parity of the ribbon eigenstates is inherited from the bulk eigenstates of Fig. 3, 4. Notably the C 3v and C 4v cases have a crucial difference: ordering of the media matters, that is, stacking the clockwise above the counter-clockwise or vice-versa leads to two different interface types for the triangular case whilst for the square case it is irrelevant which ordering is taken (see Fig. 5). The differences in the interfaces in the C 3v versus C 4v cases form a key distinction that impacts upon energy transport around corners and splitting. Full scattering simulations performed using the commercial finite element package COMSOL [31] are shown in Fig. 7, these are for the even ZLM relevant for TM polarisation and for the C 3v case; very similar ZLMs are found for C 4v (not shown). The excitation is a line source just outside the leftmost edge of interface and a very clear ZLM is excited that can be identified, on each ribbon, with the eigenstate from the ribbon. One feature that is also evident is the long-scale wave envelope with a wavelength that alters with frequency; this can be described via an effective medium approach [36] as applied to edge states [29].

Energy transport around a sharp bend and energy splitting
We now consider two related, yet distinct, problems: redirecting energy around a sharp bend using just topology, see Fig. 8(a), which is the topological alternative to the photonic crystal waveguides pioneered in [6], and a three-way energy splitter (see Fig. 8(b)). Both panels in Fig. 8 are for inclusions placed around a square, the C 4v case; the C 3v case, despite creating a clear ZLM as in Fig. 7, is incapable of supporting a ZLM along the vertical interface as those interfaces, c.f. Fig. 5, do not have non-zero Berry curvature. The square arrangement of inclusions has the very useful property that the vertical interfaces are exactly the same as the horizontal ones; the added benefit of having two reflectional symmetries in the C 4v case is now evident there is now non-zero Berry curvature along both interfaces and both support the same ZLMs. It is this insight that allows for the design of the splitter and allows transport around the bend.
We now want to optimise the transport properties. First, we may wish to minimise backscatter from the junction. To do so we note that the Fourier separation between the Dirac point of the unperturbed bulk dispersion curves and high-symmetry point is highly relevant for the transmission properties of the topological guide [25]; transmission improves for short wavelength, as opposed to long wavelength, envelopes, hence, for transmission post the junction, it is desirable to increase the distance between the Dirac cone and the point N. The latter holds due to the connection between the bulk and projected bandstructures [37]; the Brillouin zone reduces to one-dimension because the only relevant wavevector component is along the interface. All wavevectors are projected onto the ΓN line in Fourier space, hence if the distance between N and the Dirac cone is increased then the Fourier separation between oppositely propagating modes, along the topological guide, is increased. A mechanism to do this is to alter the system parameters; Eq. (16) demonstrates that the slopes of the A and B bands can be increased or decreased by the system parameters thereby altering the position of the band intersection. Second, we return to envelope wavelength and note that the distance from the source to the junction will play a role. If we took a finite length slab then fitting an integer, or half integer, numbers of envelope wavelengths along the lead interface gives a Fabry-Perot resonance with perfect transmission or perfect reflection from the far edge. We can use this knowledge to tune the system and the sharp bend is optimised by having a node of the long-scale envelope at the junction and so the energy is smoothly transported around the corner. Whilst for the splitter, the perfect reflection scenario concentrates energy at the junction for subsequent redistribution to the exit leads. This is the first example of a three-way splitter passively created due to its topology and as a result of the inherited protection should be less prone to backscatter and therefore forms the prime candidate for three-way splitting in time reversal systems.

Concluding remarks
We have constructed a three-way splitter, and designed for propagation around a right angle bend, using valleytronics and the presence of accidental Dirac points. Although we have concentrated upon the TM polarisation, it is clear that TE polarisation will also generate splitters using the geometrical designs here; the main change being that the odd (and not even) modes, now excited by dipoles, would be the physically relevant fields. We briefly illustrate this in Fig. 9 for a square C 4v case; here the inclusions of Fig. 2(b) are augmented by a central inclusion of radius 0.15d. The situation is almost identical to TM, except that the bandgap is smaller (the action of the central inclusion is to help create the bandgap) and the decay of the edgestate, shown in Fig. 9(c), is slower.
One crucial difference from the majority of the valleytronics literature is that we have chosen to operate on a square, and not hexagonal, lattice; the hexagon arrangement has several advantages -the Dirac points are symmetry induced, and the band-gaps obtained by gapping them can be constructed to be broad. For energy transport from one interface to another both the chirality and phase must match; this is easily achieved for the square arrangement designed here as the interfaces match. For four joined quadrants of structured hexagonal lattice there is a mismatch in phase and so you are confined to two-way splitting and thus the square splitter here appears the only option.