Localization and symmetry breaking in the quantum quasiperiodic Ising glass

TL;DR

This paper explores how quasiperiodic modulation in quantum Ising chains leads to a stable, long-range ordered glass phase at all energy densities, revealing novel localization and symmetry-breaking phenomena.

Contribution

It introduces the concept of a quasiperiodic Ising glass and analyzes its stability, phase diagram, and transition properties, extending Aubry-André duality to this context.

Findings

Existence of a stable quasiperiodic Ising glass phase.

Rich zero-temperature phase diagram with multiple ground states.

Unusual quantum Ising transition with intermediate properties.

Abstract

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.

Equations

$$\mathcal{H}$$

$$h_j$$

$$J_{j+1/2}$$

$${\sigma }_i^x{\sigma }_{i+1}^x$$

$${\sigma }_i^z$$

$$h^{\prime }$$

$$J^{\prime }$$

$${\varphi }^{\prime }$$

$${\gamma }_{2i}$$

$$\mathcal{H}=\frac{i}{2}\sum _j\left[J_{j+1/2}{\gamma }_{2j+1}{\gamma }_{2j+2}+h_j{\gamma }_{2j}{\gamma }_{2j+1}\right]$$

$$\mathcal{H}$$

$$c_{\alpha }^{†}$$

$$\{c_{\alpha }^{†},c_{\beta }^{†}\}=\{c_{\alpha },c_{\beta }\}=0,\{c_{\alpha }^{†},c_{\beta }\}={\delta }_{\alpha \beta }$$

$$⟨\psi |{\sigma }_i^x{\sigma }_j^x|\psi ⟩$$

$$G_{ij}^{\psi }$$

$$⟨\psi |{\sigma }_i^x{\sigma }_j^x|\psi ⟩$$

$${\tilde{\mathcal{H}}}_{1D}(\varphi )=\sum _{ij}{d}_i^{†}{H}_{ij}{d}_j$$

$$=i\sum _j\left[J_{j+1/2}{d}_{2j+1}^{†}{d}_{2j+2}+h_j{d}_{2j}^{†}{d}_{2j+1}\right]+h.c.$$

$${\tilde{\mathcal{H}}}_{2D}$$

$$+\frac{A_J}{2}{e}^{iQ(j+1/2)}{d}_{2j+1,k}^{†}{d}_{2j+2,k+1}$$

$$+\frac{A_J}{2}{e}^{-iQ(j+1/2)}{d}_{2j+1,k}^{†}{d}_{2j+2,k-1}$$

$$+\frac{A_h}{2}{e}^{i(Qj+\Delta )}{d}_{2j,k}^{†}{d}_{2j+1,k+1}$$

$$+\frac{A_h}{2}{e}^{-i(Qj+\Delta )}{d}_{2j,k}^{†}{d}_{2j+1,k-1})+h.c.$$

$${\mathcal{H}}_{AA}=\sum _j-t(d_j^{†}{d}_{j+1}+h.c.)-2V\cos (Qj+{\varphi }_V)d_j^{†}{d}_j.$$

$$\xi =1/\log |V/t|\sim |V/t-1{|}^{-1}$$

$${\Gamma }_L=\sum _{j=0}^{\infty }{\alpha }_j{\gamma }_{2j}$$

$${\alpha }_{j+1}$$

$$S$$

$$=\sum _{i=0}^{j-1}\log |h_i|-\sum _{i=0}^{j-1}\log |J_{i+1/2}|$$

$$S$$

$$I$$

$$J/h$$

$$A_J/h$$

$$\mathcal{O}=[|⟨0|{\sigma }_{L/4+j}^x{\sigma }_{3L/4+j}^x|0⟩|{]}_j$$

$$\mathcal{O}$$

$$=\sum _j|{\psi }_j^{\alpha }{|}^4.$$

Full text

I Introduction

Nearly sixty years ago, Anderson discovered that quenched disorder could localize quantum particles and thus prevent the transport necessary for equilibration in isolated systems Anderson (1958). The recent interest in the role of interactions Altshuler et al. (1997); Basko et al. (2006); Gornyi et al. (2005); Oganesyan and Huse (2007); Monthus and Garel (2010); Vosk and Altman (2013); Pekker et al. (2014); Pal and Huse (2010); Znidaric et al. (2008); Bardarson et al. (2012); Serbyn et al. (2013a); Bauer and Nayak (2013); Serbyn et al. (2013b); Swingle (2013); Nandkishore et al. (2014); Serbyn et al. (2014a); Iyer et al. (2013); Kjäll et al. (2014); Laumann et al. (2014); Luitz et al. (2015); Chandran et al. (2015); Tang et al. (2015); Michal et al. (2014); Agarwal et al. (2015); Vasseur et al. (2015); Singh et al. (2016); You et al. (2016); Potter et al. (2015) and rapid experimental developments in synthetic quantum systems Schreiber et al. (2015); Kondov et al. (2015); Bordia et al. (2016); Choi et al. (2016); Smith et al. (2016); Kaufman et al. (2016); Neill et al. (2016); Bordia et al. (2016) 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 Huse and Oganesyan (2013); Serbyn et al. (2013b); Imbrie (2014); Ros et al. (2015); Chandran et al. (2015); Monthus (2016); Rademaker and Ortuño (2016); the dynamics of entanglement has emerged as a unifying framework for understanding thermalization Znidaric et al. (2008); Bardarson et al. (2012); Serbyn et al. (2013a); and, the long-lived coherence of localized systems may serve as a resource for quantum information processing Serbyn et al. (2014b); Yao et al. (2015); Bahri et al. (2015).

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 Huse et al. (2013). The central idea may be illustrated in the 1D ferromagnetic Ising chain. Ferromagnetic order in the ground state is usually destroyed in excited states due to 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 Huse et al. (2013); Pekker et al. (2014) and has been observed numerically in small interacting chains Kjäll et al. (2014). Localization protection has also been argued to extend to a host of more exotic orders Huse et al. (2013); Chandran et al. (2014); Bahri et al. (2015); Potter and Vasseur (2016) and to periodically driven (Floquet) systems Khemani et al. (2016); von Keyserlingk and Sondhi (2016); Else and Nayak (2016); Roy and Harper (2016); Po et al. (2016).

