Topological supermodes in photonic crystal fiber

Topological states enable robust transport within disorder-rich media through integer invariants inextricably tied to the transmission of light, sound, or electrons. However, the challenge remains to exploit topological protection in a length-scalable platform such as optical fiber. We demonstrate, through both modeling and experiment, optical fiber that hosts topological supermodes across multiple light-guiding cores. We directly measure the photonic winding number invariant characterizing the bulk and observe topological guidance of visible light over meter length scales. Furthermore, the mechanical flexibility of fiber allows us to reversibly reconfigure the topological state. As the fiber is bent, we find that the edge states first lose their localization and then become relocalized because of disorder. We envision fiber as a scalable platform to explore and exploit topological effects in photonic networks.

Topological modes are protected because of the topological invariant that characterizes the system. This integer invariant can only change upon closing a bandgap. By directly associating the invariant with physical characteristics such as light propagation, topology can be used to endow a system with protection against fabricationinduced imperfections. The presence of topological modes relies on a few key properties, including the presence of a bandgap and the symmetries of the system. By selecting the necessary ingredients and structure, topological physics can be used to design mode profiles (12,16), robust pumping (17,18), and unidirectional transport (3,19).
For microwaves, topology can be achieved through macroscale structures in combination with resonance effects (3,5,6). However, at optical frequencies, the challenge remains to create waveguides more than a few wavelengths long in which light is both confined by a micrometer-scale structure and topologically protected. For example, arrays of silicon waveguides are limited by scattering loss in the near infrared (20), whereas for visible light, state-of-the-art planar fabrication techniques, such as using lasers to inscribe waveguides into a glass chip, are challenging to extend beyond the centimeter scale (21). Optical fiber that supports topological modes has previously been proposed (22)(23)(24)(25)(26), but experimental demonstrations have not been seen because of the impossibility of fabricating previous designs with current technology.
Topological fiber could play a crucial role in next-generation quantum networks. Optical fiber is a key enabling technology for long-haul telecommunications networks due to its exceptionally low attenuation (27). In the drive to scale up quantum networks, this low attenuation also plays a key role due to the impossibility of noise-free amplification of quantum states (28,29). In chipscale experiments, light propagating in topologically protected modes over hundreds of micrometers to millimeters has been shown to inherit robustness from the topological medium (11)(12)(13)(14)(15). Therefore, if the guidance of topological modes could be demonstrated in fiber at length scales of meters to kilometers, this would offer topological robustness within distributed classical and quantum networks. In addition, topological fiber design has the potential to protect the generation of entangled states of light from the detrimental effects of fabrication-induced variations (30).
Unlike telecommunications fiber, photonic crystal fiber (PCF) controls the propagation of light using a cross section with wavelength-scale microstructure (31,32). Tailoring the design of PCF allows unparalleled customization of its dispersion and mode spectrum (33). For example, many solid glass cores can be embedded in a periodic cladding of air holes to form a single flexible fiber with a typical diameter of 100 to 250 μm. In this fiber, the individual modes of each core can overlap and collectively form so-called supermodes. This type of multicore structure forms the platform for our topological PCF (TopoPCF) that translates the strategies for topological band engineering into fiber. Fabricating our TopoPCF allows us to observe propagation of edge states with visible light and to use bending for topological mode control. We show that bending the fiber leads to an effective disorder inside our topological state, which we can switch on and off through mechanical reconfiguration. This unique control of topology via disorder exemplifies fiber-based phenomena that are inaccessible in planar waveguide architectures.
We create a topological lattice that propagates light over meter scales governed by its micrometer-scale cross section by engineering the coupling between cores (Fig. 1A and the "Fiber fabrication" section). We focus on a canonical model exhibiting a topological invariant, the Su-Schrieffer-Heeger (SSH) chain (34), defined by identically shaped cores and alternating intercore coupling strengths. To achieve this, the design alternates between small and large air holes between neighboring cores. The two ends of the chain include an additional small air hole to create a chain of cores uniform in shape. We curl the chain up inside the fiber to form a spiral (Fig. 1, A and B) and observe that both ends of the spiral host topologically protected edge states.
Light is guided by fiber in normal modes. Consider a single optical polarization of the transverse electric field E at wavelength λ, which we denote ψ(x, y) (=E x or E y ). The eigenvalue equation for ψ(x, y) and the modal propagation constant β (that is, the longitudinal z component of the wave vector) takes a form analogous to the time-independent Schrödinger equation. In this analogy, the square of the local refractive index n(x, y) plays the role of the confining potential for light. Explicitly, the eigensystem reads where r 2 ? ; ∂ 2 x þ ∂ 2 y is the transverse Laplacian and k 0 = 2π/λ is the free-space wave number (see derivation in the "Theory of the photonic SSH chain" section). In our system, high symmetry within each core leads to minimal polarization coupling so that the mode profiles are degenerate under rotations of the electromagnetic fields. A solution of Eq. 1 is an electric field profile ψ(x, y) that propagates unchanged in time t along the longitudinal z direction as ψ(x, y)e i(ωt − βz) , with angular frequency ω = 2πc/λ in terms of the speed of light in vacuum, c.
The similarity between Eq. 1 and the Schrödinger equation inspires the design of fiber using analogies with electronic band structures. For example, the tight-binding model predicts that for weakly coupled cores, supermodes obey the matrix equation Cu = Δβu, where Δβ is the change in propagation constant for a supermode relative to the single-core propagation constant and u = (a 1 , b 1 , a 2 , b 2 , …) is a complex vector describing both amplitudes and phases of the single-core modes constituting the supermode. For the SSH chain (Fig. 1B), the coupling matrix C contains just two (nearest-neighbor) elements: We denote C 1 to be the coupling strength between cores A m and B m within the mth unit cell and C 2 to be the coupling between B m and A m+1 across unit cells.
The topological invariant ν is encoded in the phase differences between optical waves propagating in neighboring cores. When the chain is infinitely long, the supermodes u can be decomposed into an eigenbasis of supermodes defined by the relations a m+1 = e iθ a m and a m = ±e iϕ b m . For these normal modes, each core contains the same light intensity, with the phase difference ϕ between the two sublattices within a single unit cell and the phase difference θ between cores on the same sublattice and in neighboring cells, as illustrated in Fig. 1C. In other words, θ is the reciprocal-space variable along the chain. The phase differences ϕ and θ are related via the coupling matrix. In the eigenbasis, C(θ) is a 2 × 2 antidiagonal matrix with entries ζ(θ) = C 1 + C 2 e iθ (and its complex conjugate; see the "Theory of the photonic SSH chain" section). The phase ϕ is fixed as the complex argument of the function ζ(θ), plotted as circles in Fig. 1D. The map ϕ = arg ζ(θ) relates two phases and is, therefore, characterized by an integer invariant called the winding number ν. If the circle ζ(θ) circumscribes the origin, then the winding number ν of the SSH chain is equal to one, which occurs when the intracell coupling is weak, C 1 < C 2 . For strong intracell coupling C 1 > C 2 , this topological invariant is zero. Mathematically, we compute Physically, as we change θ by a full period, we note that for ν = 1, the phase difference ϕ between cores A and B undergoes a full period (2π) shift (Fig. 1D, top boxes). For ν = 0, ϕ stays close to zero (Fig. 1D, bottom boxes).
Topological band theory predicts that an integer invariant characterizing an unbounded system has profound consequences once boundaries are introduced. Because of the alternating coupling coefficients C 1 and C 2 in the SSH chain, a gap opens up (centered The topological invariant ν is encoded in the phase differences (ϕ within a unit cell and θ for neighboring cells) between optical waves propagating in neighboring cores (see main text). (D) The function ζ(θ) has argument ϕ and draws a circle in the complex plane. The circle winds around the origin for C 1 < C 2 (topological, red) but not for C 1 > C 2 (trivial, blue). Boxes: Phase difference ϕ between A cores (solid) and at the single-core propagation constant β), which has a gap size of 2|C 1 − C 2 | (see the "COMSOL simulations" section). No bulk modes are allowed to have propagation constants inside this bandgap. Nevertheless, the bandgap can host modes localized at a boundary, and the topological invariant encodes information about these boundary modes. A heuristic argument explains this connection. An invariant can only change discontinuously and only when the bandgap closes, such as at a boundary between a topological chain and a trivial cladding. This transition necessitates a closing and reopening of the gap, guaranteeing a mode to exist at the boundary. Because this argument does not rely on the boundary's geometry, the band topology endows this edge mode with protection against disorder.

