Persistent homology analysis of a generalized Aubry-André-Harper model

Wave propagation in low dimensional systems reveals many surprises thanks to wave coherence and interference effects book ; disorder_review ; loc_review ; fb_review ; AL_experiment ; Kohlert2019 ; Ganeshan2015 ; Wang2020 ; Goblot2020 ; An2021 . Disorder gives rise to Anderson localization and a complete suppression of transport; on the other hand, quasiperiodic potentials can exhibit localization transitions and critical phases supporting self-similar or fractal behaviour. Adding quantum or nonlinear interactions leads to even richer phenomena including subdiffusive spreading, many body localization, and quantum scars nonlinear ; MBL ; scars . While distinguishing localized from delocalized phases is relatively straightforward, the identification of critical and non-ergodic phases remains challenging, particularly when the relevant order parameters are unknown. Standard approaches involve studying the finite size or time scaling of simple observables over a huge range of scales Anderson_review ; Roemer2011 .

An alternative to finite size scaling is to identify relevant order parameters directly using machine learning ML_review , an in particular topological machine learning approaches are attracting growing interest as a means of identifying nonlocal order parameters Ghrist2008 ; TDA_review ; TDA_review2 ; Murugan2019 ; TDA_review3 . Specifically, the technique of persistent homology computes topological features over a range of scales, recording the scales at which features appear and are destroyed. This information is typically summarized using a persistence diagram, which forms a topological “fingerprint” of the data. Persistent homology has been applied to various physical systems, including the characterization of amorphous phases of matter and the detection of quantum phase transitions Hiraoka2016 ; Mengoni2019 ; Olsthoorn2020 ; Tran2021 ; Leykam2021 ; Spitz2021 ; Cole2021 ; Tirelli2021 ; Sale2022 ; Park2022 .

Existing studies have largely used persistent homology as part of a larger machine learning pipeline, in which persistent homology identifies a large set of potentially-relevant order parameters, which are fed into clustering or dimensional reduction algorithms to learn the system’s phase diagram. One drawback of this approach is that it is harder to understand than conventional methods, raising the question of whether persistent homology-based machine learning can be trusted to unveil new phenomena and not merely reproduce known phases. More cases studies are needed to better understand the strengths and limitations of persistent homology-based methods for identifying phase transitions.

In this article we apply persistent homology to study phase transitions in the Aubry-André-Harper (AAH) model, a prototypical model of localization and criticality AAH1 ; AAH2 ; Zilberberg ; Ganeshan2013 ; Liu2015 ; Purkayastha2017 ; AAH_experiment . We use persistent homology to obtain a compact encoding of the shape of the eigenfunctions across a range of wavefunction intensity scales in the form of persistence diagrams. To distinguish the different phases of the model we compute two summary statistics of the persistence diagrams: the root-mean-square feature lifetime and the persistent entropy. We find that while these summary statistics generally behave similarly to the conventional statistic measures (inverse participation ratio and Shannon entropy of the full eigenstate profiles), and they can be better at distinguishing the localized, extended, and critical phases, especially for the propagation dynamics.

Surprisingly, we find that persistent homology also distinguishes ordered from disordered states, revealing hidden order of the low energy eigenstates that is not detected by the inverse participation ratio or Shannon entropy. We show how this hidden order arises by performing a long wavelength expansion of the tight binding Hamiltonian, demonstrating that the quasiperiodic potential becomes transparent to low energy modes. Finally, we carry out numerical simulations of wavepacket propagation to show how the persistent homology approach can also distinguish the different phases using wavepacket dynamics.

The article is structured as follows: Sec. II presents the phase diagram of the generalized Aubry-André-Harper model and analyzes its low energy continuum limit. Sec. III briefly reviews persistent homology and shows how it can be applied to study the shape of eigenfunctions. Sec. IV applies persistent homology to study wavepacket propagation dynamics. The article concludes with Sec. V.

We consider a generalized Aubry-André-Harper model introduced in Ref. Liu2015 , corresponding to a one-dimensional tight binding chain with quasiperiodically-modulated on-site energies and couplings illustrated in Fig. 1(a). The lattice Hamiltonian is

