Photonic band structure design using persistent homology

There is growing interest in the design of synthetic gauge fields and topological band structures for light Haldane2008 ; Wang2009 ; Malkova2009 ; Khanikaev2013 ; Rechtsman2013 ; Hafezi2013 , motivated both by their exotic behaviour and potential applications as disorder-robust photonic devices topo_review ; topo_review2 . The most common approach has been to emulate well-known topological phases from condensed matter physics TI_review . However, recent advances in synthetic dimensions in photonic systems provide platforms for observing for the first time topological effects in higher dimensions synthetic1 ; synthetic2 ; synthetic3 . 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 far-field radiation profile of localized emitters hyperbolic_review , and flat dispersion relations are highly desirable for realizing strong light-matter interactions FB_review . 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) Ghrist2008 ; TDA_review ; TDA_review2 . 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 Murugan2019 ; Mengoni2019 ; Spitz2020 ; Tran2020 ; Olsthoorn2020 ; Cole_arxiv . For example, persistent homology has recently been used to characterize entanglement of multiqubit states Mengoni2019 , identify different dynamical regimes of the Bose-Hubbard model Spitz2020 , and recognize topological phase transitions in models of interacting spins Tran2020 ; Olsthoorn2020 ; Cole_arxiv . 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 a simple 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 Haldane0 ; Haldane1 ; Haldane2 ; Haldane3 . 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 Haldane_talk ; Kolodrubetz2013 ; Palumbo2018 ; Bleu2018 ; Rhim2020 ; Gianfrate2020 ; SCQ . 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” Sedrakyan2014 ; Sedrakyan2015 and multivalley dispersion relations multivalley ; multivalley2 . 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.

We now apply the above measures to the classification of the Haldane model shown in Fig. 4(a) Haldane0 ; Haldane1 ; Haldane2 ; Haldane3 , described by the Bloch wave Hamiltonian

	$\displaystyle\hat{H}(\bm{k})=$	$\displaystyle 2J_{2}\cos\phi\sum_{i=1}^{3}\cos(\bm{k}\cdot\bm{a}_{i})\hat{\sigma}_{0}+J_{1}\sum_{i=1}^{3}\left[\cos(\bm{k}\cdot\bm{\delta}_{i})\hat{\sigma}_{x}\right.$		(11)
		$\displaystyle\left.+\sin(\bm{k}\cdot\bm{\delta}_{i})\hat{\sigma}_{y}\right]+[M+2J_{2}\sin\phi\sum_{i=1}^{3}\sin(\bm{k}\cdot\bm{a}_{i})]\hat{\sigma}_{z},$

where $\bm{\delta}_{1}=(0,1/\sqrt{3})$, $\bm{\delta}_{2}=\frac{1}{2}(-1,-1/\sqrt{3})$, and $\bm{\delta}_{3}=\frac{1}{2}(1,-1/\sqrt{3})$ are displacements between neighbouring unit cells, $\bm{a}_{1}=(1,0)$, and $\bm{a}_{2}=\frac{1}{2}(1,\sqrt{3})$, and $\bm{a}_{3}=\frac{1}{2}(-1,\sqrt{3})$ are the lattice vectors, $J_{1},J_{2}\geq 0$ are nearest- and next-nearest-neighbour hopping strengths, $M$ is the sublattice detuning, and $\phi$ parameterizes the staggered next-nearest-neighbour hopping phases.

A previous study showed that when $M=\phi=0$, the Haldane model can exhibit “moat bands”, i.e. degenerate band minima forming closed lines in the Brillouin zone Sedrakyan2014 . The existence condition for moat bands when $M=\phi=0$ was obtained by rewriting $\hat{H}(\bm{k})$ as

		$\displaystyle\hat{H}(\bm{k})=J_{1}\hat{T}+J_{2}\hat{T}^{2},$		(12)
	$\displaystyle\hat{T}=T_{x}\hat{\sigma}_{x}+T_{y}\hat{\sigma}_{y}$	$\displaystyle=\sum_{i}\left[\cos(\bm{k}\cdot\bm{\delta}_{i})\hat{\sigma}_{x}+\sin(\bm{k}\cdot\bm{\delta}_{i})\hat{\sigma}_{y}\right],$

