Quantifying the robustness of topological slow light

Low-dimensional nanostructured materials can guide light propagating with very low group velocity vg. However, this slow light is significantly sensitive to unwanted imperfections in the critical dimensions of the nanostructure. The backscattering mean free path, xi, the average ballistic propagation length along the waveguide, quantifies the robustness of slow light against this type of structural disorder. This figure of merit determines the crossover between acceptable slow-light transmission affected by minimal scattering losses and a strong backscattering-induced destructive interference when xi exceeds the waveguide length L. Here, we calculate the backscattering mean free path for a topological photonic waveguide for a specific and determined amount of disorder and, equally relevant, for a fixed value of the group index ng which is the slowdown factor of the group velocity with respect to the speed of light in vacuum. These two figures of merit, xi and ng, should be taken into account when quantifying the robustness of topological and conventional (non-topological) slow-light transport at the nanoscale. Otherwise, any claim on a better performance of topological guided light over conventional one is not justified.

Slowing the speed of a light pulse down to human pace (m/s) requires complex interference effects [1] which manifest as a flat dispersion relation ν = ν(k), where k is the conserved wave vector and ν(k) the frequency. The group velocity vg of this slow light is determined by the derivative of the flat band and the slowdown factor is given by ng = c/vg, where c is the speed of light in vacuum. ng is the figure of merit for slow light and it determines the enhancement factor for diverse applications such as optical nonlinearities [2], optical switching [3], pulse delay [4], quantum optics [5], optical storage [6] and optical gain [7]. Slow light interacts much more efficiently with matter, e.g., atoms [1], photon polaritons [8], plasmons [9], spinors [10], magnons [11] or excitons [12]. A strategy to bring slow light to the nanoscale exploits optical resonances built up by nanostructuring a dielectric or semiconductor with low absorption, such as silicon at telecom wavelengths [13][14][15]. This potentially disruptive technological platform for enhanced light-matter interaction enables photonic applications ranging from nanolasing [16] and energy harvesting [17] to integrated quantum photonics [18]. Flat bands arise naturally in these systems based on the periodic modulation of the refractive index at optical or near infrared wavelengths [19] for which the group index diverges as ng ∝ (∂ν/∂k) −1 in the ideal situation. However, in real devices there is a limitation to the maximum ng achievable due to slight deviations of the fabricated parameters compared to the designed values. Even fluctuations in the nm-range [21] give rise to backscattering of the guided light, inducing a strong interference [23], a photonic manifestation of Anderson localization in low dimensions [24]. Imperfection limits slow light in conventional photonic waveguides to maximum values around ng ≈ 100 for very short waveguides with lengths L 5 µm, much lower than the ng values observed in atomic systems [1] but still sufficiently large to explore weak light-matter interaction leading to cavity-quantum electrodynamic phenomena [5]. In this Letter, we confirm that this limitation may be overcome by exploiting photonic topological effects that purely arise from engineering the lattice geometry.
Topological photonics has emerged very recently as a competitive approach for robust light transport [25], something extremely appealing for technological applications. There is a broad spectrum of possible implementations of photonic topological insulators which translate concepts from solid state to optical systems. The most robust one implements timereversal symmetry broken edge states [26] for which the timereversed state at -k does not exist, thus preventing direct elastic backscattering of the propagating mode. However, to implement this approach in photonics one requires strong mag-  [4]. (b) Distribution of the dielectric function in the structure: black corresponds to silicon and gray to air. Each valley crystal is formed as a triangular periodic lattice with a unit cell composed by two pillars with different diameters to break the spatial inversion symmetry. (c) Dispersion relation, ν = ν(k), and (d) calculated group index of the interface topological edge state. The unit cell is detiled in the inset. For reference, we design a conventional photonic waveguide obtained by leaving a row of pillars from a triangular lattice (f ) with a desertion relation and a group index plotted in (g and (h)), respectively. The black point in the group index curves denotes the frequency at which each waveguide has a group index of ng ≈ 300. The electromagnetic field intensities calculated at these frequencies and for perfect (non-disordered) waveguides are plotted in (e) and (i).
netic effects which at visible and near-infrared wavelengths are very weak and challenging to implement [27]. In timeinvariant topological insulators based on the quantum-spin Hall effect [28] and the valley-Hall effect [4], topology emanates from the breaking of particular spatial symmetries. In such imple- (c) and (d) show the ensemble-averaged electromagnetic field intensity excited by the dipole source at the same frequencies of a and b. For this calculation, over twenty configurations of positional disorder of the pillars with σ = 0.001a were averaged. (e) and (f ) illustrate the calculated ensemble-averaged electromagnetic field intensity profile along the topological and conventional waveguide axis, respectively. The localization length is extracted from the inverse of the exponential decay slope.
mentations, reciprocity imposes the existence of the counter propagating mode at -k. This time-reversed edge state carries the opposite value of a binary degree of freedom that plays the role of a pseudo-spin [28]. In this case, the key open question is whether backscattering is reduced and the answer will depend on whether the existing structural disorder preserves the pseudo-spin value or not. Another rather interesting question is whether the symmetries of the topological system which determine its topological invariants are preserved under different types of perturbations. Recent ground breaking experiments reported robustness in terms of a certain lack of structural back-reflection when precisely-shaped local defects were introduced in different topological waveguides [30]. However, this claimed robustness still needs to be systematically quantified against disorder, which is crucial to compare properly these waveguides to the state-of-the art conventional (nontopological) ones. Here, we engineer slow light in a valley-Hall waveguide [4] and we calculate its backscattering length, ξ, versus disorder and ng. Finally, we compare the results to those of a conventional photonic waveguide.
For our analysis, we focus on the parity-symmetry breaking valley-Hall approach [4]. When applied to photonic slabs, the topological edge states at an interface between valley Hall crystals of opposite K-valley pseudospin lay below the light line of the slab and are decoupled from the radiation continuum. In principle, these edge states have no intrinsic out-of plane losses and only fabrication imperfection can induce coupling to radiating modes. Even in that situation, in plane back-scattering is largely the dominant loss mechanism at large values of ng [5] so it is enough to consider the system as two dimensional to capture the physics of slow light backscattering, something not possible with other implementations of topological photonics [32]. Therefore, we implement two dimensional simulations instead of the more computationally expensive threedimensional slab. Fig. 1(a) displays an illustration of a fully two-dimensional valley-Hall waveguide formed at the interface of two topologically different valley crystals with bandgaps induced by breaking spatial inversion symmetry. The valley crystals are created with a triangular lattice where the unit cell is formed by two circular silicon pillars surrounded by air with different diameters d1 = 0.4a and d2 = 0.2a, where a is the lattice constant. The section of the waveguide is plotted in Fig. 1(b). The dielectric constants for silicon and air are taken here as εSi = 12 and εair = 1, respectively. For d1 = d2, the system preserves the C6v symmetry and supports a symmetryprotected gapless band structure between the first and second lowest energy bands for transverse magnetic polarized light. When the pillars have a different diameter, the system breaks the spatial inversion symmetry opening a bandgap between these two bands [4]. Fig. 1(c) shows the calculated dispersion relation of the waveguide and in Fig. 1(d) we plot the group index of the topological edge mode with vanishing group velocity at the cutoff frequency ν = 0.2955(c/a), which corresponds to ν = 177 THz for a = 500 nm. For reference, we use a standard photonic crystal waveguide obtained by leaving out a row of pillars in a triangular lattice of silicon pillars surrounded by air. The diameter of the pillar in this case is set to d = 0.4a with the same lattice unit as the topological waveguide, as shown in Fig. 1(f). The guided mode of this conventional waveguide presents an ideally vanishing group velocity, or diverging group index, when approaching the cutoff frequency of the waveguide at ν = 0.3372(c/a), as shown in Fig. 1(g) and (h). For the sake of clarity, we bring the attention to the fact that the cutoff of the topological and conventional waveguide lie on the X and Γ point, respectively. The ideal spatial field-intensity distributions in both waveguides, i.e., in absence of any perturbation, and at frequencies corresponding to ng = 300 in both cases are plotted in Figs. 1(e) and (i) for reference, evidencing a similar level of light confinement. Using silicon pillars surrounded by air instead of the usual air holes in silicon enables us to flatten the dispersion relation of both topological and conventional guided modes successfully without the need of local perturbations as in Ref. [33] or progressive interfaces as done in Ref. [34]. The parameters to obtain flat bands in the valley-Hall waveguide are detailed in the supplementary material.
The backscattering length, ξ, is the average ballistic propagation distance along the waveguide in absence of any other major loss mechanism [5]. A slow-light waveguide becomes virtually useless when L ξ. Interestingly, ξ is governed by the density of optical states of the waveguide, at least for a weak perturbation of the wavegudie [5,35], as ξ ∝ DoS −1 = (∂ν/∂k). Intuitively, a larger density of optical states induces a larger probability of scattering, thus reducing the value of ξ. As ng ∝ (∂ν/∂k) −1 = DoS, both the group index and the backscatter-ing length are intrinsically linked to each other via the DoS, at least in conventional photonic crystal waveguides [36]. Despite substantial theoretical work on ξ in non-topological electronic [37][38][39] and photonic transport [5,36,40], this parameter has only been explored recently in topological waveguides [41] although ignoring ng. It is important to remark the fact that the backscattering length is a universal parameter in onedimensional transport [42,43] and any mesoscopic-transport observable in these systems depends only on ξ, regardless of the microscopic details of the medium considered. Therefore, two waveguides with a different nanostructured pattern but with the same ξ share the same universal mesoscopic-transport properties, regardless of it being silicon pillars surrounded by air or air inclusions in silicon. We calculate ξ in perturbed topological and conventional photonic crystal waveguides as: where I is the finite-element solution of the electromagnetic field intensity emitted by a dipole at frequency ν, x is the distance from the dipole position along the waveguide and the brackets indicate the statistical ensemble average over different configurations of positional disorder. Our simulation domain has a length L = 200a in the x direction, eleven unit cells on each side of the waveguide in the y direction and it is all surrounded by perfectly-matched layers to mimic an open system. To simulate the effect of fabrication imperfection, we randomize the position of the pillars around their ideal value according to a normal distribution which standard deviation σ is our measure of disorder (more details in the supplementary material). The electromagnetic field intensity excited by a dipole source oscillating at the leftmost edge of the waveguide is shown in Fig. 2. The oscillation frequency is chosen such that both the topological (left) and the conventional (right) waveguide would enable slow light transmission with ng = 300 in the absence of any imperfection. The figure shows calculations corresponding to varying configurations of positional disorder for which the positions of the pillars are randomized with a fixed σ = 0.001a. As shown in Fig. 2(a), the excited Bloch modes in the topological waveguide are just slightly perturbed, which reveals a rather weak backscattering for this large group index value. However, strong backscattering interference prevents light transport in the conventional waveguide, as revealed by the different examples plotted in Fig. 2(b). As ξ is a statistical parameter, this requires an ensemble-average calculation of many (ideally all) different disorder configurations. This ensemble average is not easily accessible in experiments because it is difficult to discriminate backscattering from other loss sources, particularly absorption [20]. In numerical experiments, as the one performed here, ξ is easily obtained by computing the position-dependent field excited by the dipole source in the conditions described above for several structural configurations with the same nominal amount and type of disorder. Figs. 2(c) and (d) show the intensity pattern after ensemble-averaging the electromagnetic field excited by this emitter over twenty different realizations, the envelope of which decays exponentially from the position of the source with a sufficiently well-averaged slope, as plotted along the waveguide axis in Figs. 2(e) and (f) for ng = 300.
The backscattering length ξ calculated for the topological and conventional waveguides for different ng is plotted in Fig. 3. ξ is given for a fixed amount of positional disorder σ = 0.001a, which corresponds to fluctuations of ≈ 0.5 nm for a = 500 nm. This is a realistic measure of the residual imperfection resulting from a state-of-the art fabrication process [21,22]. The figure reveals the fact that, at low disorder levels, the topological waveguide considered here and based on a geometry-engineered topological phase is more robust than standard conventional slow-light waveguides when the full phase is randomized, thus mimicking the effect of imperfection in real systems. As shown in Fig. 3, topological waveguides suffer much less backscattering than conventional ones with a backscattering length comparable to the waveguide length for large group index values. Even at large values of the group index, ng 1000, the interface edge mode mode of the topological waveguide is slightly perturbed when compared to the strong backscattering suffered by the conventional one, as shown in the different configurations plotted in Fig. 2(a). The backscattering mean free path for various amounts of disorder and a fixed ng = 100 is plotted in Fig. 4. Below a critical level of disorder, σc ∼ 0.002a, the topological waveguide outperforms the conventional one. Above this level, the probability of a pseudo-spin flip upon scattering of the Bloch mode increases, leading to decrasing ξ. This behaviour was confirmed for various values of the group index at lower disorder levels. We attribute this to a topological protection which disappears at a sufficiently high degree of disorder, in this work at σc = 0.002a, as shown in Fig. 4. Above σc, the topological waveguide suffers an even stronger backscattering than the conventional one. We understand this behaviour above σc as due to the effective mass of the propagating topological Bloch mode in absence of any perturbation [36]. A larger effective mass (flatter bands) induces stronger backscattering [6,40], which is what we observe in Fig. 4 above the critical level of disorder. In our case, the effective mass of the topological waveguide is almost an order of magnitude larger than the conventional one, which underlines even more the potential of the topological protection observed below σc and for very large values of ng, the main result of this Letter.
In conclusion, the robustness of diverse designs of topological photonic waveguides based on pseudo quantum-spin Hall and valley-Hall effect against structural back-reflection has been claimed in recent experiments [30,[45][46][47], where the sole inclusion of local defects in perfect topological phases was taken into account. However, one of the key elements to quantify the robustness of topological photonic platforms is to randomize the full topological phase to mimic real imperfection due to the fabrication process, which is the real limitation in slow light transport in state-of-the art photonic crystal waveguides. We have shown how to quantify the robustness of a topological edge state against white noise on its structural parameters. Calculating the backscattering length linked to the group index enables us to do that, as these two parameters are related to each other through the density of optical states. We conclude that current proposals of topological photonic phases relying on the breaking of different parity symmetries as the valley-Hall effect are quantitatively, by almost five times, more robust than standard conventional waveguides with small disorder levels, although this protection is lost at higher imperfection amount. We expect this protection to manifest at even larger disorder levels if we take into account the much larger effective masses of the ideal topological edge modes, which plays a dominant role for light localization at frequencies close to the cutoff one. The analysis carried out here is based on a particular system of silicon pillars surrounded by air but the approach is completely generic and can be implemented for any other single-mode waveguide with arbitrary design. Nonetheless, our results point to the importance of the group index ng for a quantitative comparison of the robustness of light transport. The large values of the group index calculated here, ng 1000, affected by a relatively weak backscattering provide a promising platform for highly efficient strong light-matter interaction [48] where photon transmission over hundreds of microns is relatively backscattering-free. Our analysis, in combination with the calculation of the topological invariants applying topological quantum chemistry [49] provides a necessary tool to test the validity of any translation of topological effects from solid-state physics to photonics. Future work to evaluate topological invariants of different topological implementations will provide additional insight on the relationship between the backscattering length evolution with disorder level and the protection granted by topology.

VALLEY-HALL TOPOLOGICAL PHOTONICS.
In time-invariant solid-state systems, topological effects are attributed to mechanisms that break the parity symmetry. In the quantum-spin Hall effect, the spin-orbit coupling plays the role of the external magnetic field to induce the topological protection. In that case, the spin up and down of the electron exhibit chiral and anti-chiral integer quantum Hall effect, respectively. To translate these mechanisms from solid state to photonics, we need to find degrees of freedom that mimic that of the electron spin in order to build up a pseudospin for light [1]. Alternatively, one can implement a valley-Hall effect [2] which exploits the angular rotation of the electron wave function in the K or K valleys of the band structure generating an intrinsic magnetic moment [3] analogous to that produced by the electron spin.
The strategy to implement a valley-Hall topological photonic phase exploits just a single polarization (transverse magnetic in our case) in a lattice where the inversion symmetry is broken in the unit cell. This is fairly more simple than current implementations of the quantum-spin Hall effect for photonics, which require the folding of the dispersion relation. The idea behind this implementation [4] is to create an interface between two different valley-Hall phases with different symmetry-breaking geometries which supports highly-confined topological edge-state guided modes. All the details can be found in the seminal paper written by Ma and Shvets [4].

FINITE ELEMENT CALCULATIONS AND DISORDER.
Most of recent experiments in topological photonics explore slab geometries, i.e., material with a subwavelength thickness. In the particular case of the valley-Hall topological approach [4], the edge states appear below the light line, and are therefore perfectly guided. Disorder is a coupling mechanism between these modes and the radiation continuum, however, in the slow-light region of a band, in plane scattering dominates largely over out-of plane radiation losses [5]. Based on this, we neglect the out-ofplane dimension and implement a two-dimensional geometry, capturing the physics of the problem at a reduced numerical cost. We perform a finite-element calculation using commercially available software to find the solution of the Maxwell equations to this two-dimensional problem. To obtain the variation of the backscattering length in a wide range of group index values, we engineer slow light in a valley-Hall waveguide by creating a lattice formed of silicon pillars surrounded by air. Our design is a valley-Hall waveguide based on a primitive cell of triangular symmetry of size a composed by two dielectric silicon pillars with a dielectric constant εSi = 12 surrounded by air εair = 1, as detailed in Fig. S1(a). The two pillars of the primitive cell have different diameters (d1 = 0.4a and d2 = 0.2a) to break the inversion symmetry. The coordinates of the center of the small and large pillar within the primitive cell are {0, 4 √ 3a/6} and {0, √ 3a/3}, respectively, taking the down-right corner as origin and the x and y coordinates as indicated in Fig. S1(a). For this configuration, the transverse magnetic field presents a frequency band gap between 0.2740 and 0.3727 in normalized units, i.e., in units of a/c. The gap is shaded in red in Fig. S1(c). Finally, we build the valley-Hall waveguide by interfacing two valley-Hall crystals with the same geometry and opposed symmetry, as displayed in Fig. S1(e). In the calculations performed here, the pillars with larger diameter are closer to each other and their centers separated by a distance equal to √ 3/3a. In this configuration, we obtain a single topological edge mode and presents the cutoff at the X symmetry point. We have attempted to engineer slow light in a configuration of air inclusions in silicon without so far success.
For reference, we design a (conventional) non-topological photonic waveguide with a similar configuration based on a primitive cell of triangular symmetry of size a composed by a single dielectric pillar (εSi = 12) of diameter d = 0.4a surrounded by air (εair = 1), as shown in Fig. S1(b). The transverse magnetic field presents a frequency band gap for this lattice between 0.2740 and  Fig. S1(d). The conventional waveguide is built by leaving out a row of pillars creating a line defect in the x direction, as detailed in Fig. S1(f).
We introduce disorder in both waveguides by displacing the position of all the pillars independently a random distance from their position in the perfect lattice, ∆r. This distance is normally distributed with a standard deviation σ = ∆r 2 and ∆r = 0, where brackets indicate the ensemble average over all configurations of random fluctuations.
For the calculation of the unperturbed waveguide, we simulate a numerical domain formed by a single unit cell in the x (propagation) direction and 22 cells in the y direction. As boundaries of our simulation domain we terminate the waveguide with Bloch periodic conditions in the x direction and perfectly-matched layers in the y direction. For the calculation of the disordered waveguides, we expand the numerical domain in the x (propagation) direction to 200 unit cells.
To simulate an open system where light can leak out, we introduce perfectly-matched layers surrounding the whole computational domain. This mimics open boundaries for which electromagnetic waves are transmitted with minimal reflection. Despite the presence of such perfect-matched layers, slight reflections occur at the interfaces, which results in the generation of extremely broad Fabry-Pérot resonances whenever the system size L < ξ. In all our calculations, we implement a mesh with a ratio 150 elements per a 2 area.

CONVERGENCE OF THE RESULTS
To test the convergence of our results, we first test the robustness of the calculated waveguide eigenmode frequencies for an increasing number of elements in the mesh until we observe a saturation in their values. We start with 18 elements per unit area (a 2 ) and we finish with 200 elements per unit area. We observe that the values of the calculated frequencies do not change significantly when we have 150 or more elements per unit area. We fix the mesh with a ratio of elements per surface of the waveguides at 150 elements per a 2 area for all our calculations and, in particular, for all the disordered configurations.
In addition, we use two methods to evaluate the group index and to further test the convergence of our results, in particular the value of the group index. We do perform this security check to ensure that we have the same values of the group index for the eigenmode calculation with periodic boundary conditions and for an extended numerical domain calculation terminated by perfect-matched layers. In one case, we calculate the dispersion relation ν(k) of a single perfect unit cell with periodic boundary conditions. Then, we obtain the ng by fitting the calculated dispersion to a sixth order polynomial p(k). Finally, we derive the Supplementary Figure S2: Group index calculated for a valley-Hall (a) and a conventional (b) waveguide. We use two different methods to obtain these values with two different boundary conditions. First, by fitting the calculated eigenmode solution for a unit cell with periodic boundary conditions six order polynomial. We obtain the group index by simply deriving its value it from the fitting function (black solid line). Alternatively, we simulate a full perfect waveguide with two different lengths and perfect-matched layers at both terminations. Residual reflection for these terminations induce Fabry-Pérot resonances. We obtain the value of the group index directly rom the free-spectral range of these resonances, plotted in empty and solid dots for L = 200a and 400a, respectively. value of the group index ng from this polinomial function as: The values of the group index obtained with this method for the valley-Hall and the conventional waveguide are plotted as a solid black line in Figs. S2(a) and (b), respectively.
The alternative method relies on calculating the eigenmodes of unperturbed realizations with the same calculation configuration that we use to extract the backscattering mean free path. We simulate waveguides with two different lengths L = 200a, 400a and we terminate them in the x (propagation) direction by perfectly-matched layers. As these layers are not 100 % reflectionless and allow some minimal but sufficient reflection, we do observe clear Fabry-Pérot modes corresponding to the modes of an effectively closed cavity with same length. We calculate the value of the group index from the free-spectral range between two Fabry-Pérot modes, ∆νFSR = νn+1 − νn, as: where we approximate the frequency value at which ng is obtained simply by the mean value νn = νn+ν n+1 2 between the resonances. The group index calculated using this method for the two different waveguide lengths, L = 200a and 400a, is plotted in Fig. S2(a) and (b) as empty and solid dots, respectively. The good correspondence between the values of the ng extracted using both methods is a strong convergence test of the meshes employed throughout the whole manuscript. More importantly, it ensures us that the calculated values of the group index used to characterize the backscattering mean free path in the main text are correct and robust.

DETERMINATION OF THE EFFECTIVE PHOTON MASS AND GROUP INDEX
The effective mass of the Bloch mode is defined as m eff = ∂ 2 ν ∂k 2 −1 , where ν and k are the frequency and wave vector. By fitting the band dispersion close to the band gap edge to a second order polynomial, we can extract m eff from the second order fitting coefficient. In the supplementary Figs. S3(a) and (b), we obtain the effective mass corresponding to the propagating Bloch mode of interest for the conventional and valley-Hall waveguides, respectively. In particular, we obtain m Valley = −7.7760π/(ac) and m conventional = 1.1239π/(ac). To obtain a good fit, we limit the data to a small range around the symmetry point. The effective mass of a Bloch mode determines to great extent the backscattering properties of the system [6]. In photonic lattices, the backscattering is enhanced by a high density of optical states, as in flat bands or band edges. In general, these are the spectral regions where backscattering-induced localization appears for lower amounts of disorder. At the band edge of a photonic crystal or at the cutoff frequency of a guided mode it is possible to obtain an effective Schrödinger equation where the random potential derives from the structural disorder of the system and the effective kinetic term, ν(k) ≈ k 2 (∂ 2 ν/∂ 2 k), presents a minimum value for flat bands thus giving rise to bound states and Anderson localization. This is a possible physical picture to understand this process intuitively. Based on their values and close to the cutoff frequency as we are in our calculations, we should expect a much stronger backscattering in the valley-Hall edge state than in the conventional waveguide for the same amount of disorder and group index value.
Supplementary Figure S3: Zoom of the dispersion around the slow light maxima. (a) for the valley-Hall and (b) for the conventional waveguide. The blue solid lines correspond to second order fits used to determine the effective photon mass as described in the text, while red lines correspond to a higher order polynomial fit used to extract ng