RESULTS AND DISCUSSION
We fabricated a PCF from silica glass by the stack-and-draw process. Our fiber contains 12 cores coupled to form an SSH chain, as illustrated in Fig. 1A. Using simulations, we estimate the coupling strengths within the fiber to be C 1 = 166m −1 and C 2 = 352m −1 for light with a wavelength of 700 nm. We verify that our fiber supports chain-edge supermodes through a combination of simulations and experiments shown in Fig. 2. In Fig. 2 (A and B), we numerically solve Maxwell's full vector equations, finding the propagation constants and field profiles of the supermodes supported by our fiber. These supermodes can be separated into two classes: bulk states (labeled 2 to 11) and exponentially localized edge modes (1 and 12). The bulk states predominately contain contributions from cores that make up the body of the chain and have propagation constants that lie in the bulk bands, above and below a topological bandgap (see the "COMSOL simulations" section). Two topological edge modes are found within this bandgap. These edge modes exist at Δβ ≈ 0, are highly localized to the end of the chain, and are robust to disorder smaller than the bandgap. In experiments (Fig. 2C), light is injected separately into each of the 12 cores by butt-coupling from a single-mode PCF. For each input condition, we show the intensity distribution after light has propagated through the fiber. When light couples into cores 2 to 11, it excites a combination of bulk modes. When exciting bulk modes, light couples across all cores in the body of the chain and, after propagation, exits as a distribution of intensities across the bulk cores [blue in Fig. 2 (C and D)]. By contrast, light coupled into the two edge cores predominantly excites the topological modes and remains confined to the edge [red in Fig. 2 (C and D)]. Notably, the topological edge states exhibit an asymmetry between the A and B sublattices, showing an alternating intensity profile characteristic of the SSH chain. This intensity profile arises because of the interaction between the topological edge modes and the constraints enforced by sublattice symmetry.
To confirm the topological origin of the edge states in Fig. 2, we directly measure the winding number ν (21,35,36). We augment the method from (21), with the notable difference that instead of changing the propagation length, we scan over the optical wavelength λ using the setup in Fig. 3A (see movies S1 and S2 for the wavelength-dependent output and the "Characterization" section). With light injected into a core in the bulk, we scan the optical wavelength and observe a change in the intensity distribution at the output of the fiber. This change sees light move from cores belonging to the A sublattice to the B sublattice. In the "Weighted intensity difference and the winding number" section, we connect this change in sublattice intensity distribution with the winding number that characterizes the system, giving us access to an experimental observable defining the distinct topological invariants. At the output, the topological invariant ν can be computed as twice the weighted intensity difference I d (λ), also called the mean chiral displacement (35), between the A and B sublattices averaged over all wavelengths λ where I A = |a m | 2 and I B = |b m | 2 are the intensities in each sublattice and the sum runs over all unit cells m (see derivation in the "Weighted intensity difference and the winding number" section).
We choose a wavelength range sufficiently broad to capture the characteristic oscillations of I d (λ). Furthermore, we average I d (λ) over a wavelength subrange such that the quantity ⟨I d (λ)⟩ λ remains invariant under any renumbering of the unit cells. To do so, we pick the wavelength interval over which the measured unweighted intensity difference ∑ m (I A − I B ) vanishes, as explained in the "Weighted intensity difference and the winding number" section. We then use the SD of this unweighted intensity difference as an unambiguous choice for the error. We plot 2I d (λ) in Fig. 3B and measure 2〈I d (λ)〉 λ to be 1.01 ± 0.33, confirming bulk topology and bulk-boundary correspondence. The simulations agree with our experimental conclusions and additionally demonstrate that Bulk-boundary correspondence and the localized edge modes are a consequence of the topological invariant characterizing the system. This invariant must always exist, provided the chain displays sublattice symmetry and experiences sufficiently small disorder (that is, smaller than the photonic bandgap). In contrast to a topologically trivial system, the edge modes in our fiber always exist, provided that the invariant remains well defined. Thus, we can think of the edge modes as inheriting the robustness of the topological invariant. This robustness reveals itself when disorder is introduced to the coupling strengths. If the couplings in the chain contain disorder, topologically trivial modes can be perturbed and even destroyed, whereas topological edge states remain immune. For cases that respect sublattice symmetry, the topological mode remains protected as long as the ratio C 1 /C 2 is less than unity. By contrast, if the disorder violates sublattice symmetry, then the protection of the edge modes is weaker.
When on-site disorder is introduced into our chain, the sublattice symmetry is broken. Under these conditions, the topological edge mode loses its guaranteed protection and can now be destroyed while the ratio C 1 /C 2 is less than unity. However, for the mode to be destroyed, the on-site disorder must be sufficiently large when compared to the size of the bandgap (37). In the "Preserving sublattice symmetry" section, we explore the effects of onsite disorder that exist in our fiber. We find that for realistic structural imperfections, the strength of disorder is less than 3.2% of the photonic bandgap. For these small values of the disorder, the edge modes remain identical to the ideal case that respects sublattice symmetry.
Our fiber platform gives rise to topological phenomena inaccessible in short, rigid waveguides and solid-state devices. One such phenomenon is the reconfigurable control of topological states through geometric bending, only accessible because of the unique physical flexibility of fiber (Fig. 4). In the presence of a bend with radius R, the propagation constant of each core changes because of this geometric perturbation by Δβ ≈ 2πΔx/(λR) for cores separated by distance Δx [see the "Numerical propagation" section and (38,39)]. This geometric control enables us to globally reconfigure the SSH chain with a single macroscopic bend. Such reconfigurability is in contrast to metamaterial systems where individual unit cells can be adapted, but changing the entire lattice presents a challenge.
The changes in propagation constant due to fiber bending map directly onto on-site disorder in the system's coupling matrix. This mapping enables us to analytically predict the response of the topological edge mode to bending. Because on-site disorder does not preserve the chiral symmetry of the SSH chain, the topological invariant is no longer guaranteed. This loss of symmetry destroys the localizing topological protection and allows light in the edge mode to spread across the bulk modes. We observe this in the insets in Fig. 4A: As the curvature is increased, the topological edge mode becomes delocalized. Upon further bending, this trend begins to reverse and the on-site disorder grows sufficiently large to induce Anderson localization (40,41), which reconfines light within the initially excited core. To quantify this localization, we calculate the inverse participation ratio (IPR) IPR ¼ ð P k I 2 k Þ=ð P k I k Þ 2 , where the sums run over all cores k = 1 to 12. We plot the IPR as a function of fiber curvature in Fig. 4. We predict the breakdown of topological protection at the bend radius R * , at which the ratio of disorder Δβ to the average coupling strength C is near one [Δβ ≈ C = (C 1 + C 2 )/2]. As found in the "Numerical propagation" section, 1/R * ≈ 1.28λC/(2πΔx) ≈ 3m −1 , in agreement with Fig. 4B. This bending radius characterizes a simple reversible mechanism for globally breaking and restoring topological protection.
To conclude, we have used the stack-and-draw technique to make hundreds of meters of topological fiber for visible light. We observe the topology and control the corresponding edge states The relative intensity of the electric field for each supermode (doubly degenerate due to polarization). The two edge states (in red) correspond to supermodes that we label 1 and 12. (B) Comparison of the intensity profiles between bulk mode 6 (blue) and edge mode 12 (red). (C and D) Light propagation experiments: We inject light (wavelength, 636 nm) into a core (input core number) of a TopoPCF with a length of (17.0 ± 0.1) cm and image the output intensity profiles. For input cores 2 to 11, light excites bulk modes and spreads through the chain. For input cores 1 and 12, light remains localized to the topological edge states. (D) Quantitative comparison between the edge mode excited via core 12 and a bulk mode excited via core 4. The topological edge state demonstrates an asymmetry between A and B sublattices characteristic of the SSH chain.
through fiber bending, a mechanism inaccessible to on-chip waveguides. We envision topological fiber as a flexible platform that allows topological effects to be applied in a scalable way to classical and quantum photonic networks.