Localization, however, does not require disorder, as was first recognized by Azbel Azbel (1979), and Aubry and Andrè (AA) Aubry and André (1980) in the single particle context. These authors discovered that sufficiently strong quasiperiodic potentials can localize a quantum particle. Refs. Iyer et al. (2013); Michal et al. (2014); Mastropietro (2015) 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 wavevectors. Accordingly, many experiments in such systems have now observed single-particle localization Dal Negro et al. (2003); Fallani et al. (2007); Roati et al. (2008); Lahini et al. (2009); Modugno (2010); Segev et al. (2013) and, more recently, have also pushed into the interacting regime and high excitation energy densities to provide evidence for the many-body localized phase Schreiber et al. (2015); Bordia et al. (2016, 2016).

As quasiperiodic systems can show both localized and delocalized behavior already in the 1D non-interacting 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 have been implemented in ion traps Islam et al. (2011); Smith et al. (2016) and Rydberg systems Glaetzle et al. (2015); Labuhn et al. (2016). 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 to the coupling, rather than field, modulated regime. Our rigorous analytic controls 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 many-body localized phase of bosons persists in Ref. Iyer et al. (2013). We discuss both potential experimental realizations and the stability to interactions further in the conclusion, 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.

The associated zero temperature quantum phase transitions lie in two distinct universality classes. We find that weak quasiperiodic modulation is irrelevant at the clean Ising transition, so that 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 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 $\nu=1$, as at 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 wavefunction criticality on the ‘pure modulation’, $J/h=0$ axis in Fig. 1, which has been observed numerically before in Satija and Doria (1989); Cai (2014).

The non-interacting quasiperiodic models of Azbel, Aubry-André and their generalizations have been extensively studied by mathematicians and physicists over the last thirty years for a variety of reasons Simon (1982); Thouless (1983); Wilkinson (1984); Sokoloff (1985); Schellnhuber and Urbschat (1987); Hiramoto and Kohmoto (1989); Igloi (1993); Chaves and Satija (1997); Abanov et al. (1998); Jitomirskaya (1999); Vidal et al. (1999); Kraus et al. (2012); Ganeshan et al. (2015). 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 Abrahams et al. (1979). Mathematically, the models have a surprisingly rich analytic structure, exhibiting dualities, critically delocalized phases with fractal spectra and connections to higher dimensional Hofstadter-type models Aubry and André (1980); Hofstadter (1976). 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 Ceccatto (1989); Hermisson et al. (1997); Luck (1993a). 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 Satija and Doria (1989) or equivalently, in the modulated p-wave superconductor DeGottardi et al. (2013); Cai et al. (2013); Cai (2014). The AAA triality we introduce provides an analytic framework for explaining numerical observations in these works.

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, our 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.1. 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 quasi-periodic transverse field Ising model (TFIM) is:

[TABLE]

The model is illustrated in Fig. 2(a). Here, $\sigma_{j}^{\alpha}$ are the Pauli matrices with $j\in\mathbf{Z}$ running over the sites in the chain and $\alpha=x,y,z$, $Q$ is the wavevector of the modulation in units where the lattice spacing is $a=1$, and the phases $\phi$ and $\phi+\Delta$ shift the positions of the maxima of the couplings relative to the underlying lattice. The wavevector $Q$ is commensurate with the underlying lattice if $Q/2\pi=p/q$ is rational, and is incommensurate otherwise. We choose the wavevector of the modulation of the Ising coupling and the transverse field to be the same for simplicity and as this is natural if the modulation arises from the same physical source (eg. an incommensurate laser potential).

Global symmetries: The quasi-periodic TFIM has several global symmetries. The eponymous Ising symmetry is given by $G=\prod_{i}\sigma_{i}^{z}$ – this is the symmetry which breaks spontaneously in the $T=0$ ferromagnetic and localized Ising glass phases. The Hamiltonian $\mathcal{H}$ is also symmetric under complex conjugation $K$, which is anti-unitary. Finally, for special values of the modulation phases $\phi$ and $\Delta$ the model is symmetric under reflections across site $k$, $j\to k-j$ or bond $k+1/2$, $j\to k+1/2-j$. For example, at $\phi=0$ and $\Delta=0$, $\mathcal{H}$ is symmetric under $j\to-j$.

Ising duality: Under the duality transformation

[TABLE]

$\mathcal{H}$ maps onto another incommensurate TFIM $\mathcal{H}^{\prime}$ with the role of the field and bond couplings interchanged (and different boundary conditions). Formally, the duality maps

[TABLE]

The duality swaps paramagnetic and ferromagnetic phases but leaves the dynamical nature of the bulk single-particle excitations (spectrum and wavefunction localization) invariant. This can be seen most easily from the fermionization (see below) of $\mathcal{H}$ and its dual $\mathcal{H}^{\prime}$, whose non-interacting 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:

[TABLE]

The $\gamma$ operators are Hermitian and satisfy the canonical anti-commutation relation $\{\gamma_{i},\gamma_{j}\}=2\delta_{ij}$. The transformation maps the TFIM to a quadratic Majorana chain with Hamiltonian:

[TABLE]

