Photonic band structure design using persistent homology

The machine learning technique of persistent homology classifies complex systems or datasets by computing their topological features over a range of characteristic scales. There is growing interest in applying persistent homology to characterize physical systems such as spin models and multiqubit entangled states. Here we propose persistent homology as a tool for characterizing and optimizing band structures of periodic photonic media. Using the honeycomb photonic lattice Haldane model as an example, we show how persistent homology is able to reliably classify a variety of band structures falling outside the usual paradigms of topological band theory, including"moat band"and multi-valley dispersion relations, and thereby control the properties of quantum emitters embedded in the lattice. Our method is promising for the automated design of more complex systems such as photonic crystals and Moire superlattices.


I. INTRODUCTION
There is growing interest in the design of synthetic gauge fields and topological band structures for light, motivated both by their exotic behaviour and potential applications as disorder-robust photonic devices [1,2]. The most common approach has been to emulate well-known topological phases from condensed matter physics [3]. However, recent advances in synthetic dimensions in photonic systems provide platforms for observing for the first time topological effects in higher dimensions [4][5][6]. Given the high degree of control and flexibility we now have over the design of photonic band structures and synthetic gauge fields, with many different degrees of freedom, a natural question is whether there exists interesting and useful topological phenomena falling outside the well-established paradigm of topological band theory, which is concerned with the "shape" of Bloch wave eigenstates, i.e. whether they exhibit vortices or gauge discontinuities within the Brillouin zone, independent of the energy dispersion.
Controlling the shape of the energy dispersion landscape can be equally important. For example, the shape of a medium's isofrequency contours or surfaces dictates the farfield radiation profile of localized emitters [7], and flat dispersion relations are highly desirable for realizing strong lightmatter interactions [8]. Both of these examples are based on the band structure at a specific energy. However, perfectly flat Bloch bands are idealizations never achievable in practice, and real systems inevitably exhibit losses, meaning there will be a finite energy resolution. These considerations motivate the use of "fuzzy" tools for characterizing the topology and abstract "shapes" of photonic media.
The characterization of topological properties in the presence of finite resolution or noise can be carried out using persistent homology, which is a tool from the field of topological data analysis (TDA) [9][10][11]. Persistent homology characterizes the "shape" and connectivity of high-dimensional data sets by studying how their topological properties change as a characteristic distance scale is varied. Several studies have started to apply persistent homology and other techniques from TDA to physics problems [12][13][14][15][16][17]. For example, persistent homology has recently been used to characterize entan-glement of multiqubit states [13], identify different dynamical regimes of the Bose-Hubbard model [14], and recognize topological phase transitions in models of interacting spins [15][16][17]. To apply persistent homology, the important first step is to choose a suitable distance metric to characterize the data.
Our aim in this study is introduce persistent homology to photonics and show how it may be useful for optimization of synthetic gauge fields and discovery of novel classes of photonic band structures which merit more detailed studies. As an example, we consider the characterization of the Haldane model, a tight binding model describing a honeycomb array of waveguides or resonators in the presence of a synthetic gauge field [18][19][20][21]. We show how to use persistent homology to characterize the shape of its isoenergy contours and the connectivity of the Bloch function eigenstates in the vicinity of this energy, described by the quantum distance [22][23][24][25][26][27][28]. As an application of this approach, we will show how these two features can be used to optimize the lifetimes and radiation profiles of embedded quantum emitters.
The outline of this article is as follows: In Sec. II we provide a brief introduction to the technique of persistent homology and discuss how it may be applied to band structure optimization by using suitable distance metrics. Sec. III applies persistent homology to analyze the band structure of the Haldane model, identifying parameter ranges exhibiting "moat band" [29,30] and multivalley dispersion relations [31,32]. In Sec. IV we show how the considered distance metrics can be used to predict properties of quantum emitters embedded in the lattice. Sec. V concludes with discussion of possible future directions and applications of persistent homology in photonics.