Theory of the photonic SSH chain
Light propagating in a waveguide is described by Maxwell's equations in matter with no free charges or currents, constant permeability μ ≈ μ 0 , and spatially varying permittivity ϵ/ϵ 0 = n 2 (x, y) Using the curl-curl identity and the expansion of the first equation in Eq. 4, we obtain the vector wave equation The final term may be ignored in our system, as the variation in the refractive index is negligible everywhere the electric field is nonevanescent (i.e., in the glass). Translation invariance in the z direction along the fiber allows us to write solutions to Eq. 5 as superpositions of modes of the form E ¼ ½cðx; yÞê�e iðvtÀ βzÞ , with angular frequency ω (which is the same in glass as in air), propagation constant β, and a unit vectorê pointing perpendicular to the fiber.
Substituting these modes into Eq. 5 gives the scalar wave equation where r 2 ? ; ∂ 2 x þ ∂ 2 y , which is Eq. 1 in the Introduction with k 0 = ω/c. Let n g (n a ) be the refractive index of the glass (air) and δn = n g − n a their difference. In our multicore fiber, the air hole geometry can be discretized with hexagonal platelets to obtain an effective refractive index n(x, y), which evaluates to n g on the cores, n g − f 1 δn around the small air hole flanking the core, and n g − f 2 δn around the larger air holes everywhere else. Here, f k ¼ p 2 ffi ffi 3 p d 2 k =Λ 2 is the air hole area fraction in terms of hole diameter d k and fiber pitch Λ (see Fig. 5).
We expand ψ as a supermode ψ = ∑ j u j ψ j in terms of single-core modes ψ j . At the first approximation, we can solve for ψ j on the jth core in the absence of other cores. That is, we use Eq. 6 with a coarse-grained refractive index n j (x, y) where every ith core, with i ≠ j, is replaced by a large air hole, that is where Δn i evaluates to f 2 δn on the ith core and to zero everywhere else. Note that, because the small air hole alternates between being on the left and right sides of the core, the single-core mode on an A (B) site will be left(right)-heavy, as depicted in Fig. 5.
Writing the shift in propagation constant Δβ for a supermode compared to the single-core β, we subtract the wave equation for the individual cores from the wave equation for the supermode to , where I A, B is the normalized sublattice intensity within one unit cell and the sum is over all unit cells, m = 1…6. We vary the wavelength over a 104-nm range and propagate through 17 cm of fiber. We calculate the mean of 2I d to be 1.01 ± 0.33 over the highlighted range to confirm the system's nontrivial winding number (see main text). Error bars are set by wavelength filter resolution, and the gray confidence interval is the SD of the unweighted intensity difference ∑ m (I A − I B ). Inset: Simulation data for I d comparing fiber geometries with winding numbers 0 (blue) and 1 (red). For comparison, the topological simulation average is also plotted as a red dashed line in (B). obtain X j u j ½2βΔβ À k 2 0 ðn 2 À n 2 j Þ�c j ¼ 0 Multiplying Eq. 8 on the left by a specific core-mode 〈ψ k | while approximating n + n j ≈ 2n g and β ≈ k 0 n g , we have X where we use bra-ket notation for convenience. Since the singlecore modes are strongly localized, the left-hand side is dominated by 〈ψ k |ψ k 〉 ≈ 1, whereas sum on the right-hand side is dominated by 〈ψ j±1 |Δn j±1 |ψ j 〉 ≪ 1. Hence, for the dominant terms, we must have Δβ ≪ k 0 and which is our starting point for the tight-binding approach in the main text. These real coupling constants C kj = k 0 〈ψ k |Δn k |ψ j 〉 are symmetric because every neighboring pair of single-core modes are reflections of each other. Furthermore, this alternating left/ right asymmetry of the single-core modes allows us to define the notation C 1 = C 2m−1,2m and C 2 = C 2m,2m+1 with C 1 < C 2 .
We rewrite Eq. 10 in bra-ket notation as C|u〉 = Δβ|u〉 and the ket |u〉 = ∑ m a m |A m 〉 + b m |B m 〉 in terms of a m = u 2m−1 and b m = u 2m . Then, our coupling matrix can be expressed as Note that the second term introduces coupling between different unit cells (different m). To simplify the analysis, we use a discrete Fourier transform (DFT) to switch to a different basis enumerated by discrete angles θ k = 2πk/M for k = 1, …, M j Aðu k Þ; In this basis, C couples only between states with the same θ k where ζ(θ k ) = C 1 + C 2 e iθ k . Note that Eq. 11 contains coupling to the extraneous core A M+1 , which the DFT equates with core A 1 .
To remove these periodic boundary conditions, the standard method is to state that the boundary is sufficiently far away by sending M → ∞.
In the continuum limit, θ lives on the interval [0,2π]. We use Eq. 13 to solve the eigenproblem C|u〉 = Δβ|u〉 and obtain two solutions for every θ, given by eigenvectors where ϕ(θ) = arg ζ(θ), and eigenvalues The two bands of eigenvalues are separated by a gap at θ = π of size 2|C 1 − C 2 | around the single-core propagation constant β. Inside this gap, no supermodes can propagate through the bulk of the chain. In terms of the components of u = (a 1 , b 1 , a 2 , b 2 , …), the

