Localization and symmetry breaking in the quantum quasiperiodic Ising glass

Quasiperiodic modulation can prevent isolated quantum systems from equilibrating by localizing their degrees of freedom. In this article, we show that such systems can exhibit dynamically stable long-range orders forbidden in equilibrium. Specifically, we show that the interplay of symmetry breaking and localization in the quasiperiodic quantum Ising chain produces a \emph{quasiperiodic Ising glass} stable at all energy densities. The glass order parameter vanishes with an essential singularity at the melting transition with no signatures in the equilibrium properties. The zero temperature phase diagram is also surprisingly rich, consisting of paramagnetic, ferromagnetic and quasiperiodically alternating ground state phases with extended, localized and critically delocalized low energy excitations. The system exhibits an unusual quantum Ising transition whose properties are intermediate between those of the clean and infinite randomness Ising transitions. Many of these results follow from a geometric generalization of the Aubry-Andr\'e duality which we develop. The quasiperiodic Ising glass may be realized in near term quantum optical experiments.


I. INTRODUCTION
Nearly 60 years ago, Anderson discovered that quenched disorder could localize quantum particles and thus prevent the transport necessary for equilibration in isolated systems [1]. The recent interest in the role of interactions  and rapid experimental developments in synthetic quantum systems [30][31][32][33][34][35][36][37] have led to a deeper understanding of the full range of consequences of Anderson's original observation. The phenomenology of the localized phase is now better understood as a form of integrability with local conserved quantities [14,[38][39][40][41][42][43]; the dynamics of entanglement has emerged as a unifying framework for understanding thermalization [10][11][12]; and the long-lived coherence of localized systems may serve as a resource for quantum information processing [44][45][46].
A particularly intriguing proposal is that localization can dynamically protect long-range order in highly excited states even when such orders are forbidden in equilibrium [47]. The central idea may be illustrated in the 1D ferromagnetic Ising chain. Ferromagnetic order in the ground state is usually destroyed in excited states because of the proliferation of domain walls (an argument that goes back to Peierls). However, if quenched disorder can localize the domain walls, then the system never reaches equilibrium, and any symmetry-breaking pattern imprinted in the spin state at t ¼ 0 can persist for all time. This Ising glass order clearly exists in the transverse-field Ising chain in strong disorder treatments [8,47] and has been observed numerically in small interacting chains [19]. Localization protection has also been argued to extend to a host of more exotic orders [46][47][48][49] and to periodically driven (Floquet) systems [50][51][52][53][54].
Localization, however, does not require disorder, as was first recognized by Azbel [55], and Aubry and André (AA) [56] in the single-particle context. These authors discovered that sufficiently strong quasiperiodic potentials can localize a quantum particle. References [18,24,57] extended these results to the interacting many-body case and argued that many-body localization can persist even at high-energy density. Quasiperiodic potentials arise naturally in optical experiments using lasers with incommensurate wave vectors. Accordingly, many experiments in such systems have now observed single-particle localization [58][59][60][61][62][63] and, more recently, have also pushed into the interacting regime and high-excitation energy densities to provide evidence for the many-body localized phase [30,32,37].
As quasiperiodic systems can show both localized and delocalized behavior already in the 1D noninteracting context, they offer a well-controlled platform to study the interplay of localization and symmetry breaking. In this article, we study the effects of quasiperiodic modulation on the canonical quantum Ising chain. The most salient dynamical feature is a stable quasiperiodic Ising glass in which all excited states exhibit Ising symmetry-breaking order (red in Fig. 1). This excited-state order melts if either the ground state becomes paramagnetic or the domain-wall excitations delocalize; we find both types of transition. Remarkably, the excited-state Ising glass order parameter exhibits an essential singularity at the transition, with no signatures in the ground-state ordering.
In quench experiments, the quasiperiodic Ising glass phase appears in the persistence of arbitrary initial longitudinal magnetization (i.e., in the direction flipped by the Ising symmetry) after a short transient. This glass is accessible in current experiments in quantum optical Ising spin simulators, such as those that have been implemented in ion traps [34,64] and Rydberg systems [65,66]. Experimentally, it is better to modulate the effective spin-spin interaction (as opposed to modulating the field) by quasiperiodically modulating the positions of the trapped spins, as this is the regime most favorable to finding the Ising glass. We have accordingly focused the detailed study in this manuscript on the coupling, rather than field, modulated regime. Our rigorous analytic control extends only to nearest-neighbor spin-spin interactions, where the system can be fermionized, but we expect the Ising glass to persist in the presence of weak additional interactions, just as the quasiperiodically modulated manybody localized phase of particles persists in Ref. [18]. We further discuss both potential experimental realizations and the stability to interactions in Sec. VII.
From an equilibrium condensed-matter perspective, the zero-temperature phase diagram is interesting in its own right. There are paramagnetic (PM), ferromagnetic (FM), and quasiperiodically alternating ferromagnetic (QPFM) orders in the ground state. Moreover, the low-energy excitations exhibit extended, localized, and critically delocalized behavior depending on the strength of the quasiperiodic modulation. This leads to an array of possible combinations, which we have summarized in Fig. 1. In the literature on single-particle localization, "critical" wave functions exhibit multifractal scaling behavior; throughout the manuscript, we refer to these as "critically delocalized" in order to distinguish them from the "critical" properties associated with quantum phase transitions [67].
The zero-temperature quantum phase transitions in Fig. 1 lie in two distinct universality classes. We find that weak quasiperiodic modulation is irrelevant at the clean Ising transition, so the parabolic phase boundary in Fig. 1 exhibits quantum critical scaling with dynamic exponent z ¼ 1 and extended low-energy excitations. At strong modulation, we find a new quantum Ising transition separating the QPFM from the paramagnet. This transition exhibits dynamical critical behavior that is intermediate between that of the clean Ising transition and the infinite randomness transition that arises in the disordered model. In particular, while the correlation length diverges with ν ¼ 1, as in the clean transition, the low-energy excitations undergo a transition from critically delocalized to localized, coincident with the symmetry breaking, with an apparent exponent z ¼ 2.
Our results make use of a variety of analytic and numerical techniques. We would like to especially flag a new relative of the celebrated Aubry-André duality which we have discovered. We dub this transformation "AAA triality" as it maps cyclically among three related models. It turns out that the "self-trial" point in the quasiperiodic Ising model sits on the phase boundary between the paramagnet and QPFM, giving us analytic access to the unusual quantum critical properties. Our triality arguments explain the energy-independent wave-function criticality on the "pure modulation," the J=h ¼ 0 axis in Fig. 1, which has been observed numerically before in Refs. [68,69].
The noninteracting quasiperiodic models of Azbel, Aubry-André, and their generalizations have been extensively studied by mathematicians and physicists over the last 30 years for a variety of reasons [70][71][72][73][74][75][76][77][78][79][80][81][82]. These 1D models exhibit a single-particle localization-delocalization transition at finite modulation which mimics the metal-insulator transition in 3D disordered systems. This is in striking contrast to the disordered Anderson model in 1D, which is localized for any disorder strength [83]. Mathematically, the models have a surprisingly rich analytic structure, exhibiting dualities, critically delocalized phases with fractal spectra and connections to higher-dimensional Hofstadtertype models [56,84]. More generally, they offer a window into quantum localization without the epiphenomena associated with rare region effects in disordered systems.
There has been significant previous work on aperiodic and/or quasicrystalline quantum Ising chains [85][86][87]. These models have Ising couplings chosen from a finite set according to a recursive substitution rule, or by quasicrystalline projection. There have also been several previous studies of the zero-temperature properties in certain regions of the phase diagram of the incommensurately modulated Ising chain [68] or, equivalently, in the modulated p-wave superconductor [69,88,89]. These works were largely numerical, and many of the features of the zero-temperature phase diagram were missed. Our analysis in various limits, and particularly the AAA triality, unify and extend their results. Moreover, the quasiperiodic Ising glass order in the excited states was completely overlooked. The organization of the paper is as follows. We begin in Sec. II with a precise definition of the Ising model, a review of its fermionization, and various salient facts about quasiperiodic modulation in chains. While much of Sec. II is review, the geometric interpretation of the AA duality in two dimensions may provide an alternative perspective for many readers. In Sec. III, we derive the ground-state symmetry-breaking phase diagram. We turn to the localization properties of the low-energy excitations in Sec. IV. With these basic properties in hand, we discuss the zero-temperature quantum critical behavior in Sec. V. We investigate the properties of the excited-state Ising glass order in Sec. VI and its melting transition in Sec. VI A. We conclude with a discussion of the role of interactions, possible experimental realizations, and other open questions.