which implies its energy eigenvalues take the form $\omega_{1,2}=\pm J_{1}\sqrt{T_{x}^{2}+T_{y}^{2}}+J_{2}(T_{x}^{2}+T_{y}^{2})$. It follows that the eigenvalues exhibit a degenerate minimum along the contour $\sqrt{T_{x}^{2}+T_{y}^{2}}=J_{1}/(2J_{2})$. Since $\sqrt{T_{x}^{2}+T_{y}^{2}}\leq 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 Sedrakyan2014 .

When $M,\phi\neq 0$ the simple condition for the existence of a moat band no longer holds. We apply the persistent homology approach discussed in Sec. IIB to systematically determine the fate of the moat band for $M,\phi\neq 0$. To do so, we set the reference energy as $\omega_{0}=\mathrm{min}_{\bm{k}}\;[\omega_{1}(\bm{k})]$, considering the shape of the lowest energy eigenvalues of $\hat{H}(\bm{k})$. We sample the Brillouin zone on a grid of $N\times N=30\times 30$ points and consider features at a fixed frequency resolution $\epsilon_{\omega}=0.1J_{1}$. Figs. 4(b-f) show examples of the band structure in different phases identified using persistent homology, exhibiting different shapes of the lowest energy modes characterized by the Betti numbers plotted in Fig. 5.