eigenvectors read
which forms the starting point of our definition of topological invariant ν in the Introduction, as defined in Eq. 2 and illustrated in Fig. 1 (C and D).

COMSOL simulations
We numerically solve the full Maxwell's vector equations defined in the "Theory of the photonic SSH chain" section to find modes of the form Eðx; y; zÞ ¼ e Eðx; yÞe iðvtÀ bzÞ , including the propagation constants β and transverse field profiles e Eðx; yÞ of the normal modes supported by the fiber cross section. We confirm that the coupling to the polarization of light is sufficiently weak for modes to take the form E ¼ ½cðx; yÞê�e iðvtÀ bzÞ for arbitrary in-plane polarizationê, where we take the time variable t to be fixed. We solve this equation for a finite cross section selected to match the parameters of the fabricated fiber. Using COMSOL Multiphysics, we run finite-element analysis simulations with varying meshes, materials, and boundary conditions. The COMSOL design process is illustrated in fig. S1. In all our simulations, a bipartite mesh was used, with the core sections featuring a denser mesh than the surrounding cladding. The boundary conditions remain fixed for all simulations: A perfect electric conductor boundary was defined at the perimeter of the cladding. We model the refractive index of silica by implementing in COMSOL a wavelength-dependent Sellmeier equation with the refractive index of the air holes taking a constant value of 1.
After specifying materials, mesh sizes, and boundaries, the spectrum of normal modes is calculated. For our TopoPCF, these modes are plotted in Fig. 6. The electric field data from these modes are plotted as |E(x, y)| 2 /max(|E| 2 ) to show relative field intensities and compare to experimental observations. In Fig. 2A, we directly plot the simulated modes, and in Figs. 2B, 3B, and 4B, we integrate the intensity across each light-guiding core. We use this integration to normalize the supermode vector u such that the intensities sum to one.
To simulate fiber bending, the spatial profile n(x, y) of the refractive index of the material was modified, which allows normal modes to be calculated in a bent fiber using our two-dimensional (2D) model of the cross section. Following the work in (38,39), the refractive index was modified to be n bend = n silica (1 + q/R eff ), where R eff = 1.28R bend and q is the perpendicular distance to the bend axis, to accommodate the change that light would experience in a bent fiber. As highlighted in (38), an effective bending radius R eff was used to account for both geometric and stress-optic effects. A selection of the intensity plots found through bending simulations can be found in Fig. 7.
The edges of our SSH chain define the system's boundary conditions by determining whether the end of the chain cuts a topological unit cell. For systems with an even number of cores, symmetry necessitates that both ends feature the same boundary conditions. However, if the chain features an odd number of cores, then one side will cut a topological unit cell and subsequently lose its edge mode. Once the system features an odd number of cores, there is no way to remove the topological edge (while respecting sublattice symmetry). The topological edge state can simply be moved from one end to the other by adjusting the coupling strengths and, thus, unit cell definitions. In Fig. 8, we show through numerical simulations the difference under boundary conditions and the presence of a single edge mode. Figure 8A shows the propagation constants of an 11-core chain that both starts and ends on a core belonging to the B sublattice. The key feature in Fig. 8A is the single mode lying in the bandgap. We plot the intensity distribution of this mode in Fig. 8B to show that only the topological edge mode at the center of the fiber now remains. Figure 8C shows that core 11 (the outer end core of the chain) now contributes to bulk modes of the system and has been shifted out of the bandgap. While this change destroys a single topological edge mode, sublattice symmetry changes are needed to destroy the topology everywhere in the chain.

