Topological states in disordered arrays of dielectric nanoparticles

We study the interplay between disorder and topology for the localized edge states of light in topological zigzag arrays of resonant dielectric nanoparticles. We characterize topological properties by the winding number that depends on both zigzag angle and spacing between nanoparticles in the array. For equal-spacing arrays, the system may have two values of the winding number $\nu=0$ or $1$, and it demonstrates localization at the edges even in the presence of disorder, being consistent with experimental observations for finite-length nanodisk structures. For staggered-spacing arrays, the system possesses richer topological phases characterized by the winding numbers $\nu=0$, $1$ or $2$, which depend on the averaged zigzag angle and disorder strength. In a sharp contrast to the equal-spacing zigzag arrays, staggered-spacing arrays reveal two types of topological phase transitions induced by the angle disorder, (i) $\nu = 0 \leftrightarrow \nu = 1$ and (ii) $\nu = 1 \leftrightarrow \nu = 2$. More importantly, the spectrum of staggered-spacing arrays may remain gapped even in the case of a strong disorder.


I. INTRODUCTION
Topological photonics emerged as a novel platform to realize robust optical circuitry protected against disorder [1]. The initial study of topological effects in photonics was inspired largely by direct analogies with similar effects in condensed matter physics for topologically nontrivial energy bands of electrons. In photonics, the concept of topological phases provides non-traditional approaches in the search for innovative designs of advanced photonic devices as well as introduces novel physical effects and their applications [2]. Spatially localized topological edge states have been predicted and realized in a variety of photonic systems, including topological lasers [3], single-photon transport [4], and topological solitons [5]. The study of disorder in topological systems is very important because robustness against disorder is extensively employed as a measure of topologically protected features. Anderson localization is a well-known effect of disorder in condensed-matter physics [6], which originates from interference of coherent waves, and it has been extended from electronic to other waves such as electromagnetic waves [7,8], acoustic waves [9], and matter waves [10,11]. Strong spatial localization of light due to disorder can be realized in periodic or quasiperiodic dielectric structures, and it is represented by photonic band gaps. For the systems with nontrivial topological phases, the interplay of topology-induced spatial localization (i.e. topological edge states) and disorder-induced Anderson localization becomes a very important problem, especially when disorder is no longer weak. Here, for the first time to our knowledge, we address this important problem by studying the effects of disorder in topologically nontrivial zigzag arrays of coupled subwavelength optical resonators.
Many topological photonic systems which have been explored so far are based on waveguide geometries, and the size of their building blocks is therefore larger than the wavelength of light. A zigzag array of coupled nanoresonators has emerged as the first nanoscale system with nontrivial topological states of light [12], and its topology-driven effects have been observed experimentally via lasing from excitonic-polaritonic nanoscale cavities [3] and third-harmonic generation from nonlinear arrays of dielectric nanoresonators [13].
In this paper, we study the interplay between disorder and topology in the zigzag arrays of resonant dielectric nanoparticles. We consider four types of topological arrays with equidistant and non-equidistant spacing of particles, as illustrated in Fig. 1, and analyze their topological phases and disorder-induced topological phase transitions. We also present experimental results for Mieresonant dielectric nanoparticles which confirm the basic predictions of the effect of disorder on topological properties, however being limited by the length of arrays. Remarkably, for staggered-spacing arrays, we observe that the topological spectral gap can survive even in the presence of strong disorder, the effect that was never observed before for any topological system, to the best of our knowledge.
The paper is organized as follows. In Sec. II, we introduce our model of a zigzag array of dielectric nanoparticles. In Sec. III, we study the topological properties of an ideal (regular) array and its disordered modifications. In Sec. IV, we use the zero-energy localization method to probe the disorder-induced topological transitions. In Sec. V, we generalize our model of the topological zigzag array from the equal-spacing to staggered-spacing case. Finally, Sec. VI concludes the paper.