II. GENERAL PROPERTIES OF THE MODEL
The Hamiltonian of the one-dimensional quasiperiodic transverse-field Ising model (TFIM) is The model is illustrated in Fig. 2(a). Here, σ α j are the Pauli matrices, with j ∈ Z running over the sites in the chain and α ¼ x, y, z; Q is the wave vector of the modulation in units where the lattice spacing is a ¼ 1; and the phases ϕ and ϕ þ Δ shift the positions of the maxima of the couplings relative to the underlying lattice. The wave vector Q is commensurate with the underlying lattice if Q=2π ¼ p=q is rational, and it is incommensurate otherwise. We choose the wave vector of the modulation of the Ising coupling and the transverse field to be the same for simplicity and since this is natural if the modulation arises from the same physical source (e.g., an incommensurate laser potential).
Global symmetries.-The quasiperiodic TFIM has several global symmetries. The eponymous Ising symmetry is given by G ¼ Q i σ z i -this is the symmetry that breaks spontaneously in the T ¼ 0 ferromagnetic and localized Ising glass phases. The Hamiltonian H is also symmetric under complex conjugation K, which is antiunitary. Finally, for special values of the modulation phases ϕ and Δ, the model is symmetric under reflections across site k, j → k − j or bond k þ 1=2, j → k þ 1=2 − j. For example, at Ising duality.-Under the duality transformation, Here, H maps onto another incommensurate TFIM H 0 with the role of the field and bond couplings interchanged (and different boundary conditions). Formally, the duality maps The duality swaps paramagnetic and ferromagnetic phases but leaves the dynamical nature of the bulk single-particle excitations (spectrum and wave-function localization) invariant. This can be seen most easily from the fermionization (see below) of H and its dual H 0 , whose noninteracting Hamiltonians agree precisely up to a translation by half a unit cell so that their single-fermion modes are identical up to this half-translation.
Fermionization.-The TFIM has a well-known fermionic representation, which we review and extend to the quasiperiodic case here. The Jordan-Wigner transformation introduces a pair of Majorana fermion operators for each spin-1=2: The γ operators are Hermitian and satisfy the canonical anticommutation relation fγ i ; γ j g ¼ 2δ ij . The transformation maps the TFIM to a quadratic Majorana chain with Hamiltonian [see Fig. 2 The dynamical and symmetry-breaking properties of the TFIM Hamiltonian H follow from the properties of the single-particle Hamiltonian H defined by Hermiticity of H requires that H ij is an imaginary antisymmetric 2L × 2L matrix, where L is the number of spins. Thus, the eigenvalues of H come in AEe pairs corresponding to complex conjugate eigenmodes, ψ andψ. Labeling the L positive-energy eigenmodes by the index α, we can diagonalize H into the familiar form where The complex fermion operators c, c † satisfy the usual anticommutation relations In the ground state, the paramagnetic phase corresponds to the topologically trivial phase of the Majorana chain, while the ferromagnetic phase maps to the topologically nontrivial phase. The simplest way to detect the topological phase of the Majorana chain is with open boundary conditions, in which case the topologically nontrivial phase possesses a pair of zero-energy Majorana modes localized at the boundaries of the chain [90]. The many-body groundstate space is accordingly doubly degenerate, as the fermionic mode defined by the two zero-energy Majorana operators can be occupied or unoccupied at zero cost. We use this approach to extract the ground-state phase diagram of the quasiperiodic model in Sec. III.
The symmetries of Eq. (1) appear in the fermionic language as follows. The global Ising symmetry operator G ¼ Q j σ z j maps to the fermionic parity operator G ¼ Q j ð−iγ 2j γ 2jþ1 Þ, while the symmetry under complex conjugation forces H to be bipartite. The action of the Ising duality on the Majorana chain shifts all the site labels by a half: j → j − 1=2 as mentioned above.
As all eigenstates of H correspond to Slater determinant states of the fermions γ, they all satisfy Wick's theorem. This allows evaluation of the spin-spin correlation function, as a Pfaffian of the fermionic Green function restricted to the diagonal block from 2i þ 1 to 2j. While this representation is not easy to use analytically, it allows straightforward numerical computations of the exact correlation functions in large systems (e.g., up to L ¼ 1000 in this work). Evaluating these correlators at large separation ji − jj allows us to numerically extract the magnetization as where M ψ i is the magnetization of spin i in state jψi. The 2D model.-The incommensurate TFIM has a second useful representation in terms of a 2D model of noninteracting complex fermions. To derive this, consider first a complex fermionic model with the same singleparticle Hamiltonian as Eq. (8) [91]: Above, d i destroys a fermion at site i, and d i , d † j satisfy the usual complex anticommutation relations. On varying the phase ϕ between ½0; 2πÞ, we generate a family of distinct 1D Hamiltonians. Treating ϕ as the momentum along an extra dimension y and inverting the Fourier transformation, we obtain a 2D tight-binding model: whose spectrum at fixed y momentum ϕ reproduces the spectrum of the 1D modelH 1D ðϕÞ; see Fig. 2 (14) is off-diagonal. The localization properties of the excitations of the 1D TFIM map onto those of theH 2D in the Landau gauge at fixed y momentum ϕ. By construction, if the eigenstates of the 2D model are delocalized in the x direction, then the eigenstates of the corresponding 1D model are extended, while if the states of the 2D model are localized in the x direction, then the 1D model is localized. As we will see, the 2D picture is a surprisingly useful geometric aid for identifying localized, extended, and critical phases of the excitations.
Aubry-André model.-We need several properties of the original Aubry-André hopping chain in the analysis of the quasiperiodic TFIM. The AA model has the following Hamiltonian [56], For V > t, the single-particle states are localized at all energies, while for V < t, they are extended at all energies. At the critical point, V ¼ t, the states exhibit multifractal properties. Furthermore, the localization length diverges at the critical point with exponent ν ¼ 1. More precisely, Many features of the phase diagram follow Aubry-André duality. This duality corresponds to a π=2 rotational symmetry of the associated 2D model. The 2D model is that of a particle hopping on an anisotropic square lattice with flux Q per plaquette and hopping strength t and V in the x and y directions, respectively; at t ¼ V, this is the Hofstadter model, whose fractal character is well known [84]. The π=2 rotation swaps t and V and accordingly swaps the localized and extended phases of the original 1D model; clearly, t ¼ V is self-dual. Moreover, the rotation ensures that the localization properties are energy independent for any t and V, as we will see in more detail in Sec. IV from an analysis of the characteristic polynomial.
Lack of Aubry-André duality.-The incommensurate TFIM clearly lacks the Aubry-André rotational symmetry, as Fig. 2(c) does not map onto itself (up to swapping couplings) under rotation by π=2. While one can embed the 2D model into a larger class of 2D models that have the requisite vertical bonds, these would, in general, need to be staggered in the y direction; thus, they do not correspond to 1D incommensurate chains. The π=2 rotation in the larger model space therefore does not define a duality map on the incommensurate TFIM. This is both a blessing and a curse: The incommensurate TFIM exhibits rich phenomena not allowed by AA duality, such as energy and Q-dependent mobility edges, but it is correspondingly harder to rigorously analyze. In Sec. IV, we see that various limits of the quasiperiodic model have higher rotational symmetry when viewed in 2D. These will give us analytic control of the phase diagram in those limits.
Parameter regime we study.-Except for special lines in coupling space, the ground-state phase diagram has not been studied. The full dynamical phase diagram lives in an eight-dimensional space, parametrized by h, J, A h , A J , Q, ϕ h , Δ, and the excitation energy e; after normalizing the units of energy, there are still seven dimensionless parameters. Where simple enough, we will provide general expressions that apply to the entire phase diagram. However, as seven-dimensional phase diagrams are unwieldy, we mostly focus on the manifold defined by A h ¼ 0, which is already rather interesting. In units where h ¼ 1, the relevant parameters controlling the phase diagram are then J, A J , Q, ϕ, and the excitation energy e (note that Δ has no effect when A h ¼ 0). Most of the analytic features we derive hold for any incommensurate Q and ϕ, but we focus our numerical results (notably for spectra) on Q=2π ¼ ð ffiffi ffi 5 p þ 1Þ=2, the golden mean. By Ising duality, our analysis also produces the phase diagram at A J ¼ 0.

III. GROUND-STATE SYMMETRY BREAKING
The quasiperiodic Ising model has three ground-state phases: a paramagnetic phase, a ferromagnetic phase, and a quasiperiodically alternating ferromagnetic phase. The two latter phases spontaneously break the Ising symmetry so that the spin-spin correlation function hσ x i σ x j i has longrange order. In the simple FM, all of the spins magnetize in the same direction (although the magnetization need not be spatially uniform). In the QPFM, on the other hand, neighboring spins align or antialign quasiperiodically because of the presence of antiferromagnetic links where J j < 0. This magnetization pattern is analogous to that of a ground-state Ising spin glass as it has domain walls built in. However, the ground-state order is unfrustrated because the Hamiltonian is gauge equivalent to a purely ferromagnetic model with J j > 0 [92]. We therefore refer to both Ising ordered phases as "ferromagnetic".
The clean model (A J ¼ A h ¼ 0) spontaneously breaks the Ising symmetry for jJj > jhj. On general grounds, the corresponding gapped ferromagnetic and paramagnetic phases should persist in the presence of small incommensurate modulation A h , A J . To find the phase boundaries in general, we recall that the Ising symmetry-breaking phase corresponds to the topological phase of the fermionic representation, Eq. (7), which famously hosts zero-energy Majorana modes, Γ L and Γ R , bound to the left and right ends of an open chain [90]. To detect ground-state symmetry breaking, we look for normalizable boundary modes.
For simplicity, we focus on the left edge of a semiinfinite chain. As the Hamiltonian Eq. (7) is bipartite, only connecting the even and odd sublattices, the zero mode localized at the left edge must have the following form: (The right mode would be localized on the odd sublattice.) Substituting in the eigenvalue equation ½H; Taking a logarithm, we obtain In order for the zero mode to be normalizable, S j must decrease sufficiently rapidly with j at large distances.
In the absence of modulation, S ¼ j log jhj=jJj clearly grows linearly with j; the sign of log jhj=jJj thus determines the normalizability of Γ L . Modulation at irrational wave vector Q causes S j to fluctuate with j. However, on distances long compared to Q −1 , we expect S j to still have a linear trend as the sum averages over the modulation. This trend can be extracted formally by averaging over the phase, The sign of I determines if the zero modes are normalizable (I < 0) or not (I > 0). The Ising phase boundary is at I ¼ 0, which is shown by the black curve on the phase diagram at A h ¼ 0 in Fig. 1. We evaluate I in Appendix A. The contour lies at In Appendix G, we present an alternative derivation of the phase boundaries by approaching the incommensurate Q limit through a series of longer and longer commensurate wave numbers 2πp=q → Q. This provides a more rigorous treatment of the role of incommensuration. The two derivations are in perfect agreement.
Numerical evaluation of long-range order in the ground state agrees with the phase boundary determined from the above analysis, Eq. (23). In Fig. 3, the three panels present the average magnitude of the magnetization squared on representative cuts through the phase diagram, calculated at for open chains of length L ¼ 1000. More precisely, we detect the ground-state order through the correlator, which is nonzero in both the FM (A J < J) and QPFM (A J > J) phases. Here, O gs is an analog of the more familiar Edwards-Anderson order parameter defined using the square of the spin-spin correlation function rather than the absolute value. The square brackets ½· j indicate averaging over a small spatial window (j ¼ −w; …; w) to smooth out the quasiperiodic modulation of the magnitude. In the thermodynamic limit L → ∞, we expect where ½jMj is the average of the absolute value of the site magnetization and w is the width of the spatial averaging window.
We end with a few comments. First, the general derivation reproduces the critical point of the clean A J ¼ A h ¼ 0 model. Second, Ising duality immediately implies that the phase diagram in the A h =J, h=J plane at A J ¼ 0 looks identical to that in Fig. 1 after swapping the ferromagnetic and paramagnetic phases. Third, whether the ground state breaks Ising symmetry or not is independent of Q, so long as Q is incommensurate, as only the average log jhj=jJj over a period matters to the longdistance behavior. Fourth, the fluctuations of S relative to its linear trend are Oð1Þ at distance j. This is significantly less than the Oð ffiffi j p Þ fluctuations that develop in the case of usual disorder with independent random couplings. Finally, the cusp at A J ¼ J ¼ 2 is real and indicates that the transition on the parabolic boundary between the FM A. CHANDRAN and C. R. LAUMANN PHYS. REV. X 7, 031061 (2017) 031061-6 and PM and on the vertical boundary between the QPFM and PM have different character. We return to this in greater detail in Sec. V.

IV. LOCALIZATION OF EXCITATIONS
The dynamical phase diagram of the quasiperiodic TFIM follows from the properties of the fermionic excitations described by Eq. (7). The main feature of interest for the stability of excited-state order is the spatial extent of the single-particle wave functions-in other words, whether they are extended, critically delocalized, or localized. In general, at any given point in the phase diagram, these properties are energy dependent and the system may exhibit mobility edges. They can be studied numerically quite effectively. However, there are many special lines in the coupling space with enhanced symmetry, which provides analytic control and energy independence. On the A h ¼ 0 plane represented in Fig. 1, both the axes, A J ¼ 0 and J ¼ 0, and the large coupling limit, J → ∞, have such enhanced symmetry. These limits will be sufficient to asymptotically characterize all of the features in Fig. 1. We study these limits in the subsections below before turning to numerics to support the bulk of the phase diagram.
The simplest diagnostic of localization is the inverse participation ratio defined for a given eigenstate α as In finite-size studies, the scaling of the IPR in a given energy window, IPR ∼ 1=L γ , detects the dynamical phase.
Formally, these phases correspond to spectra that are absolutely continuous, pure point and singular continuous (fractal), respectively. We use both diagnostics in the following analysis.
A. Clean limit A J → 0 In the absence of the incommensurate modulation, the model reduces to the usual nearest-neighbor Ising model. It has extended excitations for all parameter values (J=h) and at all energies.
The usual Ising critical point at J=h ¼ 1 has gapless extended excitations at all energies. As we argue in more detail in Sec. V, the parabolic ground-state phase boundary extending from the clean critical point [J > A J , Eq. (23)] lies in the same universality class; thus, we expect the lowenergy excitations to remain extended all along this boundary. However, mobility edges are allowed at higher energy. These features are visible in the numerical data shown in Figs. 5(a) and 5(b).
In this regime, the ground state is very close to the ideal ferromagnet with all the spins pointing in the þx or −x direction. The excitations are domain walls, The Hamiltonian for a domain wall is This is simply the AA model, Eq.  24) and (35), respectively. The spatial averaging window is of size w ¼ 5, and we sample 50 eigenstates in the excited-state average. In the first panel, both the ground-and infinite-temperature order turn on at the same coupling as the excitations localize across the transition. In the second panel, the ground-state order is unchanged even as the excited-state order turns on because of the localization of the domain walls. In the third panel, the ground-state order turns on, but the excitations are extended so that the excited-state order remains zero.
of the extended to localized phase boundary at large J in Fig. 1. One can also see the emergence of the AA model in this limit by considering the associated 2D model. In the large J ≫ A J , h limit, one pairs the fermionic sites across the horizontal J links and obtains an effective square lattice with horizontal links of strength h and vertical links of strength A J .
This regime is in the bottom left corner in Fig. 1. The ground state is paramagnetic and is very close to the product state j↑↑↑ Á Á Á ↑i. The excitations are spin flips. Analogously to Sec. IV B, the effective Hamiltonian for the spin flips is given by up to a constant. This off-diagonal Aubry-André model admits the AA duality and accordingly has energyindependent localization properties. The model has been previously studied by Thouless and co-workers [71,93], who found two dynamical phases: extended for A J < J and critically delocalized for A J > J. The same authors also computed the multifractal properties at the transition A J ¼ J. Asymptotically, this coincides with the diagonal dashed line near the origin of Fig. 1.
As an aside, the critically delocalized phase for A J > J is a special feature of the A h ¼ 0 plane. At small A h and Δ ≠ 0, the critically delocalized phase is destroyed and the states for A J > J are fully localized [94].