See Fig. 2(b).

The dynamical and symmetry-breaking properties of the TFIM Hamiltonian $\mathcal{H}$ follow from the properties of the single-particle Hamiltonian $H$ defined by

[TABLE]

Hermiticity of $\mathcal{H}$ requires that $H_{ij}$ is an imaginary anti-symmetric $2L\times 2L$ matrix where $L$ is the number of spins. Thus, the eigenvalues of $H$ come in $\pm e$ pairs corresponding to complex conjugate eigenmodes, $\psi$ and $\overline{\psi}$. Labeling the $L$ positive energy eigenmodes by the index $\alpha$, we can diagonalize $\mathcal{H}$ into the familiar form

[TABLE]

where

[TABLE]

The complex fermion operators $c,c^{\dagger}$ satisfy the usual anti-commutation relations

[TABLE]

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 non-trivial phase. The simplest way to detect the topological phase of the Majorana chain is with open boundary conditions, in which case the topologically non-trivial phase possesses a pair of zero energy Majorana modes localized at the boundaries of the chain Kitaev (2001). The many-body ground state 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=\prod_{j}\sigma^{z}_{j}$ maps to the fermionic parity operator $G=\prod_{j}(-i\gamma_{2j}\gamma_{2j+1})$, while the symmetry under complex conjugation forces $H$ to be bi-partite. The action of the Ising duality on the Majorana chain shifts all the site labels by a half: $j\to j-1/2$ as mentioned above.

As all eigenstates of $\mathcal{H}$ correspond to Slater determinant states of the fermions $\gamma$, they all satisfy Wick’s theorem. This allows evaluation of the spin-spin correlation function,

[TABLE]

as a Pfaffian of the fermionic Green function

[TABLE]

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 (eg. up to $L=1000$ in this work). Evaluating these correlators at large separation $|i-j|$ allows us to numerically extract the magnetization as

[TABLE]

where $M^{\psi}_{i}$ is the magnetization of spin $i$ in state $|\psi\rangle$.

The 2D model: The incommensurate TFIM has a second useful representation in terms of a 2D model of non-interacting complex fermions. To derive this, consider first a complex fermionic model with the same single particle Hamiltonian as Eq. (8) 111This is often referred to as a ‘doubled’ model of the Majorana fermions since it corresponds to two non-interacting copies. We denote the Hamiltonians of ‘doubled’ fermion models by $\tilde{\mathcal{H}}$:

[TABLE]

Above, $d_{i}$ destroys a fermion at site $i$ and $d_{i},d_{j}^{\dagger}$ satisfy the usual complex anti-commutation relations. On varying the phase $\phi$ between $[0,2\pi)$, we generate a family of distinct 1D Hamiltonians. Treating $\phi$ as the momentum along an extra dimension $y$ and inverting the Fourier transformation, we obtain a 2D tight binding model:

[TABLE]

whose spectrum at fixed y-momentum $\phi$ reproduces the spectrum of the 1D model $\tilde{H}_{1D}(\phi)$, see Fig. 2(c). The Hamiltonian $\tilde{H}_{2D}$ describes a translation invariant hopping model with uniform flux $Q$ piercing each two-site unit cell in Landau gauge. The flux associated with hopping cycles within a unit cell depends on $\Delta$. There are no vertical hops in the 2D model (ie. from $(j,k)$ to $(j,k\pm 1)$) in Fig. 2(c) because the 1D Hamiltonian in Eq. (14) is off-diagonal.

The localization properties of the excitations of the 1D TFIM map onto those of the $\tilde{H}_{2D}$ in the Landau gauge at fixed y-momentum $\phi$. 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 will need several properties of the original Aubry-André hopping chain in the analysis of the quasiperiodic TFIM. The AA model has the following Hamiltonian Aubry and André (1980),

[TABLE]

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. Further, the localization length diverges at the critical point with exponent $\nu=1$. More precisely,

[TABLE]

Many features of the phase diagram follow Aubry-André duality. This duality corresponds to a $\pi/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 Hofstadter (1976). The $\pi/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 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 $\pi/2$. While one can embed the 2D model into a larger class of 2D models which have the requisite vertical bonds, these would in general need to be staggered in the $y$-direction and thus they do not correspond to 1D incommensurate chains. The $\pi/2$ rotation in the larger model space therefore does not define a duality map on the incommensurate TFIM.

This is both blessing and 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 will 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, even the ground state phase diagram has not been studied. The full dynamical phase diagram lives in an eight-dimensional space, parameterized by $h$, $J$, $A_{h}$, $A_{J}$, $Q$, $\phi_{h}$, $\Delta$ and the excitation energy $e$; after normalizing the units of energy, there are still 7 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 turns out to already be rather interesting. In units where $h=1$, the relevant parameters controlling the phase diagram are then $J$, $A_{J}$, $Q$, $\phi$ and the excitation energy $e$ (note that $\Delta$ has no effect when $A_{h}=0$). Most of the analytic features we derive hold for any incommensurate $Q$ and $\phi$, but we focus our numerical results (notably for spectra) on $Q/2\pi=(\sqrt{5}+1)/2$, the golden mean. By Ising duality, our analysis also produce the phase diagram at $A_{J}=0$.

III Ground state symmetry breaking

The quasiperiodic Ising model has three ground state phases: a paramagnetic (PM)phase, a ferromagnetic (FM) phase and a quasiperiodically alternating ferromagnetic (QPFM) phase. The two latter phases spontaneously break the Ising symmetry and the spin-spin correlation function $\langle\sigma^{x}_{i}\sigma^{x}_{j}\rangle$ has long-range order. In the simple FM, all of the spins magnetize in the same direction (although the magnetization is not spatially uniform), while in the QPFM, the spatial mean magnetization is zero due to the presence of antiferromagnetic links where $J_{j}<0$ We refer to both Ising ordered phases as “ferromagnetic” as the ordering even in the QPFM is unfrustrated (ie. gauge equivalent to a $J_{j}>0$ model).