Preserving sublattice symmetry
To confirm that our results are not due to localized disorder, we must ensure that all cores are sufficiently similar. Local structural uniformity ensures that all of the light-guiding cores have the same propagation constant. Uniform propagation constants are required to make a direct mapping to the tight-binding model. In our system, all cores have the same immediate environment (up to rotational symmetry). However, because of imperfections induced in fabrication and the perturbations that each core may experience from other cores in the fiber, this symmetry cannot be guaranteed in practice. That is, chiral and translational symmetries are only   Fig. 2A. The histogram highlights the gap between bulk bands (blue) in the middle of which we find the two topological edge modes (red). As predicted by the analytic theory, these modes exist at Δβ ≈ 0. approximate in our fiber. The lack of hard symmetry constraints means that, theoretically, a topological phase transition can occur without the full closing of the bandgap. However, the topological invariant that characterizes the system remains well defined as long as the bandgap is much larger than the variation in on-site potentials (37).
We confirm by numerical simulation that on-site disorder is two orders of magnitude smaller than the bandgap size when considering experimental parameters. We consider three changes in geometry that influence the propagation constant of a single core, shown in Fig. 9. The first two modifications could arise as a result of fabrication imperfections. The two different sizes of air holes immediately adjacent to a core could be over-or underinflated. This would directly change the propagation constant and, hence, the effective on-site potential. We model in COMSOL the effects of changes in air hole diameter of up to 5%; however, we typically observe less than 1% change in the geometric parameters of fabricated fiber, as demonstrated in (30). We find that with a 1% change This edge mode stays fixed to the core at the center of the fiber because, at that edge of the chain, the topological unit cell is not cut. That is, the boundary conditions at the center are the same as in the 12-core fiber that we fabricated, but the boundary conditions at the other end of the chain are different. (C) Intensity plot for a bulk mode, which now contains contributions from core 11 at the outer end of the chain. As expected, because the unit cell is cut at this end, this edge no longer hosts a topological state. to air hole radius, the core's propagation constant changes by only 3.2% of the bandgap size. In Fig. 9B, we consider how the propagation constant of a core changes when introducing a new neighboring core. When a second core is coupled to the first, the band structure splits into two supermodes. In the tight-binding approximation, the cores do not affect each other's on-site potentials and the propagation constants of the supermodes are equally split above and below the single core's propagation constant (see Fig. 9B). To estimate the shift in the single-core propagation constant from our COMSOL simulations, we then take the average of the two supermodes. This shift between the averaged supermode propagation constant and the single-core case allows us to estimate the change in on-site potential as a result of the second core. We calculate this shift of the on-site potential to be 0.8% of the system's bandgap size.
To verify the assumption that the invariant remains protected, we now insert these small changes into the tight-binding model and investigate the resultant mode structure. In Fig. 10A, we define the distribution from which our random disorder is drawn. This distribution features an SD of 1.6% change in propagation constant when compared to the bandgap. Figure 10B shows the average propagation constants for the modes of our disordered system, with the black lines showing the SD. We take 10,000 iterations of random on-site disorder from our distribution in Fig. 10A and apply these to our tight-binding model. Figure 10C shows the average edge mode intensity distribution found over the same 10,000 iterations, with the black lines again showing the SD. Figure 10 (B and C) shows no significant deviations from the ideal case without disorder and still demonstrates a complete topological bandgap and a sublattice-polarized edge mode. In Fig. 10 (B and C), we use an 11-core SSH chain in place of our 12-core experimental system. This calculation was chosen because an 11-core chain only has one topological edge mode, at Δβ = 0.