D. Large A J ≫ J, h
The 2D model in the A J → ∞ limit reduces to a collection of decoupled vertical zigzag wires (dashed orange in Fig. 2), two per original spin. Without h, it is clearly localized for the 1D model, as it decouples at every bond.
To determine the stability of the localization to small J and h, note that the pair of wires at position j have dispersion where ϕ is the y momentum [95]. The energy dispersions are shifted by y momentum Q between adjacent wire pairs. Thus, for incommensurate Q ∼ Oð1Þ, neighboring wire states at fixed momentum ϕ are typically nondegenerate, with an energy splitting of order A J . To leading order, this implies that h is typically an off-resonant perturbation, and localization in the x direction persists. This argument can be generalized to higher order using Diophantine properties of Q, and it closely mirrors arguments that can be made in the anisotropic Hofstadter problem.
Thus, we expect that all states are localized at large A J ≫ J, h to the far right of Fig. 1. We note that the 2D model does not have enhanced symmetry in this limit; thus, the analysis does not provide as complete control as it does in the Hofstadter model. E. Pure modulation J = 0 and AAA triality The analysis on the J ¼ 0 axis of Fig. 1 is significantly more involved than the previous limits, as it turns out that the model exhibits a hidden symmetry that allows us to define a "AAA triality" operation in analogy with AA duality. The existence of this triality implies that the localization properties of the wave functions are independent of energy, as in the AA model. It also enables an exact evaluation of certain properties of the secular equation, from which we find that the spectrum is critically delocalized for A J =h < 2 and localized for A J =h > 2. The selftrial point lies at A J =h ¼ 2. We note that this explains various numerical observations made on this axis in previous studies of the p-wave superconducting chain [68,69,88,89].
The AAA triality is best understood as a geometric transformation of the associated 2D model. At J ¼ 0, this model decouples into two interpenetrating honeycomb lattice "layers," each pierced by flux 2Q per hexagonal plaquette (see Fig. 4). At A J =h ¼ 2, since the layers are decoupled, the system is symmetric under rotation by 2π=3 independently in each layer, about any site. These are physical rotations in the 2D model, which need to include gauge transformations forH 2D to return to Landau gauge.
In the usual AA model, π=2 rotation in the associated 2D model leads to another model with couplings V and t swapped. Rotating both layers by 2π=3 in the honeycomb model rotates the couplings associated with the three bond angles into one another. This suggests the utility of generalizing the 2D model to have three independent couplings, A 1 , A 2 , and A 3 , corresponding to the three different bonds. In Landau gauge, this corresponds to a 1D incommensurate model with Hamiltoniañ The original 1D TFIM at J ¼ 0 corresponds to A 3 ¼ h, In the extended AAA family of models, the 2π=3 rotation mapsH AAA toH 0 AAA with cyclically permuted couplings. This is the AAA triality.
With this triality in hand, it is possible to analyze the wave-function properties of the AAA model in considerable detail. We relegate this analysis to the appendixes, as it entails a considerable calculational detour. The upshot is that the 1D incommensurate TFIM has a critically delocalizedto-localized transition, at all excitation energies, at the selftrial point, A J ¼ 2h. The multifractal properties of this transition are a subject of ongoing work [96].