The clean model ($A_{J}=A_{h}=0$) spontaneously breaks the Ising symmetry for $|J|>|h|$. On general grounds, the corresponding gapped ferro- 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, $\Gamma_{L}$ and $\Gamma_{R}$, bound to the left and right ends of an open chain Kitaev (2001). To detect ground state symmetry breaking, we look for normalizable boundary modes.

For simplicity, we focus on the left edge of a semi-infinite 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:

[TABLE]

(The right mode would be localized on the odd sublattice.) Substituting in the eigenvalue equation $[\mathcal{H},\Gamma_{L}]=0$, we obtain:

[TABLE]

Taking a logarithm, we obtain

[TABLE]

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|h|/|J|$ clearly grows linearly with $j$; the sign of $\log|h|/|J|$ thus determines the normalizability of $\Gamma_{L}$. Modulation at irrational wavevector $Q$, causes $S_{j}$ to fluctuate, but 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,

[TABLE]

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,

[TABLE]

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 wavenumbers $2\pi p/q\to 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. (III). In Fig. 3, the three panels present the average magnitude of the magnetization squared on representative cuts through the phase diagram, calculated at $Q/2\pi=(1+\sqrt{5})/2$ for open chains of length $L=1000$. More precisely, we detect the ground state order through the correlator,

[TABLE]

which is non-zero in both the FM ($A_{J}<J$) and QPFM ($A_{J}>J$) phase. The square brackets $[\cdot]_{j}$ indicate averaging over a small spatial window ($j=-w,\cdots,w$) to smooth out the quasiperiodic modulation of the magnitude. In the thermodynamic limit $L\to\infty$, we expect

[TABLE]

where $[|M|]$ 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|h|/|J|$ over a period matters to the long-distance behavior. Fourth, the fluctuations of $S$ relative to its linear trend are $O(1)$ in distance $j$. This is significantly less than the $O(\sqrt{j})$ 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 and PM and on the vertical boundary between the QPFM and PM have different character. We will 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 wavefunctions. That is, 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. These 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\to\infty$, 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 $\alpha$ as

[TABLE]

In finite size studies, the scaling of the IPR in a given energy window, $\textrm{IPR}\sim 1/L^{\gamma}$, detects the dynamical phase. In the extended phase, $\gamma=1$; in the localized, $\gamma=0$; and, in the critical, $0<\gamma<1$. Formally, these phases correspond to spectra that are absolutely continuous, pure point and singular continuous (fractal) respectively. We will use both diagnostics in the following analysis.

IV.1 Clean Limit $A_{J}\to 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. (III)) lies in the same universality class and thus we expect the low energy 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 Fig. 5(a,b).

IV.2 Large $J\gg A_{J},h$

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,

[TABLE]

Up to a constant, the Hamiltonian for a domain wall is:

[TABLE]

This is simply the AA model, Eq. (16), with $t=h$ and $V=A_{J}$. Thus, the domain walls are extended at all energies for $A_{J}/h<1$ and localized at all energies for $A_{J}/h>1$ formally at $J\to\infty$. This explains the vertical asymptote of the extended to localized phase boundary at large $J$ in Fig. 1.

Note one can also see the emergence of the AA model in this limit by considering the associated 2D model. In the large $J\gg 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}$.

IV.3 The Corner $J,A_{J}\ll h$

This regime is in the bottom left corner in Fig. 1. The ground state is paramagnetic and is very close to the product state $|\uparrow\uparrow\uparrow\cdots\uparrow\rangle$. The excitations are spin flips. Analogously to Sec. IV.2, the effective Hamiltonian for the spin flips is given by,

[TABLE]

up to a constant. This off-diagonal Aubry-Andre model admits the AA duality and accordingly has energy independent localization properties. The model has been previously studied by Thouless and co-workers Thouless (1983); Han et al. (1994), who found two dynamical phases: extended for $A_{J}<J$, and critical for $A_{J}>J$. Asymptotically, this coincides with the diagonal dashed line near the origin of Fig. 1.

As an aside, the critical phase for $A_{J}>J$ is a special feature of the $A_{h}=0$ plane. At small $A_{h}$ and $\Delta\neq 0$, the critical phase is destroyed and the states for $A_{J}>J$ are fully localized Liu et al. (2015).

IV.4 Large $A_{J}\gg J,h$

The 2D model in the $A_{J}\to\infty$ limit reduces to a collection of decoupled vertical zig-zag wires (dashed orange in Fig. 2), two per original spin. Without $h$, this 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

[TABLE]

where $\phi$ is the $y$-momentum 222Technically, the $\pm$ energy eigenstates at momentum $\phi$ are delocalized on both wires at position $j$.. The energy dispersions are shifted by $Q$ between adjacent wire pairs. Thus, for incommensurate $Q\sim O(1)$, neighboring wire states at fixed momentum $\phi$ are typically non-degenerate, 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 closely mirrors arguments which can be made in the anisotropic Hofstadter problem.

Thus, we expect that all states are localized at large $A_{J}\gg J,h$ to the far right of Fig. 1. We note that the 2D model does not have enhanced symmetry in this limit and thus the analysis does not provide as complete control as it does in the Hofstadter model.

