Bloch-like waves in random-walk potentials based on supersymmetry

Bloch's theorem was a major milestone that established the principle of bandgaps in crystals. Although it was once believed that bandgaps could form only under conditions of periodicity and long-range correlations for Bloch's theorem, this restriction was disproven by the discoveries of amorphous media and quasicrystals. While network and liquid models have been suggested for the interpretation of Bloch-like waves in disordered media, these approaches based on searching for random networks with bandgaps have failed in the deterministic creation of bandgaps. Here we reveal a deterministic pathway to bandgaps in random-walk potentials by applying the notion of supersymmetry to the wave equation. Inspired by isospectrality, we follow a methodology in contrast to previous methods: we transform order into disorder while preserving bandgaps. Our approach enables the formation of bandgaps in extremely disordered potentials analogous to Brownian motion, and also allows the tuning of correlations while maintaining identical bandgaps, thereby creating a family of potentials with ‘Bloch-like eigenstates'.


Supplementary Note 1. The localization based on SUSY transformations
Supplementary Figure 1 shows examples of SUSY-transformed eigenstates (from the 0 th to the 3 rd modes, Mψ = {ik / k 0 + ∂ x ψ 0 (x) / [k 0 ψ 0 (x)]}·ψ) in crystals and quasicrystals. In the case of the crystal with eigenstates that have highly overlapping intensity profiles, the 'bound' profile of ψ 0 (x) decreases the spatial bandwidth of each eigenstate with the contribution of {∂ x ψ 0 (x) / [k 0 ψ 0 (x)]}·ψ. However, because the eigenstates in the quasicrystal are already spatially separated, in contrast to those in the crystal, the spatial modification by ψ 0 (x) occurs in a much more complex manner ( Supplementary Figs 1e-1h).
To quantify the localization of eigenstates in SUSY-transformed potentials, we introduce the definition of the effective width 1 w eff based on the inverse participation ratio, as where I(x) is the intensity of each eigenstate. Supplementary Figure 2 shows the change of effective width for the 30 lowest eigenvalue modes by serial SUSY transformations (arrows in Supplementary Fig. 2). As shown, the effective width decreases gradually in the crystal potentials, whereas no tendency is observed in the quasicrystal potentials. However, for both crystals and quasicrystals, the cases of localized eigenstates were found from serial SUSY transformations.

Supplementary Note 2. Bloch-like waves in discrete optical systems based on coupled mode theory
The notion of discrete optical systems has played a critical role in the design of optical devices. By interpreting the actual landscape of optical potentials (permittivity and permeability) as a network of optical elements (or 'photonic molecule' 2,3 composed of Extending the discussion of the main manuscript, in Supplementary Note 2, we derive Bloch-like wave families based on supersymmetry (SUSY) in 'discrete' optical systems by applying CMT to the notion of photonic molecules 2,3 with periodicity.

Hamiltonian equation of photonic molecules based on CMT
In a 1-dimensional (or 'nearest-neighbor' coupling) problem, the elementary equation of CMT for the k th photonic atom is as follows: where ξ is a spatial axis x (or a time axis t) in the spatial (or temporal) CMT, ψ k and ρ k are the field amplitude and the self-evolution term of the k th atom, respectively (ρ k is a wavevector in a spatial CMT or a resonant frequency in a temporal CMT), and κ is the coupling coefficient between atoms.
If the photonic molecule is Hermitian without magneto-optical effects, the coupling between photonic atoms is symmetric 4 as κ ij = κ ji .

Band analysis of crystalline photonic molecules
Consider the case of 1-dimensional crystalline photonic molecules satisfying the Hermiticity Here, we investigate the binary atomic distribution in the same way as in the main manuscript (CMT parameters of identical ρ k = ρ 0 for all k, and κ 2m+1 = κ a , κ 2m+2 = κ b for m = 1, 2, …). Supplementary Figure 3 shows the schematics of binary photonic molecules, each for the waveguide-based example with a spatial CMT model 8 and the resonator-based example with a temporal CMT model 9 .
Supplementary Figure 4 shows CMT-calculated eigenspectra as a function of modal number for the discrete system shown in Supplementary Fig. 3. Due to the binary arrangement, a bandgap is achieved around the self-evolution term ρ 0 = 1. Note that the construction of the bandgap can be understood in terms of the repulsion (determined by κ a ) of even and odd parity-modal bands (formed by κ b ) 14 . By increasing the number of photonic atoms N, the eigenspectrum approaches the continuous band with a well-known cosine form 13,14 , which is the feature of the lattice having periodic couplings (Supplementary Fig.   4a). The width of the bandgap is determined by the contrast between κ a and κ b , which is the same as other photonic crystals ( Supplementary Fig. 4b).
Supplementary Figure 6 shows the eigenspectrum of each SUSY-transformed photonic molecule. Identical to the case of continuous potentials in the main manuscript, all of the spectral information in the original eigenspectrum ( Supplementary Fig. 6a) are preserved during the series of SUSY transformations (Supplementary Figs 6b,6c), including the width and position of the bandgap and eigenbands. Note that the well-known cosine form of the eigenbands, which had been believed to originate from the periodic coupling 13 , is also reproduced perfectly with disordered photonic molecules. Due to the ground-state annihilation, the lowest part of the SUSY-transformed eigenspectrum has eigenvectors localized to decoupled atoms (black arrows).

SUSY-based random-walk photonic molecules: real space design
To design the real structure corresponding to the CMT parameters used in Supplementary Fig.   5, the position and self-evolution term (a wavevector in a spatial CMT, and a resonant frequency in a temporal CMT) of each photonic atom should be determined. While the selfevolution can be easily manipulated through the design of photonic atoms, the coupling is mainly determined by the interatomic distance. The coupling coefficient is generally obtained where Δε is the perturbation of permittivity distribution and e k is the normalized field pattern of the k th photonic atom. In the weak coupling regime, based on the evanescent field overlap, the coupling coefficient can be approximated as κ ij ~ c 1 ·exp(-c 2 ·Δx ij ) where c 1,2 are platformdependent constants and Δx ij is the distance between the i th and j th photonic atoms. Two unknown constants c 1,2 are determined when κ ij for two different distances are defined, and then from c 1,2, all of the coupling coefficients in Supplementary Fig. 5 can be converted to actual physical locations. The spatial distributions of the photonic molecules (obtained from Supplementary Fig. 5) are shown in Supplementary Fig. 7, presenting the spatially disordered potential shape after the SUSY transformations.
From the results in Supplementary Fig. 7, we can now calculate the correlation of the potential shapes by applying the Hurst exponent 18,19 . Supplementary Figure 8 20 . The matrix-based SUSY randomization can also be extended into other basis systems allowing discretization, such as tight-binding analysis, plane-wave expansion methods, and density functional theory in quantum mechanics.

Supplementary Note 3. The annihilation of eigenstates from 2D SUSY transformations
Starting from the potential with the form V o (x,y) = V ox (x) + V oy (y), it can be shown that the eigenstates of the 2D potential are combinations of the eigenstates from the 1D potentials V ox (x) and V oy (y). By using the separation of variables for the 2D Schrodinger-like equation with an eigenstate ψ(x,y) = φ(x)·ϕ(y), each brace should be a constant, which is one of the eigenvalues of the 1D Schrodinger-like equation, with the potential V ox (x) or V oy (y) (green symbols in Supplementary Fig. 9).
Following the discussion in the Methods in the main manuscript, the annihilation by 2D SUSY transformations occurs not only in the ground state but also in all of the excited states that share a common 1D ground-state profile. The example of this phenomenon is shown in Supplementary Fig. 9 for the x-axis SUSY transformation.