A. A brief introduction to persistent homology
The following discussion is quite condensed and not intended to be a comprehensive introduction to persistent homology. For a more pedagogical and complete introduction we recommend the excellent Ref. [12]  other recent articles applying persistent homology to characterize quantum and condensed matter systems [13][14][15][16]. Persistent homology is a technique used to characterize sets of data {x 1 , x 2 , ..., x N } equipped with some measure of distance between them d(x i , x j ) ≥ 0. For example, the data x i could correspond to particle positions in real space characterized by the Euclidean metric d(x i , x j ) = |x i − x j |, or image data with x i encoding the locations and intensities of the individual pixels. Given the data, one can construct subsets known as k-simplices, where each k-simplex is a set of k + 1 of the points; 0-simplices are individual points {x i }, 1-simplices are edges connecting pairs of points {x i , x j }, 2-simplices are areas enclosed by three points {x i , x j , x k }, and so on, as illustrated in Fig. 1.
To characterize the shape of the data, we construct sets of simplices using the distance measure d to form simplical complexes. A k-complex is a closed set of simplices of dimension up to k. By closed, we mean that if a simplex belongs to the complex, then its edges are also in the complex. For example, if an edge {x i , x j } belongs to the complex, then its constituent points {x i }, {x j } also belong to the complex. In the case of Euclidean data a common choice is the Vietoris-Rips complex V ε , which includes all k-simplices whose constituent points all have a pairwise distance less than ε from each other.
By counting the number of simplices of each dimension in the complex we can determine its abstract "shape". For example, a complex consisting of 3 points and 2 edges describes a line, while a complex with 3 points and 3 edges forms a loop. Formally, the existence of topologically nontrivial structures of dimension k in the simplical complex is determined by computing its kth Betti number B k , which corresponds to the number of k-dimensional "holes" in the complex.
Adding a single vertex or edge to the complex can drastically alter its topological properties encoded by the Betti numbers. They key innovation provided by persistent homology which makes it useful for characterizing real world systems with noise is to study how the topology of the simplical complex varies with some characteristic scale ε, in effect studying the shape of the data over a range of scales, known as a filtration. Features which only persist for a small range of scales can be attributed to noise and discarded, while those that persist for a large range of scales are robust and meaningful.
There are several possible choice for the filtration, depending on the data to be analyzed. For the case of point cloud data, one can use the Vietoris-Rips complex V ε . However, this becomes extremely time consuming to compute for large numbers of data points. For point cloud data there are other, more efficient complexes such as the α complex [9,10]. In the case of image data (pixels on a regular grid), a standard choice of filtration is to construct simplical complexes out of all neighbouring pixels with intensities greater than the threshold value ε, thereby describing the number and shape of intensity maxima.
Given a dataset, distance metric, and choice of filtration, the persistent homology of the data can be computed using existing software libraries [33][34][35]. As long as the dataset is not too large, low dimensional topological features can be computed relatively quickly. However, computation time grows rapidly for high-dimensional features, generally requiring preprocessing to reduce the dataset size.