IV.5 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 tri-ality’ operation in analog with AA duality. The existence of this triality implies that the localization properties of the wavefunctions 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 critical for $A_{J}/h<2$ and localized for $A_{J}/h>2$. The self-‘tri-al’ 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 Satija and Doria (1989); DeGottardi et al. (2013); Cai et al. (2013); Cai (2014).

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\pi/3$ independently in each layer, about any site. These are physical rotations in the 2D model, which need to include gauge transformations for $\tilde{H}_{2D}$ to return to Landau gauge.

In the usual AA model, $\pi/2$ rotation in the associated 2D model leads to another model with couplings $V$ and $t$ swapped. Rotating both layers by $2\pi/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 Hamiltonian,

[TABLE]

The original 1D TFIM at $J=0$ corresponds to $A_{3}=h$, $A_{1}=A_{2}=A_{J}/2$. In the extended AAA family of models, the $2\pi/3$ rotation maps $\tilde{\mathcal{H}}_{AAA}$ to $\tilde{\mathcal{H}}_{AAA}^{\prime}$ with cyclically permuted couplings. This is the AAA triality.

With this triality in hand, it is possible to analyze the wavefunction properties of the AAA model in considerable detail. We relegate this analysis to the Appendices, as it entails a considerable calculational detour. The upshot is that the 1D incommensurate TFIM has a critical to localized transition, at all excitation energies, at the self-trial point, $A_{J}=2h$.

IV.6 Numerical support away from the limits

Away from the special lines discussed in the previous subsections, the localization properties of the wavefunctions are both $Q$ and energy dependent. We rely on numerics to confirm the features summarized in Fig. 1. We have extensively investigated $Q/2\pi=(\sqrt{5}+1)/2$ and checked the qualitative features for several other incommensurate wavevectors. The general features are:

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. 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. 3.

For $A_{J}/h>2$, all states are localized, consistent with the analytically proven behavior at $J=0$ (Sec. IV.5) and the analysis at large $A_{J}$ (Sec. IV.4) and large $J$ (Sec. IV.2). At large $J/h$, we observe the expected energy independent localization transition near $A_{J}/h=1$. 4. 4.

On the $J=0$ line at $A_{J}/h<2$, we confirm that all states are critical 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 zero temperature quantum phase transition between the ground state para- and ferro-magnetic 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 Fisher (1995).

At both transitions, the correlation length $\xi$ diverges with exponent $\nu=1$. From the relation $\xi\sim-1/I$, where $I$ is the integral governing the convergence of the Majorana boundary mode, Eq. (22), it is straightforward to show that:

[TABLE]

where $\delta$ is the deviation from the phase boundary at $I=0$. This is consistent with the Harris-Luck criterion Luck (1993b, a), which imposes that $\nu\geq 1/d=1$ for phase transitions in the presence of incommensurate modulation.

V.1 Clean Ising Transition

While $\nu=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

the dynamical critical exponent $z=1$ so that the gap $\Delta$ closes linearly with $\delta$:

[TABLE] 2. 2.

the low-energy excitations are extended, as these are the excitations which 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\leq J/h\leq 0.5$), the extended PM ($0.5\leq J/h\leq J_{c}/h=1.0625$), and the extended FM ($J_{c}/h\leq 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 wavefunctions are indicated qualitatively by the coloring of the states; quantitatively, the lower panel shows the scaling exponent $\gamma$ of the low energy IPR with system size,

[TABLE]

The lower left panel shows that $\gamma$ is near $1$ for the low energy 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>\approx 0.2$) are localized across the transition. Nevertheless, the gap closes linearly (center top) and the low-energy excitations remain extended nearby ($\gamma\approx 1$, in the bottom panel). Thus, the ground state symmetry-breaking transition is still in the clean quantum Ising universality class.

V.2 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.5. 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 driven 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$):

• •

The single particle gap closes on approaching the transition from the paramagnetic side with an exponent $z\approx 2$. However, the gap does not reopen on the symmetry breaking side of the transition.

• •

All excitations are localized in the symmetry breaking phase at $A_{J}/h>2$.

• •

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 $\gamma$ lies properly between [math] and $1$ and exhibits strong finite-size fluctuations.

• •

At higher energy (and $J/h>0$), the excitations can be localized even in the paramagnetic phase. We do not know whether the critical low energy 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 Chandran and Laumann , 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$; while in the latter, the excitations are localized on both sides of the transition and the dynamics are activated (roughly, $z\to\infty$). 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 decay of the boundary mode are governed by the same correlation length $\xi$ with exponent $\nu=1$ while in the infinite randomness transition these are governed by two distinct diverging length scales (with $\nu=2$ and $1$, respectively).

VI Excited state Ising glass order

We expect localized Ising glass order at all energy densities in the regions where

the ground state breaks Ising symmetry, and 2. 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\to\infty$ independent of $Q$.

Representative correlation functions: Fig. 6 shows the spin-spin correlation function $\langle\sigma^{x}_{L/2}\sigma^{x}_{L/2+d}\rangle$ 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 long-range order in the ground state is clearly detected by the blue curve which approaches a non-zero 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 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\to\infty$; instead it fluctuates on an order one scale depending on whether spin $L/2$ and spin $L/2+d$ are aligned or anti-aligned in the chosen infinite temperature state.

Frozen order parameter: The quasiperiodic Ising glass order is detected by an order parameter which generalizes the ground state order parameter, Eq. (24), to finite energy density states:

[TABLE]

Here, $[\cdot]_{E,j}$ indicates averaging over both eigenstates $|E\rangle$ in some energy window and over a small spatial window $j\in[-w,w]$ to suppress the quasiperiodic fluctuations. $\mathcal{O}_{exc}$ is non-zero as $L\to\infty$ only in states with long range 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 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 AA-like transition described in Sec. IV.2. The ground state order parameter is completely insensitive to the excited state ordering. Finally, in the leftmost panel, the ground and 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.2.