Numerical propagation
Investigating light propagation begins with defining an input vector u(z = 0), describing the complex amplitude of each individual core mode in the initial excitation. The overlap integral between this input vector and each of the 12 supermodes (u 1 , u 2 , …, u 12 ) gives the contribution of each supermode present in the initial excitation. By definition, these supermodes propagate along the fiber with only a change in phase (at any fixed time t): u 1 (z) = u 1 (0)e −iβ 1z . By multiplying each supermode by its contribution and phase change, we find the final description of the supermode at a distance z from the input location. After finding the final vectors describing each supermode at the desired distance, we sum these contributions to find the total field distribution in the supermode basis. To return to the core basis, we then perform a second set of overlap integrals for this final supermode vector with each individual core vector, returning a vector containing the contribution of each core to the final state. Multiplying this final state by its complex conjugate returns a vector I(z) = (|a 1 | 2 , |b 1 | 2 , |a 2 | 2 , |b 2 | 2 , …) of relative intensities, describing how much light is in each core in the chain at distance z.
In addition to plotting the propagation of light, the benefit of this numerical propagation is that it allows simulation of the wavelength dependence of the sublattice intensity difference I d (λ) used in Eq. 3. By finding the supermodes for wavelengths in the range shown in Fig. 3B, numerical propagation can be used to model the output of the fiber after propagating through a given distance z. This semi-analytical method allowed us to make use of the detailed simulation results without computationally expensive 3D models, while giving more accurate predictions than a purely analytical approach based on the tight-binding model.
Following the "COMSOL simulations" section, the change in propagation constant is We modify the coupling matrix C using diagonal elements corresponding to Δβ for each core. In general, each core's Δβ will have a different value, which, in addition, depends on the bending axis, therefore playing the role of (diagonal) on-site disorder. To separate the general effect of bending from the specific geometry of our design, we now bend this fiber over a range of different bending directions and average the results.
We consider 40,000 different bending directions for each point on the x axis and find the eigenvectors for every system. The IPR values for these eigenvectors are plotted as a 2D histogram in fig.  S2. Here, we see that at low curvature, there are two populations of modes with different IPRs: the localized topological edge states with IPR ≈0.6 and the bulk states with IPR ≈0.15. As we are interested in the breakdown of topologically protected states, we now discard the bulk states and focus on the effects displayed by the edge states. We select the two eigenmodes corresponding to topological edge modes (in every bend realization) and average their IPRs over all 40,000 points at each step, leaving us with the average effect of a bend on the IPR of a topological edge mode. Averaging the variances of each set of 40,000 bends and then taking the square root allow us to fit an average SD for the entire distribution, as shown in the curve labeled "tight-binding" in Fig. 3B.
To estimate the point at which the topological edge mode delocalizes, we consider the case when the on-site disorder becomes equal to the average coupling strength (Δβ ≈ C = (C 1 + C 2 )/2). At wavelength λ = 500 nm, the critical curvature corresponding to this disorder is found to be: 1=R � bend � 1:28lC=ð2pΔxÞ � 3m À 1 , which is the value quoted in the "Results and discussion" section.