B. Characterization of energy bands using persistent homology
Let us now return to physics. In order to apply persistent homology to physical systems, we first need to identify the relevant data set {x i } and a suitable distance metric to characterize the data. Several different approaches have recently been proposed in the physics literature [12][13][14][15][16]. Here our aim is to use to persistent homology to characterize energy bands of periodic photonic media.
We will consider for simplicity Bloch wave spectra of photonic waveguide lattices obtained from the Schrödinger eigenvalue problem [1] whereĤ(k) is the lattice's Bloch Hamiltonian, k is the Bloch momentum (restricted to the first Brillouin zone), |u n (k) is the Bloch function, with the bra-ket notation encoding the internal (e.g. sublattice or spin) degrees of freedom, n is the band index, and ω n is the frequency of the mode. We note that Maxwell's equations for electromagnetic modes of periodic photonic media can be cast into a similar form to Eq. (1), making our approach equally applicable to the design and optimization of photonic crystals. Numerical solutions of Eq. (1) typically sample k on a regular grid in the Brillouin zone, providing a discrete set of data (the eigenvectors and eigenvalues) to which we can apply persistent homology. Keeping with the spirit of successful applications of machine learning to physics, our aim is to employ simple distance metrics with well-established importance that can be easily interpreted.
First, we consider classifying the energy eigenvalues of a band ω(k) using a similar approach to that used for characterizing image data. That is, we sample ω on a regular k grid (e.g., triangular, cubic, etc.) spanning the Brillouin zone. We then define a filter function f (k) = |ω(k) − ω 0 |, where ω 0 is some reference frequency of interest. Persistent homology can characterise the shape of the set of momenta k whose Bloch mode frequencies are within ε ω of the reference frequency ω 0 , i.e.
In particular, persistent homology can tell us the range of frequency resolutions ε ω over which topological features of the lattice's isoenergy contours persist. An important consideration is that we inevitably sample k on a discrete grid of length N; hence there will be noise induced by the discretization. To be precise, neighbouring grid points will have an energy difference δ ω ≈ v g · δ k, where v g = ∇ k ω is the local wave group velocity, and δ k = 2π/(N) is the grid spacing (using units such that the lattice period a = 1). Thus, topological features of the isoenergy surfaces are only meaningful if they persist over a range of energies larger than δ ω = max κ (v g )/δ k. As a crude estimate, if ∆ is the width of the band of interest, the group velocity will have a scale v g ∼ ∆ 2π . Therefore as a rule of thumb one should look for features with a minimum persistence of ∆/N to avoid discretization errors; in the following we will use δ ω = 2max κ (v g )/N as the persistence cutoff.
The calculation proceeds as follows: (1) Choose a reference energy ω 0 and grid size N. (2) Compute the persistent homology of the isoenergy contours, using the energy resolution ε ω as the filter function. (3) Discard features with persistence less than δ ω. (4) Sum the number of remaining features existing at each ε ω of interest to obtain the Betti numbers, which characterize the shape of the isoenergy surface at that resolution, i.e. its connectivity and number of holes.
As an example of this procedure, we consider the energy band of a simple square lattice shown in Fig. 2(a), described by the dispersion relation where J 1 and J 2 are the nearest and next-nearest neighbour hopping strengths, respectively. We take the reference energy to be ω 0 = 0, close to the band centre, and J 2 = 0.2J 1 . The dispersion ω(k) is shown in Fig. 2(b).
Applying persistent homology to this two-dimensional dispersion "image", we obtain the energy resolution-dependent Betti numbers B 0,1 (ε ω ) shown in Fig. 2(c). For small ε there is a sharp peak in the number of distinct clusters B 0 , which is induced by discreteness of the momentum grid and highly sensitive to the chosen grid size N. From this plot we see that thresholding, i.e. consideration of features only with a minimum persistence, is essential to obtaining meaningful Betti numbers.
Applying the minimum persistence threshold in Fig. 2(d), we obtain two "phases" with different 1st Betti numbers B 1 , which we counts the numbers of loops (or holes) in the data. Fig. 2(e) illustrates the first phase, considering modes with energies |ω| < 2J 1 , where 2J 1 is the chosen energy resolution. B 0 = 1 indicates the modes form a single connected cluster, as can be readily verified. B 1 = 3 means that the image hosts three inequivalent non-contractible loops, as shown. We note that one needs to be careful with counting non-contractible loops here, due to the periodic boundary conditions of the Brillouin zone.
As the energy cutoff ε ω is increased, the band minimum at ω ≈ −3J 1 becomes included in the filtration, destroying one of the holes and resulting in the image shown in Fig. 2(f). The loss of a hole corresponds to B 1 decreasing by 1. Note that that the removal of the remaining hole in the centre of the image does not affect B 1 ; due to the Brillouin zone's periodic boundary conditions, any loop encircling this hole can be continuously deformed to a point by first expanding it to the Brillouin zone edges at k x = ±π to create a pair of loops, and then annihilating the pair at the k y = ±π Brillouin zone edge.