VI.1 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.2. 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 Chandran and Laumann .

At large $J/h$, the ground state magnetization $M^{0}_{i}$ is very close to one 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 $\xi$, there are $\sim\xi$ 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,

[TABLE]

where $N_{<i}$ is the number of domain walls to the left of site $i$. As the domain walls are non-interacting, we have

[TABLE]

where $n_{\alpha}$ is the occupation of domain wall eigenstate $\alpha$ 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

[TABLE]

and variance

[TABLE]

Here, $P^{\alpha}_{<i}$ is the probability that the domain wall in state $\alpha$ 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 $\delta N_{<i}$ to estimate the reduction of $|M_{i}|$. In the localized regime, only those eigenmodes $\alpha$ with localization centers within a distance $\xi$ of $i$ contribute to these fluctuations as $P^{\alpha}_{<i}$ approaches [math] or $1$ further away.

At large $\xi$, $\delta N_{<i}$ becomes a Gaussian distributed random variable with variance $\propto\xi$. The magnitude of the magnetization is thus reduced by,

[TABLE]

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 $\xi\sim\delta^{-1}$ (see Sec. IV.2). This leads finally to an essentially singularity in the excited state order at the transition,

[TABLE]

with $a^{\prime}$ a $\delta$-independent constant.

The quasiperiodic Ising order parameter $\mathcal{O}_{exc}$ of Eq. (35) should approach $|M|^{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 Huse et al. (2013), we expect that the glass survives the introduction of weak interactions, so long as the localization length $\xi$ 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 rare-region effects and Griffiths’ phases Vosk and Altman (2013); Potter et al. (2015); Agarwal et al. (2015); Gopalakrishnan et al. (2016); Khemani et al. (2016). 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 non-interacting case. Reducing the amplitude of modulation leads to delocalization of the domain walls which may be accompanied by a ground-state symmetry-breaking transition. The divergence of the localization length $\xi$ of domain walls leads to an essential singularity in the excited state order. However, as interactions are likely to delocalize the system before $\xi$ diverges, we expect the essential singularity to be cutoff 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 wavefunctions 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’ Chandran and Laumann .

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 Cardy (1996) nor disordered enough to obviously flow to infinite randomness under real space renormalization Fisher (1995). Perhaps a coarse-graining treatment could be developed in the semi-classical limit of $Q\to 0$, as has been done for the AA model Wilkinson (1984).

Finally, we would like to 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 Islam et al. (2011); Smith et al. (2016), Rydbergs trapped in optical tweezers Glaetzle et al. (2015); Labuhn et al. (2016), chains of trapped ions undergoing the zig-zag transition Enzer et al. (2000); Shimshoni et al. (2011) or ultracold atoms undergoing a staggering transition in a tilted lattice potential Simon et al. (2011). The simplest way to apply an incommensurate modulation in these systems is to modulate the position of the atoms/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 zig-zag 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 Imry and Ma (1975), which we leave for future work. In the latter Ising simulators, modulating the position can directly modulate couplings without introducing Ising odd terms (ie. effective longitudinal fields). We expect these platforms to be able to realize the quasiperiodic Ising glass directly.

Acknowledgements.

We would like to thank V. Oganesyan, D. Huse, A. Polkovnikov and A. Kaufman for helpful discussions. We thank the Kavli Institute for Theoretical Physics (KITP) in Santa Barbara for their hospitality during the early stages of this work and the National Science Foundation (NSF) under Grant No. NSF PHY11-25915 for supporting KITP. C.R.L. acknowledges support from the Sloan Foundation through a Sloan Research Fellowship and the NSF through Grant No. PHY-1656234. Note that any opinion, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF.

Appendix A Ground state phase boundary

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):

[TABLE]

Using contour integration on the unit circle $z=e^{i\theta}$,

[TABLE]

where $z_{\pm}=-\frac{J}{A_{J}}\pm\sqrt{\left(\frac{J}{A_{J}}\right)^{2}-1}$. 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,

[TABLE]

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,

[TABLE]

Below, the contours are vertical. The functional form of $I$ then follows immediately from explicit integration of (22) on the axes:

[TABLE]

Appendix B General commensurate analysis

In the appendices, we provide a detailed analysis of the spectral properties of the excitations in the TFIM with commensurate modulation $Q=2\pi p/q$, for $p$ and $q$ coprime integers. The incommensurate TFIM then arises by taking the limit $p,q\to\infty$ such that $Q/2\pi$ approaches an irrational value. For example, to study the incommensurate TFIM modulated with a wavevector corresponding to the golden mean, $Q/2\pi=(1+\sqrt{5})/2$, one could take the sequence $p_{n}/q_{n}=F_{n}/F_{n-1}$, where $F_{n}$ is the $n$’th Fibonacci number. Thouless Thouless (1983); Han et al. (1994) emphasized this approach to studying the original Aubry-André model, pointing out the $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. Han et al. (1994) to the TFIM.