where $V_{1}$ and $V_{2}$ are the on-site and coupling quasiperiodic modulation strengths, respectively, $t$ is the coupling in the absence of modulation, and the modulation has a frequency $Q$ and phase $k$. We consider the same parameters as in Ref. Liu2015 : $Q=(1+\sqrt{5})\pi$, $t=1$, $0\leq V_{1}\leq 4$, $0\leq V_{2}\leq 2$, $k=0$, and a lattice size of $N=300$ sites. We apply periodic boundary conditions to the wavefunction, $\psi_{n+N}=\psi_{n}$. Note that when $N$ and $Q$ are incommensurate, the periodic boundary conditions result in an effective defect in the modulation between sites $1$ and $N$. Similar to the original AAH model, there are no mobility edges and all eigenstates of Eq. (1) are either localized, extended, or critical depending on the choice of parameters Liu2015 .

First we will analyze the eigenstates of Eq. (1) with energy $E$, $\hat{H}\phi_{n}=E\phi_{n}$, focusing on the properties of their probability density distribution $|\phi_{n}|^{2}$. Later, we will consider the propagation dynamics generated by Eq. (1), i.e. $i\partial_{z}\psi_{n}(z)=\hat{H}_{n}\psi_{n}(z)$, focusing on the case of light propagation in optical waveguide arrays, where one prepares an initial state $\psi_{n}(0)$ whose propagation along the longitudinal ($z$) direction is analogous to time evolution in the Schrödinger equation. We will apply a similar approach to analyze $z$-dependent intensity profiles $|\psi_{n}(z)|^{2}$; we will use probability density and intensity interchangeably in the following.

The localization of the normalized eigenstates ($\sum_{n}|\phi_{n}|^{2}=1$) is typically studied using the inverse participation ratio (IPR),

which measures the reciprocal of the number of strongly excited sites Anderson_review ; Roemer2011 . In the large $N$ limit we have $\mathrm{IPR}\propto N^{-d}$, where $d$ is the fractal dimension of the eigenstate. $d=0$ for localized eigenstates, e.g. exponential localization $|\phi_{n}|^{2}\propto\exp(-|n|/\xi)$ yields $\mathrm{IPR}\sim 1/\xi$; $d=1$ for extended eigenstates, e.g. Gaussian random matrix ensemble eigenstates satisfy $|\phi_{n}|^{2}=O(1/N)$, with $\mathrm{IPR}=3/N$. On the other hand, critical eigenstates have $0<d<1$.

Another measure of an eigenstate’s localization is the Shannon entropy of its probability distribution entropy_0 ; entropy_1 ; entropy_2 ,

which vanishes ($S=0$) when the probability density is concentrated at a single site, and attains its maximal value ($S=\ln N$) for delocalized states with modes distributed over all $N$ sites. For these simple profiles, one can obtain $S\sim-\ln(\mathrm{IPR})$.

Figure 1(b) plots the inverse participation ratio of the ground state of $\hat{H}$ as a function of $V_{1}$ and $V_{2}$. $\mathrm{IPR}$ clearly distinguishes the extended and localized phases, but shows a poor contrast between the extended and critical phases, requiring finite size scaling to unambiguously distinguish them. This is because $\mathrm{IPR}$ is mainly sensitive to the number of strongly-excited sites, but misses information contained in the weak tails of the eigenstates. This can be addressed by considering other moments of the probability distribution, i.e. the generalized inverse participation ratios $\mathrm{IPR}_{q}=\sum_{n}|\phi_{n}|^{2q}$, corresponding to the $q$th moment of the probability distribution, to obtain a spectrum of fractal exponents $d_{q}$ (the multifractal spectrum Anderson_review ; Roemer2011 ). We note that the Shannon entropy provides similar information to $\mathrm{IPR}_{1}$, and that the relation $S\sim-\ln(\mathrm{IPR})$ can be violated for eigenstates exhibiting multifractality.