C. Characterization of Bloch functions
Persistent homology allows us to also classify more complicated objects than energy band structure "images". In particular, topological phases of light and matter are associated with non-trivial topological properties of the Bloch function eigenstates |u n (k) , which are usually computed using the Berry connection or curvature [1]. Here we will consider the abstract shape of the Bloch functions characterized by the quan-tum distance d, The quantum distance between Bloch functions at infinitesimally-separated momenta corresponds to the quantum metric tensorĝ, with components The study of interesting observables related to properties of the quantum distance and quantum metric currently attracts broad interest [22][23][24][25][26][27][28]. For example, in graphene the large quantum distance between counter-propagating modes at the Dirac cones is responsible for the suppression of backscattering in the presence of long range (valley-conserving) disorder [22]; a small quantum distance means there can be efficient scattering between the two states by scalar (e.g. spinpreserving) disorder. We use the quantum distance to define a graph, with each vertex corresponding to a Bloch function with some momentum k. To apply persistent homology, we define two points k, k to be connected when d(k, k ) < ε d . The Betti numbers of the resulting graph can then be used to characterize the "shape" of the eigenstates, such as the existence of loops or multiple disconnected clusters. One can either consider how the number of features changes as a function of ε d , or consider the features at a particular ε d .
As with the energy eigenvalues, sampling of the Brillouin zone on a discrete grid will induce discretization noise for small quantum distance cutoffs ε d . For neighbouring grid points, d(k, k + δ k) ≈ δ k ·ĝ(k) · δ k. Therefore to avoid discretization errors one should consider features that persist for 2π N √ĝ . We note that under this definition, the existence of nontrivial shapes of the energy eigenvalues is a necessary (but not sufficient) condition for the Bloch functions to have nontrivial shapes, since the Bloch functions are typically continuous functions of k. Furthermore, while the energy eigenvalues share the periodicity of the reciprocal lattice, i.e. ω n (k) = ω n (k + G) for any reciprocal lattice vector G, the quantum distance generally does not share this periodicity. Thus, d(k, k + G) can be nonzero. In the following we will compute topological features of the Bloch functions in the first Brillouin zone.
As an example, we consider the connectivity of the Bloch function eigenstates in the one-dimensional Su-Schrieffer-Heeger (SSH) model shown in Fig. 3, which is described by the Bloch Hamiltonian [1] H(k) = (J 1 + J 2 cos k)σ x + J 2 sin kσ y , where J 1,2 are coupling strengths andσ i are Pauli matrices. The eigenvectors of Eq. (6) are θ (k) = Arg(J 1 + J 2 cos k − iJ 2 sin k), and lie on the unit circle, spanning the entire circle in the nontrivial phase J 2 > J 1 , limited to −π/2 < θ < π/2 in the trivial phase. The quantum distance between eigenstates at k and k can be evaluated exactly as attaining its maximal value of 1 only in the nontrivial phase. We sample the Brillouin zone on a uniform grid k = k n = 2πn/N, using N = 30 grid points. Figs. 3(c,d) show the Betti numbers obtained from the two phases, without using any persistence threshold. Similar to the case of the energy eigenvalue images, discretization-induced noise appears for small filtration values ε d as N-dependent peaks in B 0 . To systematically remove this noise, we introduce a persistence threshold based on the quantum metric, which in the SSH model takes the form Using a threshold of 2.5 ∑ n 2π N 2 g kk (k n ) (i.e., proportional to the average quantum distance), we obtain smooth curves dis-tinguishing the trivial and nontrivial phases in Fig. 3(e,f). In particular, the winding of the Bloch functions in the nontrivial phase is characterized by B 1 = 1, which persists over a large range of scales ε d . We note that as ε d → 1 each point becomes connected to all others, destroying the hole and resulting in B 1 = 0. Thus, to characterize the shape of the Bloch functions one should use an intermediate value of ε d .