F. Near the clean Ising transition
As discussed in more detail in Sec. V, the ground-state transition between the PM and the FM is described by the clean Ising critical theory for A J =h < 2. Thus, we expect that low-energy excitations in the vicinity of this transition are extended. This constrains the boundary of the fully localized region (red in Fig. 1) to stay strictly above the parabolic ground-state transition curve. As this transition terminates at the cusp ðA J =h; J=hÞ ¼ ð2; 2Þ, the low-energy states need not be extended beyond this point.
The red-blue boundary in Fig. 1 is the simplest way to satisfy these constraints and connect to the large J=h limit of Sec. IV B. However, the precise location depends on Q and can only be determined numerically.

G. Numerical support away from the limits
Away from the special lines discussed in the previous subsections, the localization properties of the wave functions are both Q and energy dependent. We rely on numerics to confirm the features summarized in Fig. 1. We have extensively investigated Q=2π ¼ ð ffiffi ffi 5 p þ 1Þ=2 and checked the qualitative features for several other incommensurate wave vectors. The general features are as follows: (1) The extended states on the clean (A J ¼ 0) axis persist in the presence of small modulation A J . With increasing A J , the high-energy states localize before the lower-energy states. (2) Near the Ising transition along the phase boundary above the diagonal A J ¼ J, the gap closes and reopens linearly, and all excitations are extended up to a finite energy above the gap. See the first two columns of Fig. 5 for representative spectra and IPR behavior and Sec. V for more discussion. (3) For A J =h > 2, all states are localized, consistent with the analytically proven behavior at J ¼ 0 (Sec. IV E) and the analysis at large A J (Sec. IV D) and large J (Sec. IV B). At large J=h, we observe the expected energy-independent localization transition near A J =h ¼ 1. (4) Numerics confirm the qualitative trend of the full localization boundary in the A J =h < 2 regime: that the boundary approaches a vertical asymptote at A J =h ¼ 1 for large J=h and that it remains above the Ising ground-state transition, ending at the point ðA J =h; J=hÞ ¼ ð2; 2Þ. (5) On the J ¼ 0 line at A J =h < 2, we confirm that all states are critically delocalized by calculating the scaling of the IPR. At high energy, the states localize on the introduction of J > 0. However, we find that the lowest energy states (above the gap) continue to exhibit critical IPR scaling throughout the triangle below the diagonal A J ¼ J (see representative data in Fig. 5). This behavior defines the purple region in Fig. 1.

