Topological crystalline protection in a photonic system

Topological crystalline insulators are a class of materials with a bulk energy gap and edge or surface modes, which are protected by crystalline symmetry, at their boundaries. They have been realized in electronic systems: in particular, in SnTe. In this work, we propose a mechanism to realize photonic boundary states topologically protected by crystalline symmetry. We map this one-dimensional system to a two-dimensional lattice model with opposite magnetic fields, as well as opposite Chern numbers in its even and odd mirror parity subspaces, thus corresponding to a topological mirror insulator. Furthermore, we test how sensitive and robust edge modes depend on their mirror parity by performing time dependent evolution simulation of edge modes in a photonic setting with realistic experimental parameters.

Introduction -Symmetry and topology are two fundamental mathematical tools in the classification of states of matter in condensed matter physics. Recently, intensive research interests have been focused on the role of symmetry in the classification of topological states, ever since the discovery of time reversal invariant topological insulators 1-9 , in which time reversal plays an essential role. A large variety of symmetry protected topological states have been identified theoretically for different symmetry classes and dimensions 10,11 In contrast, the corresponding material realization in experimentally feasible systems of these new topological states has only been limited in several symmetry classes, mainly for time reversal invariant topological insulators [12][13][14][15][16][17][18][19][20] , the quantum anomalous Hall insulators 21,22 , topological superconductors [23][24][25][26][27][28] and topological mirror insulators [29][30][31][32][33][34] . Therefore, searching for new topological systems for symmetry protected topological states is vital for the development of this field.
Topological crystalline insulator (TCI) phases 29 , topological phases that are protected by crystalline symmetry, can exist in a large number of crystal structures with different space group symmetries and a classification of TCIs in different point symmetry groups and space symmetry groups have recently been the focus of research [31][32][33][34][35][36][37][38] . Compared to quantum spin-Hall-type topological insulators, which can only occur for spinful fermions, topological crystalline insulators can exist in both fermionic and bosonic systems. This fact opens up the possibility to realize TCIs in various bosonic systems, including photonic, phononic, magnonic and cold atom systems. As we describe below, photonic systems in particular have generated a great deal of interest as a probe for topological physics.
The realization by Haldane and Raghu 39,40 that gyromagnetic photonic crystals could exhibit non-trivial topological invariants opened the door to the new field of topological photonics 41 in which the propagation of photons in a dielectric structure is protected in a similar sense to electrons in a crystal lattice. The first experimental realization of this phenomenon was made in the group of Soljačić 42,43 for the microwave regime. However, scaling the wavelength down to the optical regime (in order to realize topological states in optical devices) was not possible using this mechanism due to weak magnetic response in that frequency regime. Other mechanisms were proposed [44][45][46][47] , and finally experimental demonstrations were made 48 in a system based on evanescentlycoupled helical waveguides; as well as in two-dimensional coupled ring resonators 49 . While photonic topological protection is conceptually similar to that of electronic topological protection (after all, this phenomenon comes down to non-interacting wave dynamics), photonics offers unique advantages and potential applications. To name a few examples, photonic systems can be designed directly by fabrication (allowing any desired lattice structure to be realized); experiments can be carried out at room temperature, meaning any emergent devices can be brought to application more realistically than those that require very low temperatures; and the robustness associated with topological protection could be of use in an array of devices that rely on the flow of light (e.g., sensors, optical interconnects, electrooptic modulators, isolators, among others).
Very recently, a prediction was made that topological photonic crystals with surface states protected by the glide symmetry could be realized for microwave photons 50 in a macroscopic ferrimagnetic structure. In this work, we propose the optical realization of crystalline symmetry protected boundary modes in photonic crystals phase using an entirely distinct mechanism from the microwave. The systems is quasi-one-dimensional, where an auxiliary parameter, φ, of which the Hamiltonian is a function, is used in place of a second spatial dimension. This is reminiscent of Ref. 51, in which a family of one-dimensional systems (defined by a pump parameter) were used to realize a photonic topological edge state and pump that are mapped to the integer quantum Hall effect. We note that another work 52 has predicted an analogue of a quantum spin Hall system that requires C 6 crystalline symmetry; however the edge always breaks that symmetry and so the topological edge states are not protected (and are gapped). The realization of crystalline symmetry protected topological states in the optical regime allows for a novel paradigm in exploring topological crystalline phases, opening up the possibility to explore the relationship between topological photonic crystals and nonlinear/interacting bosonic systems and quantum optical effects such as multi-photon quantum walks.
Tight-binding model -We start from a simple tightbinding model to illustrate our main idea and then simulate the system in a more realistic situation. We consider a quasi-one-dimensional (1D) chain along the x direction with each unit cell consisting of four sites, denoted by α = 1, · · · , 4, as marked by "a" in Fig. 1A. The positions of four sites are chosen to preserve y-directional mirror symmetry with respect to the line denoted by "b" in Fig. 1A. Bound states can be induced by lowering on-site energy, forming a potential well on each site. We denote these states as |s α , and the corresponding creation and annihilation operators as c † α and c α , respectively. Thus, the tight-binding Hamiltonian for this system is given by where α, β = 1, · · · , 4 and m = 1, 2, · · · , M denoting the unit cell index.V m describes the on-site energy and the hopping between sites within the m-th unit cell.T m is for the hopping between two nearest neighbor unit cells m and m + 1. In the {|s α } basis, the HamiltoniansV m andT m are given bŷ (2) where V m (φ) is the on-site energy of each site in the m-th unit cell, chosen to be identical. γ i denotes the hopping strength along the i = x, y direction and t, t ′ are hopping parameters between two unit cells, with γ i , t, t ′ < 0. Only vertical and horizonal nearest-neighbor(NN) intrasite hoppings, denoted by "c" in Fig. 1A, as well as horizonal and diagonal NN inter-site hoppings, denoted by "d" in Fig. 1A, are considered.
The above Hamiltonian preserves the y-direction mirror operatorM y , given bŷ which interchanges the sites 1 and 4 (2 and 3). We set Let us neglect t, t ′ , η m and δ m , and set γ x0 = γ y0 = γ 0 for the moment. In this limit, the HamiltonianV m in (Eq. 2) can be diagonalized and the eigen-energies of four eigenstates are shown in Fig. 2. We find that two degenerate states |ψ + and |ψ − with eigenenergy V 0 are well separated from other two states |ψ t and |ψ b with energies V 0 ∓2γ 0 , where we use + (−) to represent odd (even) mirror parity ofM y . Thus, in the limit η m , δ m , t, t ′ ≪ V 0 , γ 0 , we can focus on two degenerate eigenstates |ψ + and |ψ − , and project the total Hamiltonian onto the subspace spanned by these two states. The eigenstates of two degenerate states are given by |ψ + = 1 2 (+ |s 1 + |s 2 − |s 3 − |s 4 ), |ψ − = 1 2 (+ |s 1 − |s 2 − |s 3 + |s 4 ), and let us denote d ±,m and d † ±,m to be the annihilation and creation operators for the states |ψ ± on the m-th unit cell, respectively. The effective Hamiltonian after projection is given by where the summation ± is over the parity subspaces. The parameters t ±,m and V ±,m are related to the original parameters t, t ′ , δ m and V m by V ±,m = V 0 + η m ∓ δ m and As a result, the Hamiltonian (Eq. 4) is nothing but two copies of Aubry-Andre-Harper (AAH) Hamiltonian 51,53,54 with effective flux strength b and −b in the odd and even subspaces, respectively. The parameter φ is also known as the pump parameter. Here m 0 shifts the index of the first unit cell, introducing a degree of freedom which will be discussed later. As a consequence of opposite mirror parities underM y , these two copies of Hamiltonian are decoupled from each other once mirror symmetryM y is preserved. It has been shown that by correctly choosing parameters, boundary modes can exist at the end of a finite chain described by the AAH model, within fractal sets of band gaps 55 . Thus, we perform a direct calculation of energy spectrum of a finite chain for the model (Eq. 1). Here for demonstration, we choose the parameters to be 5. An extra offset between γ x and γ y is set by taking γ x0 = γ 0 + 0.025 and γ y0 = γ 0 . This offset shifts the relative energy level for odd and even bands, in order to emphasize the crossings of edge states inside the mini-gap described later. The four energy levels in one unit cell |ψ + , |ψ − , |ψ t and |ψ b expand into four energy bands as a function of the parameter φ, as shown in Here we show the zoomed-in spectrum in Fig. 3B of the mini-gap region "a". We find four energy levels in the mini-gap, denoted by |ψ L+ , |ψ R− , |ψ R+ , and |ψ L− , or "b", "c", "d" and "e" in Fig. 3B, respectively, where L and R denote the position where the wave function of the state is localized in the chain, and ± indicates the parity of the boundary modes, shown by blue (odd) and red (even) in Fig. 3B. A more realistic calculation of the wave function in this system will be given in the next section, based on simulating the full continuum problem. The crossing between the boundary modes |ψ R− and |ψ R+ , marked by "f " in Fig. 3B, is topologically protected due to opposite mirror parities between them, while the crossing between |ψ L− and |ψ R− , marked by "g" in Fig. 3B, is gapless because these two states are located at opposite boundaries of the chain. As long as the chain size is large compared with the penetration length of boundary mode, which is around twice the lattice constant for the parameters listed above, the overlapping between wavefunctions at opposite boundaries is negligible.
To test the symmetry protection by adding a mirrorsymmetry breaking term ∆H ontoV m , where and ∆ = 0.04, we find a gap opening between two edge modes ψ R− ("c") and ψ R+ ("d"), as marked by "h" in Fig. 3C. Thus, we conclude that these boundary modes are stable only when mirror symmetry is present in the system (Eq. 1).
We would like to emphasize that although our model is written in one dimension, it can be mapped to a twodimensional (2D) lattice model with complex hopping terms, corresponding to a 2D topological mirror insulator. This is a similar mapping as that made in Ref. 51. We may introduce a fictitious dimension w and extend our 1D tight-binding model to the x − w plane, forming a 2D square lattice as shown in Fig. 1C. Similar to the unit cell defined in the 1D model, each unit cell at the integer coordinates (x, w) consists of 4 sites, each having one state described by the creation operators c α † x,w where α = 1, . . . , 4. The 2D Hamitlonian can be written as whereT x(w) is the hopping matrix between sites labeled by α, β = 1, 2, 3, 4, in the adjacent unit cells along the x (w) direction, andÛ is the hopping matrix among sites in one unit cell. The detailed forms ofT x(w) andÛ are given in the Appendix. This Hamiltonian is also invariant under the mirror op-erationM y (Eq. 3) about the x-axis. We define the basis creation operators d † ± in the same way as in Eq. 4. where ± are defined for the subspaces with odd and even mirror parities as before. In this new basis, the Hamiltonian a.u. Boundary modes |ψL+ , |ψR− , |ψR+ , and |ψL− , as indicated by "b", "c", "d" and "e", are found within the mini-gap "a". "f " ("g") denotes the crossing between "c" and "d" ("c" and "e"). The light red and light blue shading indicate the mini-gap region "a" of two corresponding parities. (C) A gap opening for the crossing "f " can be induced by a mirror-symmetry-breaking term ∆H, leading to an anti-crossing "h".
takes the form where the parameter θ = 2πbx represents the phase shift during the hopping along the w direction. We notice that the form of θ corresponds to the Landau gauge for a magnetic field with the flux b in one unit cell. Thus, a Chern number can be defined in each mirror subspace for a non-zero b, as described in the appendix. The opposite signs for the phase shift in the mirror even and odd subspaces indicate that the total Chern number is canceled for the whole system. However, a mirror Chern number, as defined in Appendix, can characterize non-trivial topological property of our 2D system 32,51 . Therefore, this mapping suggests that our tight-binding model provides a realization of 2D topological mirror insulators. Numerical Simulations for photonic systems -We have now established topological mirror insulator phases in our simple 1D AAH type of tight-binding model and its relation to 2D topological mirror insulators is also illustrated above. The realization of this model in a photonic system requires more sophisticated numerical simulation, which includes continuum degrees of freedom and are based in Maxwells' equations. Next, we describe the detailed experimental setup and perform a numerical simulation of photonic lattice systems.
To realize the tight-binding model in a realistic photonic system, we consider an array of evanescently coupled elliptical waveguides in fused silica glass. The difference in refractive index inside and outside the waveguide is utilized to confine the light in the x − y plane and serve as a potential well, as shown in Fig. 1B. We follow the standard paraxial approximation for this type of waveguides, in which the Maxwells' equation of light can be simplified to a Schrödinger-type equation, namely: where ψ is the envelope function of electric field, k 0 is the wavenumber of ambient light in the medium, n 0 and ∆n are respectively the background and the waveguide deviation from background refraction index. In comparison with the Schrödinger equation, one can see that the z direction can be regarded as time, therefore the diffraction of light through the waveguides is equivalent to the time evolution of a particle in a potential determined by −∆n(x, y; z). This equation has been utilized to describe other topological phases in photonic lattices 48 , as well as a wide variety of effects in linear and nonlinear optics 56 , including (but not limited to) the prediction 57 and observation 58,59 of lattice solitons, stable photorefractive solitons; Anderson localization in optics 60 ; conical diffraction 61 ; and optical pseudomagnetism 62 . The waveguides used in this work can potentially be fabricated using the femtosecond direct laser writing technique; and are inspired by those described in Ref. 63. As shown in Fig. 1B, four elliptical waveguides reside on the corners of a rectangle to simulate one unit cell in our tight-binding model. We place these unit cells one by one to form a chain, as described in the tight-binding model. The confining potential is provided by the deviation of the refraction index inside the waveguides from the background, V = (ω/c)∆n, and the hopping term is determined by the overlap between two nearby sites, which can be controlled by adjusting the distance between them. Details of the relation between distance of different sites and energy level positions are discussed in the Appendix.
For this waveguide chain with a finite number of unit cells (M = 23), we first discuss our numerical simulation of the energy spectrum, which is shown in Fig. 4A and C for two waveguide chains with different parameter sets (See appendix for details of the choice of parameters). For simplicity, we label the waveguide chain for the energy spectrum in Fig. 4A and C as numbers I and II, respectively. Here we only focus on energy bands that correspond to |ψ ± bands in our tight-binding model. The energy spectrum of our continuum also reveals a sub-band structure with mini-band gaps for both parity bands and a detailed comparison between continuum model and tight-binding model is shown in the Appendix.
In short, we find that the continuum and tight-binding models strongly agree. For both the bands with even (red color) and odd (blue color) parities, they are split into three mini-bands, separated by mini-gaps. Within the mini-gap, boudary modes can exist for both parity bands. Examples of boudary modes are marked by a for even parity bands in Fig. 4A (waveguide configuration I) and b for odd parity bands in Fig. 4C (waveguide configuration II), both at the left edge of the chain. For these two states, the corresponding wave functions are shown in Fig. 4B and D, respectively, from which one can clearly see that both modes are highly localized with a penetration length around twice of the lattice constant. We emphasize that we need to carefully choose the value of φ in Fig. 4A (Fig. 4C), and at that φ value, only one boudary mode "a" ("b") at that edge is supported and this is important for the experiments of parity selected time evolution of edge modes described below.
We next inject a wave packet into one end of the waveguide chain and study time evolution of wave packets with different parities in this chain by solving the time-dependent Schrödinger equation, Eq. 8, numerically. Technique details of our simulation are included in the appendix. The initial wave packet is chosen to be in Gaussian form, given by where N is a normalization factor, (x 0 , y 0 ) is the center of Gaussian form and chosen to be the center of the first two unit cells. The function f describes the parity of the injected light beam and chosen to be f − (x, y) = e ikx(x−x0) for even parity mode and f + (x, y) = sin(k y (y − y 0 )) for odd parity mode. The size of the envelope function is tuned to cover the majority of the mode wavefunction. We inject an even beam ψ ini − (x, y) into the waveguide chain I with φ = 1.15π (the value of φ for the mode "a" in Fig. 4A), and the corresponding time evolution is shown in the first row of Fig. 5, from which one can see that the wave packet keeps localized after a long time evolution. In contrast, when we inject ψ ini + (x, y) into the same waveguide configuration, the light spreads into the bulk of the waveguide chain, as shown in the second row of Fig. 5. For chain II with φ = 0.25π (the value of φ for the mode "b" in Fig. 4C), the light will get delocalized with initial wave packet ψ ini − (x, y), but localized with initial wave packet ψ ini + (x, y), as shown in the last two rows of Fig. 5. Therefore, we conclude that the localization of edge modes at the boundary sensitively depends on their mirror parities, and our numerical simulation also provides an approach to probe topological crystalline protection of boundary modes in realistic experiments.