Another limitation of the inverse participation ratio and Shannon entropy is that they are not strongly sensitive to the shape of the eigenstates. For example, both measures are invariant under permutations of the site positions and fail to distinguish a mode with a single, broad peak from one with multiple narrower peaks. While similarity measures for pairs of wavefunctions including the overlap correlation function Cuevas2007 and Kullback-Liebler divergence Luitz2015 are sensitive to changes in shape such as permutations in site positions, they cannot be applied to individual eigenstates and thus they provide a relative rather than absolute measure of shape. This motivates our investigation of whether persistent homology can provide a useful and more efficient characterization of the eigenstates.

Figure 1(c) illustrates the phase diagram obtained using persistent homology analysis, which we will explain in the following section. In addition to clearly distinguishing the three phases, Fig. 1(c) exhibits a local minimum in the eigenstate measure (persistent entropy $\mathcal{E}$) along the line $V_{1}=2V_{2}\cos(Q/2)$, corresponding to the emergence of ordered eigenstates. This order can be understood by transforming to reciprocal ($p$) space, $\psi_{n}=\sum_{p}\psi_{p}e^{ipn}$, which yields (assuming $k=0$):

When $V_{1}=2V_{2}\cos(Q/2)$ the quasiperiodic potential becomes transparent to the low-energy modes, e.g. the plane wave $\phi_{p}=\delta_{p,\pi}$ forms an exact eigenstate (the ground state).

To validate this intuition, Fig. 2 shows the three lowest energy eigenstates obtained numerically for $V_{1}=1=2V_{2}\cos(Q/2)$, revealing simple standing wave forms with intensities $\sim\sin^{2}(n/N)$, $\sin^{2}(2n/N)$, $\sin^{2}(3n/N)$. The standing wave patterns can be understood by noting that the periodic boundary conditions $\psi_{n+N}=\psi_{n}$ are incommensurate with the quasiperiodic modulation, thus there is a (weak) effective defect potential in the vicinity of the $N$th lattice site; the low energy modes form standing waves that have a minimal overlap with this defect.

More generally, whenever some momentum $p$ satisfies $V_{1}\leq-2V_{2}\cos(p\pm Q/2)$ there will be a decoupling of the momentum space eigenvalue problem into different sectors, which may explain the additional scar-like features seen in Figs. 1(b,c).

In the extended phase the transparency of the quasiperiodic potential for low energy modes is missed by the standard eigenstate measures such as $\mathrm{IPR}$ and $S$, which only measure the localization of the modes and not their degree of disorder. This suggests that persistent homology may be sensitive to different features of the eigenstates, providing an additional tool for understanding complex lattice models. On the other hand, in the critical phase scars are visible in both $\mathrm{IPR}$ and $S$, since these measures can distinguish critical eigenstates from the extended states for which the quasiperiodic potential is transparent.

Most previous applications of persistent homology to physics considered point cloud data equipped with a suitable metric, such as atomic positions and Euclidean distance for characterizing the structure of materials Hiraoka2016 . This is a natural approach for dealing with sets of discrete objects. On the other hand, for analyzing individual eigenfunctions it is less straightforward to construct point clouds, requiring some way to reduce the continuous data to a discrete set of points, for example by choosing strong local maxima of the probability density distribution using some threshold Spitz2021 .

An approach better-suited to probability density or intensity data is the sublevel set filtration, which characterizes the local maxima and minima of an image or intensity profile survey . Specifically, one considers the number of distinct low intensity clusters of the profile, i.e. the number of connected regions formed by sites with $|\phi_{n}|^{2}<\epsilon$ as a function of $\epsilon$. As $\epsilon$ is continuously increased from zero, a new cluster will emerge whenever $\epsilon$ is tuned through a local minimum of the intensity. Neighbouring clusters will merge when $\epsilon$ crosses the local intensity maximum between them. Thus, the sublevel set filtration provides a stable way to characterize the critical points of an intensity or probability density profile. In effect, this allows one to count the number of distinct peaks in a profile over a range of intensity scales, with the lifetime (i.e. persistence) of each peak giving a measure of its significance. This is complementary to other approaches for studying localization in disordered systems such as fractal analysis, based on indirectly varying a characteristic intensity scale (associated with the eigenstate normalization) by changing the system size $N$ Anderson_review ; Roemer2011 .