Fiber fabrication
The TopoPCF shown in Fig. 1A was fabricated using the stack-anddraw process. The first step was to create a macroscopic preform consisting of silica rods and tubes, stacked by hand in the geometry corresponding to the fiber design. As the design required two air hole sizes, we used silica tubes with two different ratios of outerto-inner diameter. The precursor to the larger air holes in the cladding was a tube with outer/inner diameter of 25 mm/10 mm, whereas the smaller air holes, found in between every other pair of cores, were created from a 10-mm/3-mm tube. The cores were formed from a solid silica rod. Each element was drawn to a diameter of 1.05 mm and stacked into a regular triangular array with the desired layout.
The completed stack was inserted into a 25-mm/17-mm silica tube and drawn into 4-mm-diameter canes, containing the hole pattern required for the fiber. Fiber preforms were created by jacketing canes in 6-mm/4-mm tubes with brass fittings to enable pressure to be applied to the holes in the cane to prevent collapse during the fiber draw.
The fiber was drawn to a diameter of 162 μm with a high-index ultraviolet-cured polymer coating applied during the draw to prevent damage. The PCF used in this work had a pitch (the center-to-center distance for two neighboring air holes) of approximately 4.4 μm, with a large air hole diameter of approximately 1.6 μm and small air hole diameter of approximately 1.0 μm. Figure 3A highlights the key features of our experimental setup. We used a 1064-nm subnanosecond microchip laser to generate a broadband supercontinuum spanning the visible and near infrared in a length of single-core PCF. This provided white light in a single spatial mode from which a tunable wavelength range was selected by the monochromator shown in the green box in Fig. 3A. The monochromator was formed of a folded 4-f line containing a reflective diffraction grating and a mechanically adjustable slit to pass the desired wavelength range with a bandwidth of approximately 2 nm in the spectral region of interest. The pass-band was collected by a length of endlessly single-mode (ESM) PCF for delivery to the TopoPCF. The pitch of the ESM PCF was approximately 2 μm, which is less than half that of the TopoPCF, enabling the ESM PCF to selectively address individual cores in the TopoPCF by butt-coupling. After propagation, a magnified near-field image of the output of the TopoPCF was formed on the camera sensor by an aspheric lens. The ESM PCF was positioned with a three-axis flexure translation stage to couple light into each core of the TopoPCF sequentially to obtain output intensity distributions for all bulk and edge excitations shown in Fig. 2C. The total intensity  Fig. 9C, we include a normal distribution of on-site disorder in our tight-binding model. This graph shows the distribution of on-site disorder for the 10,000 instances for which data are shown in this figure. (B) Average propagation constant in the presence of disorder compared to the system with no disorder. The black lines show the SD (smaller than symbol size). We conclude that realistic experimental disorder does not change any of the conclusions of our work. (C) Intensity profile for a single edge mode, averaged over the disorder and compared to the case without disorder. The black lines show SD (which is so small compared to the average intensity that it is not visible).