I. DISCUSSION AND CONCLUSION
In conclusion, based on a tight-binding model and continuum numerical simulation, we propose an experimental setup of a photonic lattice to realize topological crystalline protection of boundary modes in photonic systems. Numerical simulation also suggests a possible experimental configuration to detect topological edge modes, of which the edge localization behavior depends on the parities of both the injected wave packets and localized edge modes. In addition, we would like to emphasize that two degenerate states realized in our system can also be regarded as "pseudo-spin" and thus provide a natural platform to construct SU(2) Landau levels 64 and spin-orbit coupling in future studies. The fact that the model was simulated with both tight-binding and more experimentally realistic continuum simulations suggests that can be implemented in a scheme similar to previously realized time-reversal-broken topological phases in photonic lattices. This will open the door to the exploration of topological crystalline protection in interacting systems based on the nonlinear optical response of the ambient medium (giving rise to the nonlinear Schrödinger / Gross-Pitaevskii equation), and entangled quantum walks in such phases (based on injecting entangled photons into the structure). A question of central importance in TCI physics is the question of disorder: since edge state protection is achieved using a symmetry that is easily broken (i.e., crystalline/mirror symmetry) by disorder, what is the nature of the robustness of these states? Will all protection simply break down, or is it preserved in an ensemble-sense? These are fundamental questions where photonics provides a unique and versatile path forwards. In this section of the appendix, we will describe the detailed form of 2D topological mirror insulator model, which is mapped from the 1D topological mirror insulator model in the main text. The 2D tight-binding Hamiltonian is given by satisfying 2γ 0 >> λ. Here γ 0 plays a similar role as that in the 1D model. γ 0 ensures the ψ t and ψ b states are away from the ψ ± bands. Here α, β = 1, · · · , 4 are indices for four sites in one unit cell. I is the identity matrix. We will show that θ is linked to a virtual external gauge field, and the parameters t and t ′ here match the definition in Eq. 2.
We can project the Hamiltonian onto the subspace of d ± basis, which is written as where ± is for odd and even subspaces, respectively. The parameter θ could be interpreted as the phase shift due to an external gauge field, where the field directions for the odd and even subspaces are opposite. We may choose the phase factor as ±θ = ±2πbx for two parity subspaces, which correponds to the Landau gauge of a magnetic field in the x − w plane in both the even and odd subspaces. Substituting the Fourier transform along the w direc- We find that H(φ) takes the same form of Eq. 4. From the Fourier transform, the phase parameter φ in the original model plays a role of the momentum in the extra dimension w.