For example, Fig. 3 presents a simple intensity profile and its resulting persistence diagram. Each point in the persistence diagram represents a distinct local minimum of the intensity profile, with the horizontal (birth scale) and vertical (death scale) coordinates corresponding to the intensity at the local minimum and closest local maximum, respectively. All points occur above the diagonal (indicated by the dashed line), because the intensity at the local minimum will always be less than that at the nearest local maximum. The distance of a point to the diagonal, i.e. the intensity contrast between the local minima and maxima, provides a measure of the significance of the feature. Here there are two long-lived features (prominent peaks), and two short-lived features created by small fluctuations of the profile at $x\approx 3.6$.

We now construct persistence diagrams for the ground states of Eq. (1) in each of its three phases. Fig. 4 shows the eigenstates and corresponding persistence diagrams. In the extended phase the probability density of the ground state is oscillatory with significant intensity throughout the entire lattice, corresponding to all features (except one) being created and destroyed within a small range of intensities.

With the increase of $V_{1}$ the generalized AAH model enters the localized phase, with the intensity localized to a very small interval, corresponding to only a single long-lived feature. Other features (induced by fluctuations in the eigenstate tails) have negligible lifetimes.

In the critical phase the probability density distribution is not completely concentrated in a certain interval, with multiple disconnected local maxima. This results in many long-lived features with a range of birth scales, consistent with the self-similar (fractal) nature of the critical phase.

The above sublevel set filtration approach of persistent homology summarizes all the critical points of an intensity profile, but cannot be directly used to quantify differences between the three phases. For example, how do we identify the transition from the extended and critical phases by inspecting the persistence of their features? One solution is to introduce a metric quantifying the distance between persistence diagrams, such as the persistence landscape landscape . An simpler alternative we will consider here is to compute summary statistics of the feature lifetime, reducing each persistence diagram to a single, more easily-interpretable number.

We will consider two summary statistics. The first is the $p$-norm of the feature lifetimes Tran2021

where $b$ is the intensity at which a feature appears (i.e. a local minimum), and $d$ denotes the intensity at the local maximum at which the feature is destroyed. $d-b$ is thus the feature lifetime. For example, $\mathcal{P}_{\infty}$ measures the lifetime of the longest-lived feature in the persistence diagram $D$, $\mathcal{P}_{1}$ is the sum of feature lifetimes, and $\mathcal{P}_{2}$ is the root mean squared lifetime. Thus, the $p$-norms resemble the generalized inverse participation ratios.

The second summary is the persistent entropy PH_entropy_stability ,

where ${\mathcal{S}(D)}=\sum_{(b,d)\in D}{|d-b|}$ is the sum of feature lifetimes in the persistence diagram D. The persistent entropy quantifies the non-uniformity of the feature lifetimes $d-b$; it is simply the Shannon entropy of the feature lifetime fractions.

Essentially, the standard approach for quantifying eigenstate localization is to compute summary statistics of their intensity distribution such as the inverse participation ratio or Shannon entropy. On the other hand, persistent homology computes and characterizes the eigenstates’ critical points (local maxima and minima). In effect, we compress the information contained in the eigenstate profiles. The point summaries not only capture the localization properties of the eigenstate, but also its level of disorder. Referring back to Fig. 1, we see that the persistent homology measures not only reproduce the known phase boundaries of the generalized AAH model, but also identify model parameters for which order is restored, i.e. along the line $V_{1}=2V_{2}\cos(Q/2)$.