First, when $M=\phi=0$ we observe that the persistent homology approach accurately reproduces the features of the energy minima obtained in Ref. Sedrakyan2014 . 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 $\phi$ 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 $\bm{K}$ points. The differing behaviour under increasing $\phi$ and $M$ is noteworthy; increasing $\phi$ 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},\phi)$ 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 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 $\epsilon_{\omega}=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 $\omega_{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 $\epsilon_{d}=0.5$, motivated by the SSH model example discussed in the previous Section. The Bloch function connectivity provides additional information about the shape of the band properties, complementary to the shape of the low energy eigenvalues in Fig. 5. Since the Bloch functions are continuous functions of the momentum $\bm{k}$, they can only form disconnected clusters if the isoenergy contours also form disconnected clusters. This explains why the $B_{0}=2$ states appearing in Fig. 6(a) are fully contained in the $B_{0}=2$ region of Fig. 5(a). Similarly for the case of $B_{1}$, isoenergy loops are neccessary for loops of the Bloch functions to appear.

Moderate $\phi$ and $M$ have qualitatively different effects on the shape of the Bloch functions; moderate $\phi$ 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 $\bm{K}$ points $\bm{K}_{\pm}=(\pm 4\pi/3,0)$. A long wavelength expansion of Eq. (11) yields the continuum Hamiltonian

	$\displaystyle\hat{H}(\bm{K}_{\pm}+\bm{p})$	$\displaystyle=-3J_{2}\cos\phi\;\hat{\sigma}_{0}+\frac{\sqrt{3}}{2}J_{1}(\mp p_{x}\hat{\sigma}_{x}-p_{y}\hat{\sigma}_{y})$
		$\displaystyle+(M-J_{2}\sin\phi[\pm\sqrt{3}+p_{x}+\sqrt{3}p_{y}])\hat{\sigma}_{z},$		(13)

Here $M$ and $\phi$ have qualitatively different effects on the Dirac mass $M\mp J_{2}\sqrt{3}\sin\phi$. When $\phi=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, $\phi$ creates two valleys with opposite masses, and hence eigenvectors residing on opposite sublattices with a large quantum distance between them, corresponding to two clusters.

Even though the correspondence between the Betti numbers and the more familiar band topological invariants such as the Chern number $C$ is not exact, the latter may be used to guide the design of Bloch function clusters and loops. For example, in trivial $C=0$ phases the Bloch functions of a band can be continuously deformed to a point, such that the quantum distance vanishes topo_review . Hence the Bloch functions are likely to form a single cluster at any energy resolution. On the other hand, such a deformation is impossible in non-trivial $C\neq 0$ phases and there always exists Bloch functions separated by a large quantum distance, which enables the formation of distinct clusters. In the case of the Haldane model, phases with non-trivial Chern number $C$ occur when the two valleys have opposite signs of the Dirac mass.

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 ($\hbar=1$) Longhi2018 ; Longhi2019 ; giant_atom ; giant_atom2

	$\displaystyle\hat{H}$	$\displaystyle=\hat{H}_{e}+\hat{H}_{L}+\hat{H}_{\mathrm{int}},\quad\hat{H}_{e}=\sum_{\alpha=1}^{N_{e}}(\omega_{\alpha}-i\gamma_{\alpha})\hat{c}_{\alpha}^{\dagger}\hat{c}_{\alpha},$		(14)
	$\displaystyle\hat{H}_{L}$	$\displaystyle=\sum_{n=1}^{N_{b}}\int d\bm{k}\omega_{n}(\bm{k})\hat{u}_{n}(\bm{k})\hat{u}_{n}^{\dagger}(\bm{k}),$		(15)
	$\displaystyle\hat{H}_{\mathrm{int}}$	$\displaystyle=\sum_{\alpha,n}\int d\bm{k}g_{\alpha,n}(\bm{k})\hat{c}_{\alpha}^{\dagger}\hat{u}_{n}(\bm{k})+\mathrm{h.c.},$		(16)

where $\hat{c}_{\alpha}$ annihilates a bound state with energy $\omega_{\alpha}$ and non-radiative decay rate $\gamma_{\alpha}$, $\hat{u}_{n}(\bm{k})$ annihilates a Bloch wave in band $n$ with momentum $\bm{k}$, and $g_{\alpha,n}(\bm{k})$ is the emitter-Bloch wave coupling strength,

$$g_{\alpha,n}(\bm{k})=J_{e}\sum_{\bm{r}}\braket{u_{n}(\bm{k})}{g_{\alpha}(\bm{r})}e^{i\bm{k}\cdot\bm{r}}\equiv J_{e}\braket{u_{n}(\bm{k})}{g_{\alpha}(\bm{k})},$$

(17)

which we have decomposed into an energy scale $J_{e}$ and a dimensionless part $|\braket{u_{n}(\bm{k})}{g_{\alpha}(\bm{k})}|=\sqrt{1-d_{\alpha,n}(\bm{k})}$, where $d_{\alpha,n}(\bm{k})$ is the quantum distance between the emitter and Bloch wave polarizations at momentum $\bm{k}$. For simplicity, we focus on the single excitation subspace, in which the emitter and Bloch wave probability amplitudes evolve according to

	$\displaystyle i\partial_{t}c_{\alpha}(t)$	$\displaystyle=(\omega_{\alpha}-i\gamma_{\alpha})c_{\alpha}+\int d\bm{k}g_{\alpha,n}(\bm{k})u_{n}(\bm{k},t),$		(18)
	$\displaystyle i\partial_{t}u_{n}(\bm{k},t)$	$\displaystyle=\omega_{n}(\bm{k})u_{n}(\bm{k},t)+\sum_{\alpha}g_{\alpha,n}^{*}(\bm{k})c_{\alpha}(t).$		(19)

We assume that only the emitters are excited at $t=0$, i.e. $u_{n}(\bm{k},0)=0$, and applying the weak coupling and Markov approximations (see e.g. Refs. Longhi2018 ; Longhi2019 ), we arrive at an effective evolution equation for the emitter amplitudes,

$$i\partial_{t}c_{\alpha}(t)=\sum_{\beta=1}^{N_{e}}\mathcal{H}_{\alpha,\beta}c_{\beta}(t),$$

(20)

where $\mathcal{H}_{\alpha,\beta}=(\omega_{\alpha}-i\gamma_{\alpha})\delta_{\alpha,\beta}+\Delta_{\alpha,\beta}$,

$$\Delta_{\alpha,\beta}=\int d\bm{k}\frac{g_{\alpha,1}(\bm{k})g_{\beta,1}^{*}(\bm{k})}{\omega_{1}(\bm{k})-\omega_{\beta}+i\gamma_{\beta}}$$

(21)

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 $\Delta_{\alpha,\beta}$ 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 $|\omega_{1}(\bm{k})-\omega_{\beta}|<\gamma_{\beta}$, such that $\Delta_{\alpha,\beta}$ becomes sensitive to the shape and connectivity of the eigenvalues and Bloch functions within this energy resolution. For example, when the emitter couplings $g_{\alpha,1}\approx 0$ at the high symmetry points of the Brillouin zone, the emitter-lattice 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 $\gamma_{e}=0.1J_{1}$ and coupled to a 6-site plaquette to form a “giant atom” giant_atom ; giant_atom2 , described by the coupling

$$\ket{g_{e,e}(\bm{k})}=\left(\begin{array}[]{c}1+e^{-i\bm{k}\cdot\bm{a}_{1}}+e^{i\bm{k}\cdot\bm{a}_{3}}\\ e^{i\bm{k}\cdot\bm{\delta}_{1}}[1+e^{-i\bm{k}\cdot\bm{a}_{1}}+e^{-i\bm{k}\cdot\bm{a}_{2}}]\end{array}\right).$$

(22)

Fig. 7(a,b) shows the radiative decay rate $\Gamma=-\mathrm{Im}[\Delta_{e,e}]$ as a function of $J_{2},M$, and $\phi$. The radiative decay rate is smallest when $J_{2}=0$, corresponding to a point-like band minumum at $\bm{k}=0$. The peaks in the decay rate coincide reasonably well with the formation of the moat bands. To show that this enhancement of $\Gamma$ 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 $\Gamma/\mathcal{N}$, where

$$\mathcal{N}=J_{e}^{2}\int d\bm{k}\left|\frac{\omega_{1}(\bm{k})-\omega_{1}-i\gamma_{\beta}}{(\omega_{1}(\bm{k})-\omega_{1})^{2}+\gamma_{1}^{2}}\right|.$$

(23)

is the the $\ket{g_{1}(\bm{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 $\omega_{1}(\bm{k})-\omega_{e}$ and the emitter-Bloch wave coupling $\braket{u_{n}(\bm{k})}{g_{e}(\bm{k})}$.

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 reliable high-throughput classification of photonic band structures. The approach may be used to save many graduate student-hours, for example if it is necessary to compute and check thousands of band structures to find an optimal lattice design subject to certain constraints. Moreoever, the method is highly flexible and can be readily generalized to higher dimensional settings synthetic1 ; synthetic2 , complex lattices with large unit cells or overlapping bands such as Moiré superlattices TBG ; Moire_photonic , continuous media such as photonic crystals photonic_crystal , and nonlinear or strongly interacting topological wave media topo_review2 .

For the simple tight binding models we considered, the small parameter space could be covered through brute-force scanning. Persistent homology is compatible with gradient-based 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 Murugan2019 ; Mengoni2019 ; Spitz2020 ; Tran2020 ; Olsthoorn2020 ; Cole_arxiv , 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 optimize the behaviour of other types of observables using persistent homology, study connections between the Betti numbers and the usual invariants of topological band theory such as the Chern number, and explore applications of persistent homology to other areas of photonics such as imaging. We also hope this approach will also attract interest in branches of physics beyond photonics, such as the design and classification of novel condensed matter and quantum systems.

This research was supported by the National Research Foundation, Prime Ministers Office, Singapore, the Ministry of Education, Singapore under the Research Centres of Excellence programme, and the Polisimulator project co-financed by Greece and the EU Regional Development Fund.

The data that support the findings of this study are available from the corresponding author upon reasonable request.