III. CLASSIFYING BAND MINIMA IN THE HALDANE MODEL
We now apply the above measures to the classification of the Haldane model shown in Fig. 4(a) [18][19][20][21], described by the Bloch wave Hamiltonian 3), and δ 3 =  3) are the lattice vectors, J 1 , J 2 ≥ 0 are nearest-and next-nearestneighbour hopping strengths, M is the sublattice detuning, and φ parameterizes the staggered next-nearest-neighbour hopping phases.
A previous study showed that when M = φ = 0, the Haldane model can exhibit "moat bands", i.e. degenerate band minima forming closed lines in the Brillouin zone [29]. The existence condition for moat bands when M = φ = 0 was obtained by rewritingĤ(k) asĤ which implies its energy eigenvalues take the form ω 1,2 = . It follows that the eigenvalues exhibit a degenerate minimum along the contour T 2 x + T 2 y = J 1 /(2J 2 ). Since T 2 x + T 2 y ≤ 3, moat bands can only exist when the next-nearest neighbour coupling strength is sufficiently strong, i.e. J 2 > J 1 /6. When a moat band is formed the ground state of hard-core bosons at low densities spontaneously breaks time-reversal symmetry [29].
When M, φ = 0 the simple condition for the existence of a moat band no longer holds. We apply the persistent homology approach introduced in the previous Section to systematically determine the fate of the moat band for M, φ = 0. To do so, we set the reference energy as ω 0 = min k [ω 1 (k)], considering the shape of the lowest energy eigenvalues ofĤ(k). We sample the Brillouin zone on a grid of N × N = 30 × 30 points and consider features at a fixed frequency resolution ε ω = 0.1J 1 . First, when M = φ = 0 we observe that the persistent homology approach accurately reproduces the features of the energy minima obtained in Ref. [29]. Note that the threshold for the emergence of a single loop is slightly larger than the analytical threshold J 1 = J 1 /6, due to the finite frequency resolution we consider, which results in a filled disk for small loop radii. Fig. 4(b) shows an example of a single-looped energy minimum (B 1 = 1), which persists for nonzero φ or M. As J 2 is increased the loop expands, generating additional loops once it intersects the edges of the Brillouin zone, characterized by Betti numbers B 1 = 4 [ Fig. 4(c)] or 2 [ Fig. 4(e)], depending on the parameter values. For large J 2 the single connected contour splits into two isolated ring-shaped energy minima centred at the two K points. The differing behaviour under increasing φ and M is noteworthy; increasing φ always results in the destruction of the loops and their replacement with a pair of isolated, point-like energy minima [ Fig. 4(d)], whereas the loops can persist in the presence of moderate M.
Interestingly, by fine-tuning a pair of parameters (J 2 , φ ) or (J 2 , M), one can even obtain a "phase transition point" at which multiple Betti numbers B 0,1 intersect [ Fig. 4(f)]. At  M (b,d), obtained by applying persistent homology to the energy eigenvalues with a resolution ε ω = 0.1J 1 . Persistent homology captures the emergence and changes in the moat band minimum known for φ = 0, and demonstrates its persistence for nonzero φ . For large φ or M and J 2 /J 1 the moat is destroyed and replaced with a pair of isolated minima. this point it is possible to obtain a variety of shapes of the low energy modes by applying small perturbations.
It is also noteworthy that for the chosen energy resolution ε ω = 0.1J 1 , the nontrivial low energy mode shapes emerge when the lower band is much flatter than the upper band, making the shapes difficult to directly identify from the images of the full band structure in Fig. 4. One could alternatively normalize the energy resolution by the width of the band of interest, or consider features at other energies. For example, if ω 0 is set to the band gap, the nonzero overlap between the upper and lower bands in Figs. 4(c,d,f) can generate multiple disconnected loops with B 0,1 > 1.
Next, in Fig. 6 we study the connectivity of the low energy Bloch functions characterized by the quantum distance. We consider features at ε d = 0.5, motivated by the SSH model example discussed in the previous Section. We observe that moderate φ and M have qualitatively different effects; moderate φ results in the creation of two disconnected clusters of eigenvectors, while under increasing M the Bloch functions only form a single cluster. To understand these differing behaviours, we note that at large J 2 the band minima are centred at the K points K ± = (±4π/3, 0). A long wavelength expansion of Eq. (11) yields the continuum Hamiltonian Here M and φ have qualitatively different effects on the Dirac mass M ∓ J 2 √ 3 sin φ . When φ = 0 the two valleys have the same mass, meaning that there is a small quantum distance between them, hence the low energy Bloch functions form a single connected cluster. On the other hand, φ creates two valleys with opposite masses, and hence eigenvectors residing on opposite sublattices with a large quantum distance between them, corresponding to two clusters. This behaviour is fully consistent with the differing topological properties of the bands characterized by the Chern number [1].