II. EQUAL-SPACING ZIGZAG ARRAYS
It was established that topologically nontrivial states may appear in zigzag arrays of nanoparticles with polarization-dependent interaction between electromagnetic modes [12][13][14][15]. Here, we follow those earlier predictions and consider the electric or magnetic polarizations of dielectric nanoparticle in the x − y plane, as shown in Fig. 2(a), coupled due to the dipole-dipole interaction. Such polarization-dependent interaction between two neighboring nanoparticles can be written in the form [15] where e (j,j ) and e (j,j ) ⊥ denote respectively the identity vectors parallel and perpendicular to the link vector r j − r j . V (j,j ) = V (j ,j) , and the coupling decays with the distance between the jth and j th nanoparticles due to the dipole-dipole interaction. It is reasonable to keep only the nearest-neighboring interaction [12][13][14][15].
(3) Here, a † j,ν and a j,ν are the creation and annihilation operators for a ν mode at the jth nanoparticle, and ψ is the zigzag angle, as indicated in Fig. 2(b). E 0 is the on-site potential which causes an overall energy shift and hence is neglected. The ratio t /t ⊥ depends on the interaction, e.g., t /t ⊥ = −2 for the dipole-dipole interaction, and t /t ⊥ = −4 for the quadrupole-quadrupole interaction. Here, we consider the case of t /t ⊥ = −2.

A. Ideal arrays
We start with the study of an ideal system without disorder. In this case, an infinite system is translationally invariant. Under periodic boundary condition, the Bloch Hamiltonian reads [14], with Here, h 0 = t 1 + e −iK , and h(K) ≡ (h x , h y , h z ) = ∆ 2 e −iK sin ψ, 0, ∆ 2 1 + e −iK cos ψ are components of an effective magnetic field with t = t + t ⊥ /2 and ∆ = t + t ⊥ . σ ≡ (σ x , σ y , σ z ) are the Pauli matrices.
According to the classification of topological insulators and superconductors, this system belongs to the so-called BDI class [16], meaning a certain "Cartan label" given to the corresponding symmetric space inÉlie Cartan's classification scheme dating back to 1926. The topological invariant of such systems is presented by the winding number, According to bulk-boundary correspondence, the nontrivial topological invariant in a bulk is always accompanied by the existence of topological edge states at the boundaries. Under open boundary condition, the energy spectrum with respect to the zigzag angle is shown in Fig. 2(c). The spectrum can be either gapless or gapped, where in-gap zero-energy states appear only in the gapped case, implying topologically nontrivial properties. To confirm this, we calculate the winding number as a function of the zigzag angle; see Figs. 2(d-e), respectively. Indeed, in the gapless case the winding number vanishes, while in the gapped case the winding number equals to 1. Our results are consistent with earlier studies [14], where topologically nontrivial phase was predicted for |ψ − π/2| < ψ thre with the threshold value We notice here that the original winding number (6) is ill-defined in the gapless case. When the energy bands are gapless, the trajectory of det Q in the complex plane crosses the original point where its angle is ambiguous; see Figs. 3(a-b). However, a tiny perturbation without breaking the chiral symmetry does not change the topological property. To fix this problem, in numerical calculations we need shifting detQ(K) and redefine ν = ν mod 2 to avoid this ambiguity [14]. For completeness, we also show the energy bands and the trajectory in the complex plane for the gapped case; see Figs. 3(c-d). det Q is always nonvanishing, and the winding number is 1.