V. GROUND-STATE QUANTUM PHASE TRANSITIONS
In this section, we focus on the properties of the zerotemperature quantum phase transition between the groundstate paramagnetic and ferromagnetic phases. The phase transitions above and below the diagonal A J ¼ J in Fig. 1 are qualitatively distinct. Above the diagonal (A J < J), the quasiperiodic modulation is irrelevant, and the transition lies in the standard 1D quantum Ising universality class. Below the diagonal (A J > J), the transition lies in a new "quasiperiodic Ising" universality class, with behavior intermediate between the clean Ising critical point and the infinite randomness critical point (IRCP), which governs the disordered Ising transition [97].
At both transitions, the correlation length ξ diverges with exponent ν ¼ 1. From the relation ξ ∼ −1=I, where I is the integral governing the convergence of the Majorana boundary mode, Eq. (22), it is straightforward to show that where δ is the deviation from the phase boundary at I ¼ 0. This is consistent with the Harris-Luck criterion [87,98], which imposes that ν ≥ 1=d ¼ 1 for phase transitions in the presence of incommensurate modulation.
FIG. 4. The 2D model associated with the incommensurate TFIM at A h ¼ 0, J ¼ 0 decouples into two interpenetrating layers (blue and orange) of a honeycomb lattice. Each hexagonal plaquette is pierced by flux 2Q. We have adjusted the geometric embedding of the lattice points relative to Fig. 2 in order to emphasize the symmetry associated with rotation by 2π=3. All parallel bonds carry the same coupling magnitude (though the hopping phase depends on the choice of gauge). The incommensurate TFIM at J ¼ 0 is embedded into a larger AAA family of models with couplings

A. Clean Ising transition
While ν ¼ 1 at both transitions, the dynamical properties are quite distinct. Above the diagonal, the phase boundary connects to the standard clean Ising transition at J=h ¼ 1, A J =h ¼ 0. At small A J , as the Harris-Luck criterion is marginal, it is natural to conjecture that weak quasiperiodic modulation is (marginally) irrelevant. This would imply that, in the vicinity of the phase boundary, (1) the dynamical critical exponent z ¼ 1, so the gap Δ closes linearly with δ, (2) the low-energy excitations are extended, as these are the excitations that mediate the phase transition. Both of these expectations are borne out numerically quite beautifully.
The top left panel of Fig. 5 shows the excitation spectrum versus J=h for a vertical cut at A J =h ¼ 0.5. The vertical cut intersects three regimes: the critically delocalized PM (0 ≤ J=h ≤ 0.5), the extended PM (0.5 ≤ J=h ≤ J c =h ¼ 1.0625), and the extended FM (J c =h ≤ J=h). The boundary zero mode associated with the FM is clearly visible to the right of the transition. As promised, the gap closes linearly, and the low-energy excitations above the gap remain extended. The extension of the wave functions are indicated qualitatively by the coloring of the states; quantitatively, the lower panel shows the scaling exponent γ of the low-energy IPR with system size, The lower-left panel shows that γ is near 1 for the lowenergy excitations for J=h > 0.5 and, in particular, near the transition. The center column of Fig. 5 presents similar data at stronger incommensurate modulation, A J =h ¼ 1.5. Note that the axes are zoomed in on the phase boundary and on low energies relative to the previous plot. The high-energy excitations (e=h > ≈0.2) are localized across the transition. Nevertheless, the gap closes linearly (center top), and the low-energy excitations remain extended nearby (γ ≈ 1, in the bottom panel). Thus, the ground-state symmetry-breaking transition is still in the clean quantum Ising universality class. The color value is given by − log IPR= log L, which varies between 0 (localized, black) and 1 (extended, copper). (Lower row) The IPR exponent α is defined by the scaling ½log IPR ∼ −α log L, where the mean ½· is taken over the lowest twentieth of the excitations. Here, α varies between 0 (localized) and 1 (extended). To indicate the finite-size trends approaching the thermodynamic limit, we fit data from both L ¼ 100, 200, 350, 500, 750, 1000 (blue) and L ¼ 350, 500, 750, 1000 (green).

B. Quasiperiodic Ising transition
Below the diagonal, the dynamics near the transition between the QPFM and PM change character rather dramatically. Along the J ¼ 0 line, the symmetry-breaking transition at A J =h ¼ 2 coincides with the transition from critically delocalized to localized excitations at all energies, as shown in Sec. IV E. There are no fully extended excitations. In analogy with the irrelevance of A J at the clean transition, we conjecture that J is irrelevant at low energies to the strong quasiperiodically modulated transition along the vertical phase boundary.
Numerically, this conjecture is borne out by the following observations (see third column of Fig. 5 for representative data along a particular cut at J=h ¼ 0.5): (i) The single-particle gap closes on approaching the transition from the paramagnetic side with an exponent z ≈ 2. However, the gap does not reopen on the symmetry-breaking side of the transition. (ii) All excitations are localized in the symmetrybreaking phase at A J =h > 2. (iii) The low-energy excitations (bottom tenth of states above the gap) in the paramagnetic phase exhibit critical IPR scaling and extreme finite-size sensitivity. See the lower row of Fig. 5. The IPR exponent γ lies properly between 0 and 1 and exhibits strong finite-size fluctuations. (iv) At higher energy (and J=h > 0), the excitations can be localized even in the paramagnetic phase. We do not know whether the critically delocalized lowenergy states are separated from these localized states by a mobility edge or whether there is a long crossover. While we leave a full analytic study of this transition to forthcoming work [96], it is clear that the quasiperiodic Ising transition is intermediate between the well-known clean and infinite randomness critical points. In the former, the critical excitations are extended on both sides of the transition with dynamical exponent z ¼ 1; in the latter, the excitations are localized on both sides of the transition, and the dynamics are activated (roughly, z → ∞). Across the quasiperiodic transition, the excitations pass from critically delocalized and gapped (with exponent z ¼ 2) to localized and gapless across the transition. Also, unlike the infinite randomness transition, the correlation functions in the quasiperiodic case do not acquire broad distributions at long distances. For example, the mean and typical decays of the boundary mode are governed by the same correlation length ξ with exponent ν ¼ 1, while in the infinite randomness transition, they are governed by two distinct diverging length scales (with ν ¼ 2 and 1, respectively).