IV. APPLICATION TO DECAY OF QUANTUM EMITTERS
In this Section we will show how our chosen metrics can be used to predict the behaviour of interesting physical observables. Namely, we will demonstrate that the shape and connectivity of the energy band minima affect the lifetime of localized quantum emitters coupled to the lattice and resonant with the band edge, by imposing constraints on the effective emitter-Bloch wave coupling strength.
We consider the coupling of N e localized quantum emitters to the Haldane lattice, described by the Fano-Anderson model (h = 1) [36][37][38][39] whereĉ α annihilates a bound state with energy ω α and nonradiative decay rate γ α ,û n (k) annihilates a Bloch wave in band n with momentum k, and g α,n (k) is the emitter-Bloch wave coupling strength, g α,n (k) = J e ∑ r u n (k)|g α (r) e ik·r ≡ J e u n (k)|g α (k) , (17) which we have decomposed into an energy scale J e and a dimensionless part | u n (k)|g α (k) | = 1 − d α,n (k), where d α,n (k) is the quantum distance between the emitter and Bloch wave polarizations at momentum k. For simplicity, we focus on the single excitation subspace, in which the emitter and Bloch wave probability amplitudes evolve according to i∂ t c α (t) = (ω α − iγ α )c α + dkg α,n (k)u n (k,t), (18) i∂ t u n (k,t) = ω n (k)u n (k,t) + ∑ α g * α,n (k)c α (t).
We assume that only the emitters are excited at t = 0, i.e. u n (k, 0) = 0, and applying the weak coupling and Markov approximations (see e.g. Refs. [36,37]), we arrive at an effective evolution equation for the emitter amplitudes, where describes the emitters' coupling and self-energy mediated by the lattice, and we have dropped the summation over the band index n, as only the band energetically-closest to the emitters will be relevant in the weak coupling regime. Note that under the weak coupling approximation J e only sets the scale of ∆ α,β and is otherwise irrelevant. When the emitter frequencies are close to the lower band edge the integral in Eq. (21) will be dominated by modes with energies within the emitter linewidth |ω 1 (k) − ω β | < γ β , such that ∆ α,β becomes sensitive to the shape and connectivity of the eigenvalues and Bloch functions within this energy resolution. For example, when the emitter couplings g α,1 ≈ 0 at the high symmetry points of the Brillouin zone, the emitterlattice coupling will be strongly suppressed unless the band minimum forms a "moat" band.
We illustrate this in Fig. 7 for the case of a single emitter tuned to the band minimum with non-radiative decay rate γ e = 0.1J 1 and coupled to a 6-site plaquette to form a "giant atom" [38,39], described by the coupling Fig. 7(a,b) shows the radiative decay rate Γ = −Im[∆ e,e ] as a function of J 2 , M, and φ . The radiative decay rate is smallest when J 2 = 0, corresponding to a point-like band minumum at k = 0. The peaks in the decay rate coincide reasonably well with the formation of the moat bands. To show that this enhancement of Γ is sensitive to the shape of the low energy modes and not merely determined by the number of resonant modes, we also show in Fig. 7(c,d) the normalized decay rate Γ/N , where is the the |g 1 (k) -independent part of Eq. (21). The normalized decay rate is largest when there are multiple loops, i.e. B 1 = 2 or 4, which optimizes the trade-off between emitter-Bloch wave mode detunings ω 1 (k) − ω e and the emitter-Bloch wave coupling u n (k)|g e (k) .

V. CONCLUSION
We have shown how a method from the field of topological data analysis, namely persistent homology, can be used to identify novel classes of periodic lattices by systematically classifying the abstract shapes of their energy dispersion and Bloch functions. By using physically-motivated choices for the minimal persistence, we were able to eliminate discretization-induced noise in the topological features, thereby enabling robust high-throughput classification of photonic band structures. The approach is highly flexible and can be readily generalized to higher dimensional settings [4,5], complex lattices with large unit cells such as Moire superlattices [40,41], continuous media such as photonic crystals [42], and nonlinear or strongly interacting topological wave media [2].
For the simple tight binding models we considered the small parameter space could be covered through brute-force scanning. Persistent homology is compatible with gradientbased optimisation approaches, so an interesting direction for future research will be to use gradient-based methods to find specific energy band topologies of interest given constraints in the lattice design.
In contrast to some other recent studies employing persistent homology in physics [12][13][14][15][16][17], rather than focusing on the construction of persistence diagrams and their similarity measures, we focused on distance metrics directly related to physical observables of interest, such as the energy resolution. This allowed us to reliably separate robust features of the band structure from topological noise induced by numerical discretization. Future work may study and optimize the behaviour of other types of observables using persistent homology. We also hope this approach will also attract interest in other branches of physics beyond photonics, such as the design and classification of novel condensed matter and quan-tum systems.