B. Disordered arrays
Next, we study the effects of disorder on the zigzag array. When the zigzag angle between nanoparticles experiences random variations, we may expect the well-known Anderson localization to occur. At the same time, in the topologically nontrivial case characterized by the winding number ν = 1, we expect topological edge states, which are localize exponentially at one or two edges. Since the zigzag-angle disorder can be realized in experiment [13], we study the interplay between topology-induced edge localization and the Anderson localization. In the following, we keep t = 2, t ⊥ = −1 fixed and assume that disorder is applied to the zigzag angle, ψ = ψ 0 + δψ i , with averaged zigzag angle ψ 0 and δψ i = W i , where i is a random number distributed uniformly within [−π/2, π/2]. Such disorder affects only the coupling, and it will not break the chiral symmetry. As the range of zigzag angles ψ is limited, see Fig. 5(a), the strength of random angles δψ i should be limited within δψ ∈ [−π/2, π/2]. Thus, the strength of disorder in this context is restricted to The winding number (6) is obtained generally in the translation-invariant systems, and usually the straighforward definition fails for disordered systems. If disorder does no break the essential discrete symmetry, i.e. the chiral symmetry, the winding number can be still extracted via a real-space formula [17,18]. The calculation of real-space winding number relies on the homotopically equivalent flat band HamiltonianQ =P + −P − , in which are the projector onto the subspace with positive or negative energies, respectively. It can be found thatQ shares the same eigenstates with HamiltonianĤ [16,19]. For chiral-symmetric systems, Q can be composed by Here Γ A and Γ B are respectively the projectors onto the two sublattices, labelled by A and B. Then, the winding number ν can be calculated in real space via the formula [17,20,21] where X is the position operator. This formula remains valid even in the presence of disorder [17,18,20,22], and it has been used widely in the studies of disorder and topology [23][24][25][26][27].
In our system, the gapless regime results in ambiguity when we calculate the winding number. However, when the disorder is present, this quantity is shown to be valid even if the gap is closed [17,28]. We calculate the phase diagram of the winding number in the parameter space (ψ 0 , W ); see Fig. 4. As turns out, the topology of such a system is nontrivial for large area where W > 0. The appearance of a small area of the trivial phase can be understood as the finite-size effect, and it will tend to disappear in the thermodynamic limit L → ∞. This analysis suggests that the system may become topological when disorder of the zigzag angle is added. In the regular case (W = 0), the gapless region |ψ 0 − π/2| > ψ thre is topologically trivial. Once the disorder is added, it immediately changes to topologically nontrivial. Hence, the case W = 0 can be understood as a critical point of the phase transition for the region |ψ 0 − π/2| > ψ thre . When disorder is present, the system is pushed towards a topological phase.

IV. ZERO-ENERGY LOCALIZATION LENGTH
Above, we have found that our system is topologically nontrivial when the angle disorder is included. We notice that, for one-dimensional chiral-symmetric systems, a topological transition is usually accompanied by a divergence of the localization length at E = 0 [17,29,30]. The Lyapunov exponent, which is related to the inverse of the localization length, will vanish at the point of the topological transition and then it changes its sign. This is an important quantity to reveal the phase boundary.

A. Transfer matrix approach
First, we start with the transfer matrix method and employ it for our system. According to the theory of Anderson localization [6,31], a one-dimensional lattice system becomes an insulator as long as any even weak disorder appears, and the electron wavefunction decays exponentially characterizing a localized state, where ξ is called the localization length. By employing the transfer-matrix technique [32,33], we calculate the localization length of our one-dimensional system. The discrete Schrodinger equation for the eigenstates of our model reads where ϕ x and ϕ y are the wavefunctions of the x and y polarized modes at ith site, respectively. t i,i+1 σ (σ = x, y) is the coupling strength between i-th and (i + 1)-th particles, Equation (10) can be written in a matrix form is the so-called transfer matrix between i-th and (i + 1)th sites; we use the fact that t x t y − t 2 xy = t t ⊥ ), and also apply a constraint det(M i ) = 1. Thus, we are able to know the L-th wave function of the entire system by multiplying iteratively the transfer matrices, The Lyapunov exponent λ, that is the inverse of localization length λ ≡ 1/ξ, can be calculated from the eigenvalue problem, Here, we employ the efficient numerical method of Ref. [34,35] to calculate the localization length. The averaged zigzag angle is ψ0 = π/2, which is the same setup as Ref. [13]. (e) Edge localization ratio extracted from experimental observations as a function of disorder angle.