VI. EXCITED-STATE ISING GLASS ORDER
We expect localized Ising glass order at all energy densities in the regions where (1) the ground state breaks Ising symmetry, and (2) all the excitations are localized. This region is indicated in red in Fig. 1. Note that the boundary of the region above the diagonal (J > A J ) depends on Q. Although we do not have an explicit functional form for this Q dependence, the analysis in Sec. IV constrains the boundary to be at A J =h ¼ 1 as J=h → ∞ independent of Q and to lie strictly above the parabolic ground-state boundary.
Representative correlation functions.- Figure 6 shows the spin-spin correlation function hσ x L=2 σ x L=2þd i in the ground state (blue) and a random excited state drawn from the infinite-temperature ensemble (green) at three different representative points in the phase diagram. In the left panel, the spin-spin correlation function decays rapidly in both states, confirming that the point A J =h ¼ J=h ¼ 0.5 is in the PM phase. The center panel shows the correlation functions at a point in the phase diagram where the ground state is ordered and the excitations above it are extended. The longrange order in the ground state is clearly detected by the blue curve, which approaches a nonzero value at long distance d. The green curve, on the other hand, decays quickly to zero and shows that infinite-temperature eigenstates are not ordered. This result agrees with the Mermin-Wagner-Peierls theorem that states there can be no long-range order at any finite temperature in 1D.
The right panel provides evidence for Ising glass order in a randomly chosen infinite-temperature state at a point in the phase diagram where the ground state is ordered and the single-particle excitations are localized. The green curve does not decay as d → ∞; instead, it fluctuates on an order-one scale depending on whether spin L=2 and spin L=2 þ d are aligned or antialigned in the chosen infinitetemperature state.
Frozen order parameter.-The quasiperiodic Ising glass order is detected by an order parameter that generalizes the ground-state order parameter, Eq. (24), to finite energy density states: Here, ½· E;j indicates averaging over both eigenstates jEi in some energy window and over a small spatial window j ∈ ½−w; w to suppress the quasiperiodic fluctuations. Note that O exc is nonzero, as L → ∞ only in states with longrange Ising symmetry breaking, such as the Ising glass. We plot the order parameter along several representative cuts in the phase diagram in Fig. 3 at L ¼ 500 for the ground state (blue) and for states drawn from the infinite-temperature ensemble (green). The three panels are consistent with the Ising glass order being present only in the red shaded region. In the rightmost panel, the weak quasiperiodic modulation leaves the excitations extended at all J=h; accordingly, the excited states are always paramagnetic, irrespective of ground-state ordering. In the center panel at large J=h, the excited-state order develops as the excitations localize across the AAlike transition described in Sec. IV B. The ground-state order parameter is completely insensitive to the excitedstate ordering. Finally, in the leftmost panel, the groundand excited-state order develop at the same coupling A J =h ¼ 2 because the ground-state symmetry-breaking phase transition and the localization transition of the excitations coincide, as discussed in Sec. V B.
Similar Ising glass order has previously been studied in disordered Ising chains [8,19,47]. The disordered glass and the quasiperiodic glass are diagnosed by the same eigenstate order parameters, and they show similar persistence of arbitrary initial longitudinal magnetization in quench experiments. However, the quasiperiodic realization is, in many respects, simpler to study and provides additional control knobs. First, unlike the disordered realization, there is a delocalization transition already in the noninteracting limit (discussed below). Second, sample-to-sample fluctuations are much larger in disordered systems than quasiperiodic. This is because of the outside influence in one dimension of rare spatial regions with atypically high or low effective disorder. These regions lead to broadened distributions of local magnetization, localization length of domain walls, and infinite randomness physics at the associated phase transitions. Since localization and the localized glass can take place without these features, the quasiperiodic glass is a simpler realization of localization protected order.

A. Excited-state transition
The Ising glass order is destroyed if either the ground state becomes paramagnetic or the domain-wall excitations delocalize. The central panel of Fig. 3 illustrates the latter transition at large J=h. In the leftmost panel of the same figure, on the other hand, the ground-state transition coincides with the delocalization of excitations (see Sec. V B). Below, we develop a picture of the transition of the central panel, which is entirely governed by the localization properties of the excitations. We leave more careful study of the other excited-state transition to future work [96].
At large J=h, the ground-state magnetization M 0 i is very close to 1 on each site i. In the excited states, the magnitude of the magnetization is reduced by the fluctuations of domain walls across site i. When the domain walls have typical localization length ξ, there are ∼ξ such relevant domain walls. Any domain wall localized further away merely flips the sign of M i without reducing its magnitude. Thus, the fluctuations of the parity of the domain walls to the left of i ultimately control the excited-state magnetization. Mathematically, where N <i is the number of domain walls to the left of site i. As the domain walls are noninteracting, we have where n α is the occupation of domain-wall eigenstate α and the indicator function is 1 if that domain wall is to the left of i. Within an excited eigenstate, N <i is thus a sum of independent random variables with mean and variance Here, P α <i is the probability that the domain wall in state α lies to the left of i. As the mean value of N <i only adjusts the overall sign of M i , we focus on the fluctuations δN <i to A. CHANDRAN and C. R. LAUMANN PHYS. REV. X 7, 031061 (2017) 031061-12 estimate the reduction of jM i j. In the localized regime, only those eigenmodes α with localization centers within a distance ξ of i contribute to these fluctuations as P α <i approaches 0 or 1 further away. At large ξ, δN <i becomes a Gaussian-distributed random variable with variance ∝ ξ. The magnitude of the magnetization is thus reduced by where a is related to the proportionality constant in the variance. Since in the large J=h limit, the localization transition is of AA type, we have ξ ∼ δ −1 (see Sec. IV B). Finally, this leads to an essential singularity in the excitedstate order at the transition, with a 0 a δ-independent constant. The quasiperiodic Ising order parameter O exc of Eq. (35) should approach jMj 2 at large L (and large spatial averaging window w) and accordingly inherits the essential singularity. Our numerics are consistent with this form but are inconclusive as it is difficult to distinguish an essential singularity from a shift in the apparent critical point.

VII. CONCLUSIONS
We have presented the first analytical study of localization-protected excited-state order without disorder. Incommensurate modulation of the exchange couplings leads to a large Ising glass phase in the canonical quantum Ising chain. By arguments similar to those presented in the context of disordered Ising chains [47], we expect that the glass survives the introduction of weak interactions so long as the localization length ξ of the domain walls is sufficiently short compared to the typical domain-wall density. Quasiperiodic modulation arises naturally in optical experiments; it also provides an analytic platform for further study of localization. In some ways, it is simpler than disordered localization, which is plagued by rareregion effects and Griffiths' phases [7,25,29,99,100]. Understanding the nature of the quasiperiodic localization transition, with and without interactions, may thereby cut to the heart of the phenomenon.
We also presented a theory of the melting transition for the excited-state order in the noninteracting case. Reducing the amplitude of modulation leads to delocalization of the domain walls, which may be accompanied by a groundstate symmetry-breaking transition. The divergence of the localization length ξ of domain walls leads to an essential singularity in the excited-state order. However, as interactions are likely to delocalize the system before ξ diverges, we expect the essential singularity to be cut off and the transition to qualitatively change. This is an especially interesting direction for future work. Our analysis is aided by our generalization of Aubry-André duality, AAA triality, which applies along the pure modulation axis (J ¼ 0) of the phase diagram in Fig. 1. The Ising-symmetry-breaking ground-state transition at A J =h ¼ 2 coincides with the self-trial point of this transformation. The triality requires that the wave functions are critically delocalized at all single-particle energies on one side of the transition (PM) and localized on the other (QPFM). The associated quantum critical point thus lies neither in the clean Ising universality class where the excitations are extended nor in the infinite randomness class of disordered systems where the excitations are fully localized. Intriguingly, numerics show that these critical properties persist away from the pure modulation line at low energies. This suggests that the entire phase boundary lies in this intermediate universality class, the quasiperiodic Ising class [96].
From a zero-temperature perspective, we expect the three ground-state phases to be stable to the inclusion of interactions as they are either gapped or localized. Whether interactions are irrelevant at the peculiar critical point discussed in the previous paragraph is an interesting open question. The transition is neither clean enough for a field theoretic renormalization treatment [101] nor disordered enough to obviously flow to infinite randomness under real-space renormalization [97]. Perhaps a coarse-graining treatment could be developed in the semiclassical limit of Q → 0, as has been done for the AA model [72].
Finally, we comment on the potential for studying the quasiperiodic Ising glass experimentally. Essentially, any quantum optical system that realizes a tunable Ising chain would be able to probe the excited-state Ising glass order by quench experiments. These include, for example, linear chains of trapped ions using hyperfine states as Ising degrees of freedom [34,64], Rydbergs trapped in optical tweezers [65,66], chains of trapped ions undergoing a zigzag transition [102,103], or ultracold atoms undergoing a staggering transition in a tilted lattice potential [104]. The simplest way to apply an incommensurate modulation in these systems is to modulate the position of the atoms or ions using an extra effective spatial potential, whether that be with a standing wave or with optical tweezers.
However, there are two classes of Ising model simulators: those where the Ising degree of freedom is spatial (such as at the zigzag transition) and those where it is internal (such as in the ion-trap Ising simulator). In the former, as the order parameter directly couples to spatial position, incommensurate modulation potentials locally break the Ising symmetry. This raises interesting questions regarding the Imry-Ma stability of excited-state symmetry-breaking order [105], which we leave for future work. In the latter Ising simulators, modulating the position can directly modulate couplings without introducing Ising odd terms (i.e., effective longitudinal fields). We expect these platforms to be able to realize the quasiperiodic Ising glass directly.
We include the elementary derivation of the contours of I at A h ¼ 0. Consider the differential (setting h ¼ 1 and assuming A J ; J > 0 in order to avoid writing absolute values throughout) Using contour integration on the unit circle z ¼ e iθ , where z AE ¼ −ðJ=A J Þ AE ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðJ=A J Þ 2 − 1 p . For J > A J , the pole at z þ lies inside the unit circle, while for J < A J , both z þ and z − lie outside. Explicitly, Thus, the contours of I behave differently in the regions above and below the diagonal J ¼ A J . Solving dI ¼ 0 above the diagonal leads to parabolic contours of the form Below, the contours are vertical. The functional form of I then follows immediately from explicit integration of Eq. (22) on the axes:

APPENDIX B: GENERAL COMMENSURATE ANALYSIS
We provide a detailed analysis of the spectral properties of the excitations in the TFIM with commensurate modulation Q ¼ 2πp=q, for p and q coprime integers. The incommensurate TFIM then arises by taking the limit p, q → ∞ such that Q=2π approaches an irrational value. For example, to study the incommensurate TFIM modulated with a wave vector corresponding to the golden mean, Q=2π ¼ ð1 þ ffiffi ffi 5 p Þ=2, one could take the sequence p n =q n ¼ F n =F n−1 , where F n is the nth Fibonacci number. Thouless [71,93] emphasized this approach to studying the original Aubry-André model, pointing out that q plays a role analogous to finite size in a scaling theory of the transition, as it determines the length over which the incommensurate model can be approximated by the commensurate one. Here, we review and generalize the commensurate analysis of Ref. [93] to the TFIM.
For commensurate modulation Q=2π ¼ p=q, Bloch's theorem implies that the single-particle energy spectrum of Eq. (8) is a function of the quasimomentum k x ∈ ½−π=q; π=qÞ with 2q bands. The 1D bands depend on the phase ϕ through the explicit Hamiltonian dependence on ϕ; we note that this dependence is periodic in ϕ → ϕ þ 2π=q, as the site labels can be shifted by an integer j → j − l, where lp ¼ 1ðmod qÞ to absorb such a shift. (As p and q are coprime, p has a multiplicative inverse modulo q.) It is natural to view the bands as 2D sheets over both k x and ϕ [this is the band structure of the associated 2D model, Eq. (15)], but one must remember that the bands of the actual 1D model correspond to a k x slice of the 2D bands at fixed ϕ.
Explicitly, the eigenvalue problem of Eq. (14) satisfies the difference equations: where h i ; J iþ1=2 , defined in Eq. (2), are periodic with period q. By Bloch's theorem, the solutions satisfy the twisted boundary conditions on a 2q-site chain: The In general, the characteristic polynomial depends on the phase Δ as well; we suppress this dependence as we specialize to the case A h ¼ 0 below. The analysis can be straightforwardly extended to nonzero A h . Note that CðE; k x ; ϕ h Þ is a polynomial of E of degree 2q. Since M is Hermitian, the polynomial is real (for E real): The bipartite (chiral) symmetry C, which takes ψ j →ð−1Þ j ψ j , anticommutes with M: CMC ¼ −M. Thus, the characteristic polynomial is an even function of E, and, accordingly, the eigenvalues come in AEE pairs within each k x sector. By explicit evaluation, C has the following structure: where all of the k x dependence lies in the energy-independent term. This dependence follows from the two k x -dependent terms in the explicit expansion of the determinant: which simplifies to the last term in Eq. (B6) using h j ¼ 1 and the definition In Appendix F, we evaluate this to be where T q is the Chebyshev polynomial of order q. Next, we evaluate the rest of the constant term, K 0 ðϕÞ. This is most easily accomplished by working at k x ¼ E ¼ 0, where M is antisymmetric and the determinant is the square of the Pfaffian. At h j ¼ 1, we have To make further progress, we need control of the higherorder coefficients K m ðϕÞ. For general couplings, these are not so easy to compute, although we set up some formalism exploring this in Appendix D. In the important special case of J ¼ 0, the triality discussed in Sec. IV E allows us to show that the K m are actually independent of ϕ for all m > 0; see Appendix C. With this simplification, we will be able to determine the energy-independent localization properties on the J ¼ 0 line; see Appendix E.

APPENDIX C: HIGHER-ORDER COEFFICIENTS WITH TRIALITY
In the Aubry-André models [56,71,93], the AA duality implies that CðE;k x ;ϕ;t; VÞ ¼ CðE;ϕ;−k x ;V;tÞ. Equating this order by order in E, K m ðϕ; t; VÞ ¼ K m ð−k x ; V; tÞ ð C1Þ for all m > 0. Differentiating with respect to ϕ, we see that K m ðϕÞ is independent of ϕ for m > 0. This leads to the energy independence of the AA localization transition by the logic described in Appendix E.
In the AAA models defined by Eq. (31), the triality transformation implies that CðE; k x ; ϕ; A 1 ; A 2 ; A 3 Þ ¼ CðE; Rðk x ; ϕÞ; A 2 ; A 3 ; A 1 Þ, where R is the linear transformation implementing the threefold rotation on the momentum space defined by k x , ϕ. We note that R is a geometric rotation by 2π=3, conjugated by an anisotropic scale transformation as our embedding of the honeycomb structure is into a rectangular lattice as in Fig. 2(c), rather than the geometrically symmetric embedding in Fig. 4. Nonetheless, where α, β are the appropriate matrix elements of R. Again, differentiating with respect to k x , it follows that K m is actually independent of its first argument ϕ in the AAA models.
In particular, this holds along the J ¼ 0 line of the incommensurate TFIM, which corresponds to A 3 ¼ h,