Characterization
at the output of each core was calculated from the images using Python.
To make the direct observation of the topological invariant shown in Fig. 3B, we first established a wavelength range over which the bandwidth of the input light remained approximately constant. We tuned the pass-band of the monochromator over a 104-nm range between 632.8 and 727.8 nm, and using an optical spectrum analyzer, we observed that the bandwidth of the transmitted light did not change significantly. This process also enabled the slit position of the monochromator to be calibrated with transmission wavelength. After calibration, the light was coupled from the ESM PCF into core 2 in the TopoPCF. Small wavelength steps (average step size, 2.7 nm) were made by varying the slit position in the monochromator without changing any other components of the system. At each wavelength step, the output near-field image from the TopoPCF was recorded. From these images, the total intensity difference in each unit cell was found.
The bend measurements required additional components in the optical setup. An acrylic board was laser-etched to create circular grooves with 10 different radii, corresponding to the curvatures seen in Fig. 4B. The board was then fixed in place between the input fiber and the camera to hold the TopoPCF at specific bends. After coupling light into core 12, the input was held constant while the fiber was moved between grooves in the board. For each curvature, the output end of the TopoPCF was imaged onto the camera. This process was repeated three times for two different rotational orientations of fiber to minimize any specific geometric effects from the axis of the bend. The IPR was calculated from the total intensity I k in each core, where k = 1…12, using IPR ¼ ð P k I 2 k Þ=ð P k I k Þ 2 . The results are plotted in Fig. 4B.

Weighted intensity difference and the winding number
Here, we derive Eq. 3, which allows us to compute the winding number from propagation measurements. We then justify our analysis for this experimental winding number computation, including the error estimate. We consider how the intensity distribution changes as light propagates through the system, from being input into a single A or B core in a unit cell labeled by the integer m ⋆ . We define two operators for the weighted intensity difference I d and the unweighted intensity difference e I d as where I Am = |A m 〉〈A m | is the intensity of light in the A site (and I Bm = |B m 〉〈B m | is the intensity in the B site) of the mth unit cell and m runs from 1 to 6 in our fiber. The evolution of these quantities along the fiber can be described in terms of the coupling matrix C multiplied by the propagation distance z to yield the expectation value I d (z) We take a Fourier transform along the chain as if the chain were infinitely long, as described in the "Theory of the photonic SSH chain" section. Identifying 〈A(θ)| = (1,0) and 〈B(θ)| = (0,1) in Fourier space, we replace m → i∂ θ and I Am − I Bm → σ z . We write the vector ζ = (Reζ, Imζ, 0) using complex notation, such that C = ζ·σ. Here, σ = (σ x , σ y , σ z ) is the vector of Pauli matrices We note that our definition of I d (z) depends on the labeling of the input cell m ⋆ . Supposing, for notational simplicity, that the light starts out in the B site of this unit cell, i.e. |u(0)〉 = |B m ⋆〉, we have ð0 1Þðe þim w u e þiðz�sÞz i∂ u σ z e À iðz�sÞz e À im w u Þ 0 1 Here, I w d ¼ P m ðm À m w ÞðI Am À I Bm Þ is the centralized weighted intensity difference. Note that I d ¼ I w d if the input cell is labeled m ⋆ = 0.
We first show that the winding number ν can be computed using the wavelength average of I w d ðzÞ. To start, the integrand of I w d ðzÞ can be simplified by first expanding the exponent of the Pauli vector e +iðz�sÞz ¼ cos ðj z j zÞ1 + i sin ðj z j zÞðn � sÞ ð24Þ wheren ¼ z= k z k¼ ðcos f; sin f; 0Þ in terms of the phase ϕ = arg ζ(θ). Using the product of Pauli vectors ðw � sÞðv � sÞ ¼ ðw � vÞ1 þ iðw � vÞ � s and the commutation relation ðn � sÞs z ¼ À s z ðn � sÞ because ζ lies in the xy plane, the operator in the I w d ðzÞ integrand reduces to e iðz�sÞz i∂ u σ z e À iðz�sÞz ¼ iσ z ½c1 À isðn � sÞ� where, for compactness, we abbreviated c = cos(|ζ|z) and s = sin(|ζ|z). When integrated against θ, the total derivative of a single-valued and periodic function vanishes by the fundamental theorem of calculus. For the remaining term, we note that ðn � ∂ un Þ � σ ¼ ð∂ u fÞσ z , leaving us with In our experiment, we probe λ in the range of 600 to 750 nm. From simulation data, we know that over this range, the coupling constants C 1 and C 2 are well approximated by a linear dependence on λ, and the coupling ratio C 1 /C 2 varies on the order of 1% per 100 nm. As a direct consequence, we can take, for every θ, the directionn as constant and the magnitude |ζ| as linear with respect to λ, with a slope between |∂ λ C 1 − ∂ λ C 2 | = 31 m −1 and |∂ λ C 1 + ∂ λ C 2 | = 92 m −1 per 100 nm. Hence, if we take measurements at a length z > π/31 m −1 ≈ 10 cm, then the quantity |ζ(λ) | z samples point across the entire domain [0, π] of the sin 2 , over which sin 2 has a mean value of 1/2. In our setup z = 17 cm, so we can average hI w d ðzÞi l ¼