B. Topology and localization length
In Fig. 5(b), we show the phase diagram of the inverse localization length 1/ξ with respect to (W, ψ 0 ) at E = 0 In the gapped case, the zero-energy localization length is exactly the length of localization of the edge state.
Disorder will increase the localization length; see blue curve in Fig. 5(c). In the gapless case, the zero-energy localization length diverges, and disorder will result in a decrease of the localization length; see the red curve in Fig. 5(c). We find that divergence of the localization length occurs only in the gapless region |ψ 0 −π/2| > ψ thre when W = 0, and we do not find any divergence of the localization length when W > 0. This implies that the system remains in the same topological phase apart from the region (|ψ 0 − π/2| > ψ thre , W = 0). Thus, we may conclude that the topologically nontrivial system remains topological when W > 0.
As we find the winding number of the system remains nontrivial for W > 0, the system should present topologically nontrivial phenomena. One of the features in one-dimensional topological systems is the appearance of topological edge states. In Fig. 5(d), we show the density distribution of the eigenstates with E ≈ 0 under different strengths of disorder. We find that the topological edge states with averaged zigzag angle ψ 0 = π/2 are preserved as the value of W grows. Meanwhile, with an increase of disorder strength, there appear localized states for the zigzag angle ψ 0 = 0. For a weak disorder, the localization length of the edge state is so large that the edge state is difficult to distinguish for a small number of nanoparticles. However, the edge-localization effect will become distinct when the system size is increased.

C. Experimental results
We have verified our theoretical predictions in experiment. For experiment, we fabricate many different types of disordered zigzag arrays of nanodisks placed on a glass substrate, as shown in Figs. 6(a-d). The samples are obtained with the same platform and techniques as described earlier in Ref. [13]. To map topological states, we excite zigzag arrays by femtosecond laser pulses, and observe the generation of a third-harmonic signal, being a tool for mapping strong spatially-localization distributions of light in the arrays. We set different strengths of the angle disorder, as shown in the examples of Figs. 6(ad), and also in Supplementary Fig. 10. We analyze the edge localization via the third-harmonic field in the zigzag arrays composed of eleven nanodisks with the averaged zigzag angle ψ 0 = π/2 and different angle disorder. We find that the number of distinguishable edge states decreases when the disorder strength increases, see Fig. 5(e), however some of the edge is clearly observed for most of the cases. The reduction of the edge localization observed in experiment corresponds directly to the effect predicted numerically, and it can be understood as follows. The localization length increases with the strength of disorder W for the averaged zigzag angle ψ 0 = π/2, see Fig. 5(d). When the localization length of the edge state excesses the length of the array, we are not able to distinguish it from other states, similarly observed for the trivial case when ψ 0 = 0, see Fig. 5(d). For example, the localization length of the edge state is about ξ = 10 at W = 0.5 (this corresponds to the disorder angle δψ ≈ π/4). This is almost the same length as the extension of our eleven-disks nanoparticle arrays employed in experiment. Thus, while in general we confirm the basic prediction of the theory about the character of the topology-induced localization and its interplay with the disorder, it becomes difficult to recognize the edge states and their localization for stronger disorder when W > 0.5 with finite-extent arrays.

V. STAGGERED-COUPLING ZIGZAG ARRAYS
Above, we have found theoretically that disorderinduced topological transitions may happen only at the gapless case when |ψ 0 − π/2| > ψ thre . Here, we consider a generalized model and explore richer topological phases and predict disorder-induced topological phase transitions, for the first time to our knowledge. We notice that the coupling between particles strongly depends on the distance between nanoparticles [12], and in our model discussed above all nanoparticles are equally spaced. In this section, we consider a staggered-spacing array in which the relative distance between nanoparticles alternates for odd and even bonds, as shown in Fig. 7(a). In this case, the generalized Hamiltonian readŝ where ∆ (j) = 0 for even bond, and ∆ (j) = ∆ for odd bond.