Figure 5(a,b) shows the phase diagrams obtained using feature lifetime norm $\mathcal{P}_{2}$ [Eq. (5)] and the ground state’s Shannon entropy [Eq. (3)]. $\mathcal{P}_{2}$ is poor at distinguishing the extended and critical phases, similar to IPR. On the other hand, $S$ provides similar contrast to the persistent entropy, but does not detect the ordered line. This suggests that the main advantage of the persistent homology approach may be in distinguishing order from disorder. We illustrate this further by plotting all four eigenstate measures along 1D cuts of the full phase diagram in Fig. 5(c,d). While all measures can reproduce the correct phase boundaries, the persistent entropy provides stronger contrast between the critical and extended phases and can also distinguish order from disorder in the ground state, visible as a prominent dip in $\mathcal{E}$ in Fig. 5(d).

For a given eigenstate measure the contrast between the different phases will depend on the system size $N$. For example, the generalized inverse participation ratios will scale in the large $N$ limit as $\mathrm{IPR}_{q}\propto N^{d_{q}(1-q)}$, where $d_{q}$ is the fractal dimension of the $q$th eigenstate moment Anderson_review ; Roemer2011 . Thus, $d_{q}$ can be obtained by studying the scaling of an eigenstate measure with $N$. Fig. 6 plots the finite size scaling of the measures of the lowest energy eigenstate, taking the system size $N$ to be Fibonacci numbers and corresponding rational approximants of $Q$ to minimize finite size effects. We observe that $\mathrm{IPR}$ and $(\mathcal{P}_{2})^{2}$ and $S$ and $\mathcal{E}$ exhibit identical scaling exponents, up to the accuracy of the obtained power law fits.

Both the standard and persistent homology-based measures give similarly faithful estimates for the fractal dimensions of the eigenstates by taking the large $N$ limit. On other hand, based on the offsets between the fits obtained for the different phases, it appears that, for fixed $N$, the $\mathrm{IPR}$ (persistent entropy) performs better at distinguishing the extended and critical phases than P2 (Shannon entropy). Since both approaches yield the same fractal dimensions $d_{1,2}$ we anticipate their multifractal spectra will also coincide. A more detailed study of multifractal properties probed using persistent homology is an interesting direction for future research.

Next we study how persistent homology analysis can be used to characterize the propagation dynamics. For this, we consider the evolution of a localized (single waveguide) excitation, described by Schrödinger equation $i\partial_{z}\psi_{n}(z)=\hat{H}\psi_{n}(z)$, with $\psi_{n}(0)=\delta_{n,N/2}$, where $N$ is the total number of lattice sites. We compute IPR, $S$, $\mathcal{P}_{2}$, and $\mathcal{E}$ of the $z$-dependent intensity profiles $|\psi_{n}(z)|^{2}$ in Fig. 7. To obtain statistics and uncertainty estimates, we construct an ensemble by varying the parameter $k$ appearing in Eq. (1), which is equivalent to spatial translations of the input beam.

The standard approach towards distinguishing the phases of $\hat{H}$ using the propagation dynamics is to consider the scaling of IPR, i.e. using a power law fit $\mathrm{IPR}\propto z^{-\alpha}$, where the fitted exponent $\alpha$ can distinguish localized, diffusive, and ballistic transport regimes loc_review ; AAH_experiment . In the localized phase there is an initial transient expansion until the wavepacket width is comparable to the localization length $\xi$, after which spreading stops, corresponding to $\mathrm{IPR}\sim 1/\xi$ and $\alpha=0$. In delocalized phase a single site excitation will spread indefinitely, corresponding to ballistic expansion with $\mathrm{IPR}\propto z$, i.e. $\alpha=1$. In the critical phase, the transport is approximately diffusive, with $\alpha\approx 1/2$Longhi2021 .

Consistent with these expectations, in the localized phase in Fig. 7 (orange curves) the IPR does not decay to zero but saturates at a finite value, indicating localization. Similar saturation is seen in $\mathcal{P}_{2}$. Meanwhile, $S$ and $\mathcal{E}$ grow the fastest in the extended phase and saturate at a low values in the localized phase. However, it is noticeable that there is a subtle difference between Fig. 7(a,b) and Fig. 7(c,d) after a long propagation distance: the persistent homology measures show a quantitatively larger difference between the localized phase and the extended and critical phases.