Appendix B: Numerical Simulation
In this section, we will describe the details of our numerical simulation of time dependent Schrödinger equation Eq. 8, as well as the parameter setup of our waveguide chains I and II. In the numerical simulation, the lattice constant is set to 28 µm. The separations between center of sites in the same unit cell is roughly 12 µm in x and 13 µm in y direction. The other system parameters are the wavelength (633nm); the change in refractive index, ∆n = 7.2 × 10 −4 ; and n 0 = 1.45 (for silica). The system is discretized into a grid of 60 × 120 pixels for each unit cell on average.
The shape of the potential well for each site is given by where σ ′ x = 1.9 µm, σ ′ y = 5.5 µm. The parameter V 1 = (ω/c)∆n determines the depth of the potential well. Our numerical test shows that the on-site energy is linearly proportional to this depth within the energy range that we are interested in. The overlap between the wave functions for neighboring sites, which contributes to the hopping parameter, drops exponentially as a function of distance, as shown in Fig. 6. From this numerical test, we find the hopping between neighboring sites γ dominates over the hopping along the off-diagonal direction γ diag (shown in the inset of Fig. 6), and thus we can neglect the off-diagonal hopping term between the sites 1 and 3 (2 and 4) and obtain the Hamiltonian Eq. 2. By controlling carefully the distance between different sites, four elliptical waveguides in one unit cell can indeed reproduce four eigen states with the desired energy levels as in our tight-binding model. By interpolation on Fig. 6, we can establish the relation between the parameters described in the tight-binding model and the realistic parameters, such as the position of the potential well and refraction index in each site, in the continuum simulation model. Thus we are able to map the configuration in numerical simulation to the parameters in our tight-binding model. The parameters of the two models fall in the same range if we take the unit of hopping parameters in the tight-binding model to be mm −1 .
If we focus on the structure ofV m in Eq. 2 with the condition γ x = γ y = γ 0 , the commutation relationship of V m (φ),V m (φ + π) = 0 suggests that there will be multiple boundary states within the band expanded by |ψ + and |ψ − . For an arbitrary φ, there could exist more than  Fig. 7. The values are fitted from the mapping between the tight-binding model and the continuum simulation.
one boundary mode with different energies, as shown in Fig. 7. A similar argument can be applied to numerical simulations and results in the excitation of multiple modes when the input beam is injected on the edge. To find a specific φ where there is only one boundary mode occupying a boundary, an integer m 0 setting the index of the leftmost unit cell is introduced into the tight-binding model and it does not change the topological properties of the system 51 . By carefully tuning m 0 , we are able to find the configurations with only one boundary mode of one mirror parity.
With the mapping of parameters between tightbinding model and continuum simulations, we can also identify the required wave-guide chain configurations for our continuum simulations. The dispersion spectra for the chains I and II have been given for continuum model in Fig. 4A and C. As a comparison, we also show the spectra of the chains I and II in tight-binding model in Fig. 7. The positions of boundary modes for both parities match with the continuum simulations, although band widths and energies reveal minor difference. This is due to neglecting the long-range hopping and the correction of t from the change of separations within nearby sites. Since we are only interested in boundary modes, this discrepancy is not important for our purpose. The detailed parameters used for the tight-binding model is listed in Table I.
The normalized input beam Eq. 9 is treated as the initial wave function at time=0. For each timestep dt in the simulation, the real-space wave function is multiplied by the time evolution operator e −iHdt , where the potential part of Hamiltonian H is discretized from the sum of the cubic-Gaussian potential well functions Eq. B1 and the Fourier transform of the kinetic part. The reflection of light on the boundary of the grid is manually reduced to emulate experimental conditions. The timestep is set to 5um/(c/n 0 ). defined by 51,66 where C p (φ, ξ) = Tr P p ∂P p ∂φ , ∂P p ∂ξ . (C4) The above expression is equivalent to the famious Thouless-Kohmoto-Nightingale-Nijs formula 67 . With this definition of Chern number, the mirror Chern number can be defined as which characterizes topological mirror insulator phase for a system with mirror symmetry.