APPENDIX D: GENERAL HIGHER-ORDER COEFFICIENTS
In general, in the TFIM, the coefficients K m ðϕÞ are both ϕ dependent and nontrivial to evaluate for m > 0. Although we do not need any of the following formalism for the results used in the manuscript, we summarize here a few formulas for posterity.
To calculate the higher-order terms, it is helpful to consider the explicit representation of the determinant The factors of E come from diagonal matrix elements, so the coefficient K m comes from the permutations that hold 2m sites fixed but are otherwise off-diagonal, Moreover, since M only connects adjacent sites, only those permutations that permute within each diagonal block (from i j þ 1 to i jþ1 − 1) are nonzero. At fixed i 1 ; …; i 2m , the contribution to K m thus factors into a product of the determinants of the diagonal sub-blocks of M between rows i l and i lþ1 : where the l ¼ 2m block wraps around the corners of the matrix. Since each sub-block of M is antisymmetric, the determinant is nonzero only if i jþ1 − i j − 1 is even-that is, if the i j alternate between even and odd. Explicitly, Keep tracking of the number of i's, we find that the coefficients K m alternate in sign.
where Vði; jÞ ¼ 8 > > < > > : In this form, K m looks like The magnitude of K m is given by the canonical partition sum at temperature one for a system of length 2q with 2m domains of energy Vði; jÞ.
The number of wavelengths of J that fit in the chain of length 2q is p. If p is held fixed as q → ∞, then J iþ1=2 is smooth on the lattice scale a, as the wavelength of the incommensurate field q=p ≫ a. In this limit, the continuum approximation in which a ¼ 0 should be a good starting point for a semiclassical analysis that incorporates the leading effects of a small a=q. If on the other hand, p=q approaches a nonzero constant as q → ∞ (for example, the golden mean), then q=p ∼ OðaÞ. We leave further analysis of K m ðϕÞ in these limits to future work.

APPENDIX E: LOCALIZATION
PROPERTIES AT J = 0 From Appendix C, the characteristic polynomial may be written along the J ¼ 0 axis, where K m depends on couplings A J and h but not on k x or ϕ.
This form entails a remarkable geometric property of the 2q bands in the 2D band structure: They are all approximately the same shape up to shift and scale. This implies that the ratio r of the k x and ϕ dispersion of each band is independent of band number n. Thus, if this ratio approaches zero sufficiently rapidly with q as we approach the incommensurate limit by Q ¼ 2πp=q, the total dispersion in the k x must also go to zero, even when summed over all q bands (because the total 2D spectrum is bounded). This argument allows us to show that the spectrum is fully localized for A J > 2. For A J < 2, the r approaches a constant as q → ∞, which indicates the presence of critically delocalized states, also independent of energy. See Ref. [93] for a similar analysis in a simpler model. Suppose we identify the nth zero of CðEÞ at fixed ðk x ; ϕÞ ¼ ð0; 0Þ: CðE 0 n ; 0; 0Þ ¼ 0: The dispersion of the nth band is then determined by following this zero while varying k x and ϕ. To linear order in E − E 0 n , we have 0 ¼ CðE 0 n ; k x ; ϕÞ þ ðE n ðk x ; ϕÞ − E 0 n ÞC 0 n ; ðE3Þ where C 0 n ¼ ½ð∂CÞ=∂Ej E 0 n is a constant independent of the 2D momenta from Eq. (E1). Rearranging this equation, the nth band has dispersion E n ðk x ; ϕÞ ¼ −1 C 0 n ð−1Þ q ½PðϕÞ 2 − 2PðϕÞ cosðk x qÞ þ E n ; where E n is an n-dependent shift independent of momenta. Thus, each band has the same shape up to shift E n and scale 1=C 0 n , up to higher-order corrections in the E dependence of C.
Although the absolute position and scale of each band depends, in a detailed way, on the parameters, the ratio of the bandwidth in the k x and ϕ directions is independent of n and therefore easy to compute. The bandwidth in the k x direction is given by the total variation of Eq. (E4) at fixed ϕ, Δ n;k x ðϕÞ ¼ 4jPðϕÞj=C 0 n : Maximizing over ϕ gives Δ n;k x ¼ 16ðA J =2Þ q =C 0 n ; ðE6Þ where we have assumed q is even (the result differs by an unimportant factor of 2 for odd q). Similar elementary considerations yield the maximum variation in the ϕ direction at fixed k x , where again we have assumed q is even. Taking the ratio of the bandwidths in the k x and ϕ directions, which holds for all bands n. From Eq. (E8), it follows that the incommensurate TFIM is localized at all energies in the incommensurate limit (q → ∞) for A J > 2. The total dispersion in the k x direction at fixed ϕ (summed across bands) is exponentially smaller than that in the ϕ direction at fixed k x . Since the total 2D variation is upper bounded by a function that, at most, increases linearly with q [as each band has, at most, Oð1Þ bandwidth], the variation in the ϕ direction (which is exponentially larger in q than that in the k x direction) provides the entire contribution as q → ∞. Thus, the total bandwidth in the k x direction goes to zero exponentially as q → ∞, so the 1D incommensurate TFIM for A J > 2 has a pure-point spectrum.
For A J < 2, the ratio approaches 1, and the variation of the bands is the same in both the k x and ϕ directions. In previously studied models, the 1D spectrum is critically delocalized whenever the ratio remains of Oð1Þ in the incommensurate limit [93]. We conjecture and have confirmed numerically that this holds for the TFIM as well at J ¼ 0 (i.e., that all states are critical at all energies). The coincidence of criticality and the order-one ratio of bandwidths has not been proven mathematically in any model as far as we know.

APPENDIX F: EVALUATION OF PRODUCT OF PERIODIC COUPLINGS
In this appendix, we evaluate the expression The derivation is identical to that in Appendix A of Ref. [93]; we include it here for completeness. First, PðϕÞ can be rewritten as The expression is periodic in ϕ with period 2π=q as the addition of 2π=q to ϕ yields a permutation of the terms in the product. Thus, the only terms that survive the product in Eq. (F2) contain q factors of e iϕ or q factors of e −iϕ or equal numbers of factors of e iϕ and e −iϕ . The first two cases give A J 2 q ð−1Þ pq ð2 cosðqϕÞÞ: The ϕ-independent term can be obtained by evaluating PðϕÞ minus the term above at ϕ ¼ 0: The rhs is a polynomial in J=A J with zeros at for j ¼ 0; …; q − 1. Using the definition of the Chebyshev polynomial T q of order q, T q ðxÞ ≡ cosðq arccosðxÞÞ; we see that T q ðJ=A J Þ ¼ ð−1Þ pqþ1 is a polynomial of J=A J with the same zeros as the rhs of Eq. (F4). Thus, the two must be proportional. Comparing the coefficient of ðJ=A J Þ q in the two expressions, we obtain Combining this expression with Eq. (F3) gives ðT q ðJ=A J Þ þ ð−1Þ pq cosðqϕÞÞ: ðF8Þ Using ð−1Þ pq ¼ ð−1Þ pþqþ1 when p and q are relatively coprime, we obtain Eq. (B9).

APPENDIX G: GROUND-STATE SYMMETRY BREAKING FROM THE COMMENSURATE ANALYSIS
In Appendix F, we evaluated PðϕÞ ¼ Q q−1 m¼0 J mþ1=2 when the exchange coupling is commensurate with the underlying lattice at wave number Q ¼ 2πp=q. In order to approach the incommensurate limit, we take p; q → ∞ in such a way that the wavelength 2π=Q approaches an irrational number (in units of the underlying lattice constant a ¼ 1). As PðϕÞ is independent of the ratio p=q, the limiting expression is independent of the value of the irrational wavelength 2π=Q. From the discussion in Sec. III, the behavior of PðϕÞ controls the ground-state phase diagram: If PðϕÞ increases (decreases) exponentially with q, then the system is in the ferromagnetic (paramagnetic) phase.
From the explicit expressions,