For commensurate modulation $Q/2\pi=p/q$, Bloch’s theorem implies that the single-particle energy spectrum of Eq. (8) is a function of the quasi-momentum $k_{x}\in[-\pi/q,\pi/q)$ with $2q$ bands. The 1D bands depend on the phase $\phi$ through the explicit Hamiltonian dependence on $\phi$; we note that this dependence is periodic in $\phi\to\phi+2\pi/q$, as the site labels can be shifted by an integer $j\to j-l$ where $lp=1(\textrm{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 $\phi$ (this is the band structure of the associated 2D model, Eq. (15)), but one must remember that the bands of the actual 1D model corresponds to a $k_{x}$ slice of the 2D bands at fixed $\phi$.

Explicitly, the eigenvalue problem of Eq. (14) satisfies the difference equations:

[TABLE]

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:

[TABLE]

The eigenspectrum of the infinite set of equations in Eq. (B) follow from the eigenspectrum of the finite matrix

[TABLE]

with $J_{q-1/2}=J_{-1/2}$. The characteristic polynomial of $M$ is:

[TABLE]

In general, the characteristic polynomial depends on the phase $\Delta$ as well; we suppress this dependence as we specialize to the case $A_{h}=0$ below. The analysis can be straightforwardly extended to non-zero $A_{h}$.

$C(E;k_{x},\phi_{h})$ is a polynomial of $E$ of degree $2q$. Since $M$ is Hermitian, the polynomial is real (for $E$ real):

[TABLE]

The bipartite (chiral) symmetry $\mathcal{C}$ which takes $\psi_{j}\to(-1)^{j}\psi_{j}$ anticommutes with $M$: $\mathcal{C}M\mathcal{C}=-M$. Thus, the characteristic polynomial is an even function of $E$,

[TABLE]

and accordingly that the eigenvalues come in $\pm E$ pairs within each $k_{x}$ sector.

By explicit evaluation, $C$ has the following structure,

[TABLE]

where all of the $k_{x}$ dependence lies in the the energy independent term. This dependence follows from the two $k_{x}$ dependent terms in the explicit expansion of the determinant:

[TABLE]

which simplifies to the last term in Eq. (52) using $h_{j}=1$ and the definition,

[TABLE]

In Appendix F, we evaluate this to be

[TABLE]

where $T_{q}$ is the Chebyshev polynomial of order $q$.

Next, we evaluate the rest of the constant term, $K_{0}(\phi)$. 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:

[TABLE]

To make further progress, we need control of the higher order coefficients $K_{m}(\phi)$. For general couplings, these are not so easy to compute, although we setup some formalism exploring this in App. D. In the important special case of $J=0$, the triality discussed in Sec. IV.5 allows us to show that the $K_{m}$ are actually independent of $\phi$ for all $m>0$, see App. C. With this simplification, we will be able to determine the energy independent localization properties on the $J=0$ line, see App. E.

Appendix C Higher order coefficients with triality

In the Aubry-Andre models Aubry and André (1980); Thouless (1983); Han et al. (1994), the AA duality implies that $C(E;k_{x},\phi;t,V)=C(E;\phi,-k_{x};V,t)$. Equating this order by order in $E$,

[TABLE]

for all $m>0$. Differentiating with respect to $\phi$, we see that $K_{m}(\phi)$ is independent of $\phi$ for $m>0$. This leads to the energy independence of the AA localization transition by the logic described in App. E.

In the AAA models defined by Eq. IV.5, the triality transformation implies that $C(E;k_{x},\phi;A_{1},A_{2},A_{3})=C(E;R(k_{x},\phi);A_{2},A_{3},A_{1})$, where $R$ is the linear transformation implementing the three-fold rotation on the momentum space defined by $k_{x},\phi$. We note that $R$ is a geometric rotation by $2\pi/3$ conjugated by an anisotropic scale transformation as our embedding of the honeycomb structure is into a rectangular lattice as in Fig. 2c, rather than the geometrically symmetric embedding in Fig. 4. Nonetheless,

[TABLE]

where $\alpha,\beta$ are the appropriate matrix elements of $R$. Again, differentiating with respect to $k_{x}$, it follows that $K_{m}$ is actually independent of it’s first argument $\phi$ in the AAA models.

In particular, this holds along the $J=0$ line of the incommensurate TFIM, which corresponds to $A_{3}=h,A_{1}=A_{2}=A_{J}/2$.

Appendix D General higher order coefficients

In general in the TFIM, the coefficients $K_{m}(\phi)$ are both $\phi$-dependent and non-trivial 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 formulae for posterity.

To calculate the higher order terms it is helpful to consider the explicit representation of the determinant (with $B=M-E$),

[TABLE]

The factors of $E$ come from diagonal matrix elements, so the coefficient $K_{m}$ comes from the permutations which hold $2m$ sites fixed but are otherwise off-diagonal.

[TABLE]

Moreover, since $M$ only connects adjacent sites, only those permutations which permute within each diagonal block (from $i_{j}+1$ to $i_{j+1}-1$) are non-zero. At fixed $i_{1}\ldots i_{2m}$, the contribution to $K_{m}$ thus factors into a product of the determinants of the diagonal subblocks of $M$ between rows $i_{l}$ and $i_{l+1}$:

[TABLE]

where the $l=2m$ block wraps around the corners of the matrix. Since each subblock of $M$ is antisymmetric, the determinant is non-zero only if $i_{j+1}-i_{j}-1$ is even – that is, if the $i_{j}$ alternate between even and odd. Explicitly,

[TABLE]

Keep tracking of the number of $i$’s, we find that the coefficients $K_{m}$ alternate in sign.

For $A_{h}=0$,

[TABLE]

where

[TABLE]

In this form, $K_{m}$ looks like,

[TABLE]

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$ which fit in the chain of length $2q$ is $p$. If $p$ is held fixed as $q\to\infty$, then $J_{i+1/2}$ is smooth on the lattice scale $a$, as the wavelength of the incommensurate field $q/p\gg a$. In this limit, the continuum approximation in which $a=0$ ought to be a good starting point for a semi-classical analysis that incorporates the leading effects of a small $a/q$. If on the other hand, $p/q$ approaches a non-zero constant as $q\to\infty$ (for example, the golden mean), then $q/p\sim O(a)$. We leave further analysis of $K_{m}(\phi)$ in these limits to future work.

Appendix E Localization properties at $J=0$

From App. C, the characteristic polynomial may be written along the $J=0$ axis,

[TABLE]

where $K_{m}$ depends on couplings $A_{J}$, $h$ but not on $k_{x}$ or $\phi$.

This form entails a remarkable geometric property of the $2q$ bands in the 2D bandstructure: they are all approximately the same shape up to shift and scale. This implies that the ratio $r$ of the $k_{x}$ and $\phi$ dispersion of each band is independent of band number $n$. Thus, if this ratio approach zero sufficiently rapidly with $q$ as we approach the incommensurate limit by $Q=2\pi p/q$, the total dispersion in the $k_{x}$ must also go to zero even 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\to\infty$, which indicates the presence of critical states, also independent of energy. See Han et al. (1994) for similar analysis in a simpler model.

Suppose we identify the $n$’th zero of $C(E)$ at fixed $(k_{x},\phi)=(0,0)$:

[TABLE]

The dispersion of the $n$’th band is then determined by following this zero while varying $k_{x}$ and $\phi$. To linear order in $E-E_{n}^{0}$, we have

[TABLE]

where $C^{\prime}_{n}=\left.\frac{\partial C}{\partial E}\right|_{E_{n}^{0}}$ is a constant independent of the 2D momenta from Eq. (E). Rearranging, this implies the $n$’th band has dispersion

[TABLE]

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^{\prime}_{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 $\phi$ directions is independent of $n$ and therefore easy to compute. The bandwidth in the $k_{x}$ direction is given by the total variation of the Eq. 74 at fixed $\phi$.

[TABLE]

Maximizing over $\phi$ gives,

[TABLE]

where we have assumed $q$ even (the result differs by an unimportant factor of 2 for $q$ odd). Similar elementary considerations yield the maximum variation in the $\phi$ direction at fixed $k_{x}$,

[TABLE]

where again we have assumed $q$ is even. Taking the ratio of the bandwidths in the $k_{x}$ and $\phi$ directions,

[TABLE]

which holds for all bands $n$.

From Eq. 78, it follows that the incommensurate TFIM is localized at all energies in the incommensurate limit ($q\to\infty$) for $A_{J}>2$. The total dispersion in the $k_{x}$ direction at fixed $\phi$ (summed across bands) is exponentially smaller than that in the $\phi$ 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 $\phi$ direction (which is exponentially larger in $q$ than that in the $k_{x}$ direction) provides the entire contribution as $q\to\infty$. Thus, the total bandwidth in the $k_{x}$ direction goes to zero exponentially as $q\to\infty$, so that 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 $\phi$ directions. In previously studied models, the 1D spectrum is critical whenever the ratio remains of $O(1)$ in the incommensurate limit Han et al. (1994). We conjecture and have confirmed numerically that this holds for the TFIM as well at $J=0$ (ie. 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 are aware.

Appendix F Evaluation of product of periodic couplings

In this appendix, we evaluate the expression:

[TABLE]

The derivation is identical to that in Appendix A of Ref. Han et al. (1994); we include it here for completeness. First, $P(\phi)$ can be re-written as:

[TABLE]

The expression is periodic in $\phi$ with period $2\pi/q$ as the addition of $2\pi/q$ to $\phi$ yields a permutation of the terms in the product. Thus, the only terms that survive the product in Eq. (80) contain $q$ factors of $e^{i\phi}$ or $q$ factors of $e^{-i\phi}$ or equal numbers of factors of $e^{i\phi}$ and $e^{-i\phi}$. The first two cases give:

[TABLE]

The $\phi$ independent term can be obtained by evaluating $P(\phi)$ minus the term above at $\phi=0$:

[TABLE]

Rearranging:

[TABLE]

The RHS is a polynomial in $J/A_{J}$ with zeros at:

[TABLE]

for $j=0,\ldots,q-1$. Using the definition of the Chebyshev polynomial $T_{q}$ of order $q$:

[TABLE]

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. (F). Thus, the two must be proportional. Comparing the coefficient of $(J/A_{J})^{q}$ in the two expressions, we obtain:

[TABLE]

Combining this expression with Eq. (81) gives:

[TABLE]

Using $(-1)^{pq}=(-1)^{p+q+1}$ when $p$ and $q$ are relatively co-prime, we obtain Eq. (55).

Appendix G Ground state symmetry-breaking from the commensurate analysis

In Appendix F, we evaluated $P(\phi)=\prod_{m=0}^{q-1}J_{m+1/2}$ when the exchange coupling is commensurate with the underlying lattice at wavenumber $Q=2\pi p/q$. In order to approach the incommensurate limit, we take $p,q\to\infty$ in such a way that the wavelength $2\pi/Q$ approaches an irrational number (in units of the underlying lattice constant $a=1$). As $P(\phi)$ is independent of of the ratio $p/q$, the limiting expression is independent of the value of the irrational wavelength $2\pi/Q$. From the discussion in Sec. III, the behavior of $P(\phi)$ controls the ground state phase diagram: if $P(\phi)$ increases (decreases) exponentially with $q$, then the system is in the ferromagnetic (paramagnetic) phase.

From the explicit expressions:

[TABLE]

the asymptotic forms of $P(\phi)$ at large $q$ follow:

[TABLE]

Note that $h=1$ in the expression above. Reinstating the $h$-dependence, we obtain the boundary between FM and PM ground states to be at $A_{J}/h=2$ for $0\leq J/h\leq 2$ and at $J/h+\sqrt{(J/h)^{2}-(A_{J}/h)^{2}}=2$ for $A_{J}/h<2$. Simplifying the latter expression, the critical value of the exchange coupling is $(J/h)_{c}=1+(A_{J}/2h)^{2}$ for $A_{J}/h\leq 2$. This is in perfect agreement with Eq. (III).