We repeat this procedure, scanning over the model parameters $V_{1}$ and $V_{2}$ and computing IPR, $S$, $\mathcal{P}_{2}$ and $\mathcal{E}$ after the wavepacket has propagated a sufficiently large distance. Fig. 8 illustrates the phase diagrams reconstructed from the beam propagation simulations. We see that both the conventional measures (IPR, $S$) and the persistent homology measures ($\mathcal{P}_{2}$, $\mathcal{E}$) are able to distinguish the dynamics in the three phases. The most striking difference compared to Fig. 1 is that the ordered line $V_{1}=2V_{2}\cos(Q/2)$ is not detected by the propagation dynamics. The reason for this is that along this line, the quasiperiodic potential is only transparent to the low energy modes; the single site initial state excites a superposition of all the modes. Thus, the field becomes disordered after propagating a sufficiently long distance.

To determine whether the ordered line can be feasibly observed in experiment we construct a finite size Gaussian beam with width $w=10$ and alternating phase between neighboring sites, $\psi_{n}(0)=(-1)^{n}\exp[-(n-N/2)^{2}/(2w^{2})]$, shown in Fig. 9(a). Owing to the similarity between this profile and the ground state, in the extended phase this initial state mainly excites the low energy modes. Fig. 9(b) shows two intensity profiles obtained after propagation at ordered and disordered points in the extended phase. Thanks to the suppression of the high energy modes there is a clear difference in the output intensity profiles; the profile is smooth on the ordered line, whereas elsewhere in the extended phase rapid intensity fluctuations are generated. The Shannon entropy of the output intensity profiles shown in Fig. 9(c) does not distinguish these two profiles, whereas the persistent entropy in Fig. 9(d) correctly identifies the ordered line, showing that the persistent homology method may also be useful for distinguishing ordered from disordered propagation dynamics without requiring ensemble averages.

We have studied how the persistent homology analysis can be applied to characterize the eigenstates and propagation dynamics of tight binding models. Focusing on the example of the generalized Aubry-André-Harper model, we have shown how point summaries of the persistence diagrams can be used to distinguish its three phases (extended, localized, critical), as well as distinguish ordered from disordered regimes. This is despite the fact that the persistence diagrams form a highly compressed representation of the eigenstates, i.e. containing only information about their local maxima and minima. Persistent homology-based observables such as the persistent entropy may be a useful tool for detecting novel phases and parameter regimes of complex lattice models.

We focused on the case of a one-dimensional lattice for ease of understanding. However, persistent homology can also be readily applied to higher dimensional intensity profiles, where it can not only characterize the number of high intensity clusters, but also detect higher-dimensional topological features such as holes and loops. This could be useful for distinguishing localized, extended, and compact localized modes in disordered flat band lattices such as the Lieb lattice and its higher dimensional generalizations Lieb2D ; Lieb3D .

Another interesting direction would be to investigate whether persistent homology-related observables are amenable to analytical or semi-analytical calculations. For example, in 1D tight binding lattices the energy-dependent localization length can be computed without diagonalizing the Hamiltonian by evaluating a recurrence relation satisfied by the eigenstates loc_review ; Anderson_review . Analytical results may reveal to which specific features of eigenstates persistent homology is most sensitive.

Since persistent homology captures multiscale properties of wavefunctions, the techniques developed here may be useful for characterizing the response of active photonic systems, such as supercontinuum and frequency comb generation in microring resonators, where the intensity can be measured over a huge range of scales comb . Thus, we envision the techniques could be further developed for advanced studies, which may find applications in deep-learning-based photonic design and AI-empowered nanophotonics that will continue to grow in the coming decade final1 ; final2 .

This research was supported by the National Natural Science Foundation of China (12134006), National Research Foundation, Prime Minister’s Office, Singapore, the Ministry of Education, Singapore under the Research Centres of Excellence Programme, China Postdoctoral Science Foundation (BX2021134, 2021M701790) and the Polisimulator project co-financed by Greece and the EU Regional Development Fund.