A. Ideal arrays
We start with an ideal array when disorder is absent. Chiral symmetry is still preserved in the presence of the extra term, and the topological property is still characterized by the winding number Eq. (6). The phase diagram of the winding number is obtained by varying the averaged zigzag angle and disorder strength; see Fig. 7(b). We uncover three topologically different phases characterized by the winding numbers ν = 0, 1, 2. When ∆ = 0, the generalized model becomes the original model of the zigzag array discussed above.
In Figs. 7(c-d), we show two typical energy spectra under open boundary condition and the corresponding winding numbers. When there is no zero-energy edge mode, the topological phase is trivial, and therefore ν = 0. When zero-energy edge modes appear, the zero-energy edge states are degenerate, and their numbers n edge = 2, 4 correspond to ν = 1, 2, respectively. This reflects the well-known bulk-edge correspondence.

B. Disordered arrays
Next, we study the effects of disorder in zigzag angles. In Fig. 8 (see the top row). For each fixed ∆, we find a narrow yellow line in the diagram with an extremely large localization length, indicating a boundary of a topological phase transition. We calculate also the phase diagram of the winding number using Eq. (8), as shown in the bottom row in Fig. 8. The boundary of the phase transition coincides with the one shown in the diagram of the localization length. Besides, we notice that the region of trivial winding numbers shrinks to W = 0 when ∆ van-ishes. This is consistent with the case of ∆ = 0 discussed above, where the system is topologically nontrivial in the presence of angle disorder for any W > 0. As stated above, we observe three topologically distinct phases characterized by ν = 0, 1, 2. Besides the disorder-induced topological transition between ν = 0 and ν = 1, there is also the transition between ν = 1 and ν = 2 when ∆ < 0. The phase diagram of ∆ < 0 is similar to that of ∆ > 0, see the last column in Fig. 8. Next, we compare two examples of ∆ > 0 and ∆ < 0 in Fig. 9. We identify a disappearance of two-fold degenerate edge states in Fig. 9(a), corresponding to the transition from ν = 1 to ν = 0. In Fig. 9(b), we find that two-fold degenerate edge states are transformed into four-fold degenerate edge states, which characterises the transition from ν = 1 to ν = 2. A change in the number of zeroenergy edge states is a result of a topological transition, and the number of zero-energy edge states coincides with the winding number in Figs. 9(c-d).
When the deviation |∆| is small, the spectral gap will be closed when disorder becomes stronger. Intuitively, this happens because the disorder becomes dominating, and it "wipes away" the details of the staggered-spacing distribution. However, in both these cases of Fig. 9, we observe that the spectral gap does not vanish for a strong disorder. This result differs from a common belief that stronger disorder will force a spectral gap to close [17], or topological transitions occur due to disappearance of energy gaps. Our finding may help understanding better the nature of disorder-induced topological transitions.

VI. CONCLUSION
We have studied the effect of disorder on topological properties of zigzag arrays of nanoparticles with polarization-dependent interaction. For equal-spacing arrays, even in the presence of angle disorder, we have found that the system will remain in the topological regime supporting edge states, and this feature was found to be in agreement with experiments for finite-extend arrays of silicon nanoparticles with an introduced disorder. For staggered-spacing arrays, we have found richer topological phases and observed that the angle disorder may induce topological phase transitions. Remarkably, in this latter case the spectral gap remains open even for a strong disorder, and we believe this system provides the first example of a topologically nontrivial system with disordered parameters. We present in Fig. 10 results on experimentally accumulated statistics of field distributions in zig-zag chains with angular disorder varying from δψ = 5 • to δψ = 90 • . All chains are 11 disks long, and their configurations are schematically shown in left-side columns. The field distributions along the chains are studied via nonlinear imaging: the chains are excited by short-pulse, high peak power laser pulses at 1590 nm wavelength, and generation of third harmonic signal is imaged onto a camera. The corresponding distributions of the third harmonic are shown as false-colour images in the right-side columns. In the chosen experimental configuration, the topological state is expected to occur at the bottom left edge of the chains. Experimental realisations in which the formation of the edge state was not observed are highlighted with red colour. As the level of disorder increases, the localisation length of the edge states increases reaching and further surpassing the length of the 11-disk chains. Large localisation length in highly disordered zig-zags reduces the probability of edge states observation in finite-length chains.