Non-power-law universality in one-dimensional quasicrystals
Attila Szabó Ulrich Schneider





Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, United Kingdom


    We have investigated scaling properties of the Aubry–André model and related one-dimensional quasiperiodic Hamiltonians near their localisation transitions.
    We find numerically that the scaling of characteristic energies near the ground state, usually captured by a single dynamical exponent, does not obey a power law relation. Instead, the scaling behaviour depends strongly on the correlation length in a manner governed by the continued fraction expansion of the irrational number β describing incommensurability in the system. This dependence is, however, found to be universal between a range of models sharing the same value of β. For the Aubry–André model, we explain this behaviour in terms of a discrete renormalisation group protocol which predicts rich critical behaviour.
    This result is complemented by studies of the expansion dynamics of a wave packet under the  Aubry–André model at the critical point. Anomalous diffusion exponents are derived in terms of multifractal (Rényi) dimensions of the critical spectrum; non-power-law universality similar to that found in ground state dynamics is observed between a range of critical tight-binding Hamiltonians.


Non-power-law universality in one-dimensional quasicrystals
    Ulrich Schneider
    
===========================================================





§ INTRODUCTION




Quasiperiodic structures, which are long-range ordered without being periodic, represent a rich and fascinating middle ground between ordinary periodic crystals and disordered systems. They were first discovered among aperiodic tilings of the plane, the best known of which is the fivefold symmetric Penrose tiling <cit.>. Interest in quasiperiodicity within the physics community was sparked by the discovery of quasicrystals by Shechtman <cit.> and the equivalence between Landau levels on two-dimensional lattices and a one-dimensional quasiperiodic chain <cit.>.

Recently, quasiperiodic structures became popular in ultracold atom experiments as a proxy for random potentials  in the study of disordered quantum gases, Bose glasses, and many-body localisation, as they can conveniently be realised by superimposing two incommensurate optical lattices <cit.>. 

Quasiperiodic tilings also lie at the heart of recent results in the study of quantum complexity, such as the proof of the undecidability of the spectral gap <cit.>.

Quasiperiodicity gives rise to a range of unusual behaviour including critical spectra and multifractal eigenstates away from phase transitions <cit.> and localisation transitions at a finite modulation of the on-site potential <cit.>. 

In this paper, we investigate localisation transitions of one-dimensional quasiperiodic systems, in particular the tight-binding Aubry–André model, also known as the Harper model <cit.>:

    H = -J ∑_j( a_j^† a_j+1 + H.c.) - Jλ∑_jcos(2πβ j) a_j^† a_j

and related models. Here β∉ℚ and λ are the incommensurate wave number and dimensionless amplitude of the on-site energy modulation, respectively, J is the hopping matrix element,
and a^†_j is a bosonic creation operator on the jth lattice site. Since the integer part of β is irrelevant, we assume β<1.

This model is known to undergo a localisation transition at λ=2 for any irrational value of β <cit.>: below this critical value, all eigenstates are extended while above it, they are exponentially localised. This is a consequence of Aubry duality: under the Fourier transform

    b_k = 1/√(𝒩)∑_n exp(2π iβ kn) a_n,

(<ref>) turns into another Aubry–André Hamiltonian in momentum space with λ changed to 4/λ and all energies rescaled by a factor of λ/2 <cit.>: λ=2 is the fixed point of this transformation. 


It is well known that the spectra of one-dimensional quasiperiodic Hamiltonians are hierarchical <cit.>, meaning they contain a hierarchy of progressively smaller gaps. In the case of tight-binding models, the spectrum is bounded and its entire structure is governed by the continued fraction expansion of the incommensurate ratio β <cit.>,

    β   = 1/n_1+β_1 = 1n_1+1n_2+β_2 = … = 1n_1+1n_2+1⋱
       =:[0;n_1,n_2,n_3,…],

where n_k are integers and the irrational residuals β_k are between 0 and 1. 
The hierarchical spectrum of these Hamiltonians
can be constructed as the limiting case of periodic superlattices with increasing periods N_k  described by rational approximants of β≈ M_k/N_k = [0;n_1,…,n_k], as  discussed in detail below. In going from the kth-order superlattice to the (k+1)st, each band of the spectrum is split into n_k+1 new ones <cit.>, see Fig. <ref>. In this manner, the periods of these approximant superlattices, N_k, act as `microscopic length scales' of the problem: the structure of the spectrum and eigenstates of the Hamiltonian on length scales around each N_k is controlled solely by the coefficient n_k. As a consequence, the spectrum is self-similar if and only if the continued fraction expansion of β is periodic. Furthermore, its hierarchy is topologically protected under smooth deformations between different models sharing the same value of β <cit.>.


In a continuous phase transition, the correlation length ξ diverges at the transition point. In conventional  disordered or crystalline systems, the effect of microscopic structure becomes immaterial once ξ is much larger than all microscopic scales of the system. Therefore, their behaviour near the phase transition is described by scale-invariant functions, that is, power laws <cit.>. 

In quasiperiodic systems, however, such a scaling regime is never reached due to the increasingly large `microscopic' length scales N_k discussed above. Instead, the behaviour of the system  is governed by scaling properties of the critical spectrum and eigenstates at length scales close to ξ as it diverges, which in turn depends on the coefficients n_k.  While the connection between the structure of the spectrum and the length scales of the system has tacitly been known, its effects on phase transitions were not discussed, mostly because all numerical and most analytical studies focused on β's of particularly simple continued fraction expansions, such as the golden mean ϕ^-1=(√(5)-1)/2=[0;1] <cit.>. (The overbar denotes a periodic continued fraction, e.g., [0;1,2,3] = [0;1,2,3,2,3,…].)


In this paper, we explore some consequences of this non-power-law critical behaviour on the localisation transition of the single-particle Aubry–André model (<ref>) using exact diagonalisation and renormalisation group arguments. In particular, we investigate the critical dynamics of the model for different values of β and demonstrate that power-law behaviour emerges only when the continued fraction expansion of β is periodic.

Section <ref> reviews the origins of hierarchical spectra in quasiperiodic systems and presents a renormalisation group treatment of the Aubry–André model based on Ref. Suslov1982. In Sec. <ref>, we discuss the scaling of energy scales near the ground state; Sec. <ref> deals with fractal properties of the spectrum and quench dynamics at critical points. In both cases, we find  equivalent behaviour for different models sharing the same β. We understand this equivalence as a novel kind of universality, distinct from power-law thermodynamic universality, but similarly protected by symmetries of the underlying systems. Conclusions are presented in Sec. <ref>.





















§ STRUCTURE OF THE SPECTRUM AND EIGENSTATES






 §.§ Structure of the critical spectrum







We consider a one-dimensional quasiperiodic tight-binding system characterised by the incommensurate ratio β. For simplicity, we assume that the continued fraction terms of β, n_k, are all very large; however, the qualitative structure of the spectrum remains the same for all n_k≥ 2  <cit.>[11]The case of n_k=1 is special. It implies β_k-1>1/2 which can be replaced with 1-β_k-1 without changing the resulting structure: the first continued fraction term of this number is, however, greater than 1.. 

Now, as discussed in Sec. <ref>, the structure of the spectrum can be described in terms of a sequence of periodic superlattices described by M_k/N_k=[0;n_1,…,n_k], which are the closest rational approximations of β in the sense that <cit.>

    |N_kβ-M_k|<|Nβ-M|  ∀ M,N∈ℤ, 0<N<N_k.

At the first step of this protocol, M_1/N_1=1/n_1: Bloch's theorem applies to the superlattice of period n_1, resulting in a spectrum consisting of n_1 subbands with continuous dispersion. At the next step, the period of the superlattice and thus the number of subbands is N_2=n_1n_2+1≈ n_1n_2 [Not all basis states of the system are accounted for in this approximation. This is compensated by the appearance of an additional kth order subband in the middle of each (k-2)nd order subband (see Fig. <ref>), also described by the incommensurate ratio β_k <cit.>. The approximation, however, provides a good description of states near the edges of the spectrum <cit.>.]. Since the approximation to β changes very little, the spectrum is still dominated by the n_1 first-order bands, each now split into n_2 narrower subbands, see Fig. <ref>. As further continued fraction terms are taken into account, more and more narrow subbands are formed, each time by splitting existing subbands into n_k new ones. 





The formation of this hierarchical structure can be understood in terms of a discrete renormalisation group procedure <cit.>. Creating first-order subbands can be taken as renormalising length scales by a factor of β^-1≈ n_1: the new `lattice sites' correspond to approximate Wannier states located at each minimum of the quasiperiodic modulation (see Fig. <ref>). Since the modulation period is incommensurate to the lattice spacing, the jth renormalised lattice site will have a phase shift 2πβ_1 j compared to the original lattice sites. This results in a quasiperiodic modulation of incommensurate ratio β_1 in the effective Hamiltonian of each subband. These Hamiltonians can now be renormalised by a factor of β_1^-1≈ n_2, giving rise to second-level subbands modulated with the new incommensurate ratio β_2: repeating such steps indefinitely constructs the entire spectrum. It can be shown <cit.> that for the Aubry–André model with n_k≫1, the renormalised on-site potential remains purely sinusoidal, that is, the effective Hamiltonian of each subband is of Aubry–André form with incommensurate ratio β_k at the kth step.








 §.§ Behaviour in the extended and localised regimes



At λ=2 in the Aubry–André model, the hierarchical structure of the spectrum described above is manifest at all energy scales and therefore all length scales. Away from the critical point, the correlation length ξ of the system becomes finite. On this length scale, either the potential or the kinetic energy term of (<ref>) becomes irrelevant, resulting in an either absolutely continuous spectrum and extended states for λ<2 or a dense point spectrum and exponentially localised states for λ>2. The crossover between the critical spectrum and the extended or localised spectra can be demonstrated using the thermal entropy of a single particle in a canonical ensemble:

    S    = -k_B∑_i p_i log p_i    
        p_i    = e^-E_i/(k_BT)/∑_j e^-E_j/(k_BT) .

At temperature T, exp(S/k_B) is a measure of the number of states up to ∼ k_BT above  the ground state. S(T) is plotted for different values of λ≤2 in Fig. <ref>. At criticality, each renormalisation step defines a new energy scale resulting in an infinite staircase structure. For λ≠2, such stairs persist down to energy scales corresponding to lengths on the order of ξ. Below this scale, the stairs smooth out and the scaling of entropy with temperature approaches that expected for an unmodulated tight-binding chain.





In the localised phase λ>2, ξ is normally identified with the localisation length of the wave function envelope which can be calculated without detailed analysis of the wave function <cit.>: in the Aubry–André case <cit.>,


    ξ (λ>2) = (logλ/2)^-1

for all eigenstates and all values of β. Due to Aubry duality, the structure of the Aubry–André spectrum for modulation amplitudes λ and 4/λ is identical save for an overall rescaling <cit.>. This implies that the length scale ξ where the crossover happens in the two cases is the same, giving

    ξ (λ) = |logλ/2|^-1.







The crossover between critical and extended or localised behaviour is also manifest in the structure of the wave functions. At λ=2, nontrivial structure appears at all length scales: away from criticality, this structure is only manifest up to length scales ≈ξ beyond which the density distribution is either dominated by exponential decay or becomes uniform. This is demonstrated for the ground state in Fig. <ref> which also confirms the localisation length given by (<ref>).









 §.§ Analytic scaling theory





In this section, we present a full renormalisation group treatment of the Aubry–André model based on Ref. <cit.>. This treatment becomes exact in the limit when all continued fraction terms of β are large, that is, β,β_1,…≪1.


As in Sec. <ref>, we start by approximating β with 1/n_1, that is, we consider the following periodic Hamiltonian:

    H = -   J ∑_j( a_j^† a_j+1 + H.c.) 
    
        -    Jλ∑_jcos(2π/n_1 (j-ϕ)) a_j^† a_j,

where ϕ is a well-defined global spatial offset. 
The spectrum of the resulting periodic lattice splits into n_1  subbands each of which gives rise to  Wannier states with a spacing of n_1 lattice sites (red in Fig. <ref>). Since n_1 is large, beyond-nearest-neighbour couplings between these new states are vanishingly small, and thus each subband can be well described by a new tight-binding model with dispersion

    E(κ,ϕ) = E_0(ϕ) - 2J'(ϕ) cosκ,

where J' is the hopping between two neighbouring Wannier states, κ is the renormalised quasimomentum, and E_0 is the mean energy of the subband. In principle, both E_0 and J' depend on the phase ϕ in (<ref>). As we will discuss below, the variations of E_0 are on the order of J', which is exponentially small. Similarly, for large n_1, the variations of J' are exponentially smaller than J' itself, thus they can safely be ignored in any effective theory.


We now determine the dependence of E_0 on ϕ in this periodic approximation by applying the Aubry duality transformation (<ref>). It is important to note that, as β is rational, (<ref>) only generates n_1 distinct reciprocal space modes. In order to make the transformation unitary, (<ref>) is replaced by a Hamiltonian acting on n_1 lattice sites with twisted periodic boundary conditions:

    H = -   J∑_j=1^n_1(a_j^† a_j-1 e^iκ/n_1+ H.c.)
    
        -    Jλ∑_j=1^n_1cos(2π/n_1 (j-ϕ)) a_j^† a_j,

where a_0 = a_n_1. Applying the duality transformation (<ref>) to (<ref>), it becomes

    H = -Jλ/2   [∑_ℓ=1^n_1(b^†_ℓ b_ℓ-1 e^2π iϕ/n_1 + h.c.) 
       +4/λ∑_ℓ=1^n_1cos(2π/n_1(ℓ-κ/2π)) b^†_ℓ b_ℓ].

That is, the duality transformation exchanges the quasimomentum κ and the offset ϕ. By (<ref>), the energy eigenvalue of the dual Hamiltonian depends on its quasimomentum 2πϕ as a simple cosine the amplitude of which is taken independent of κ. Therefore, E_0 also has a cosine dependence on ϕ:

    E_0(ϕ)    = E_0 - (Jλ)' cos(2πϕ)
        
    
        (Jλ)'    = λ/2× 2J'_dual.

Combining (<ref>) and (<ref>), the dispersion relation of the periodic approximation is finally given by

    E(κ,ϕ)    = E_0 -2J'cosκ -(Jλ )' cos (2πϕ).



In the quasiperiodic system however, 1/β differs from n_1 by a small irrational number β_1, therefore, the jth minimum of the potential is shifted ϕ_j=jβ_1 away from a lattice site (see Fig. <ref>). Equation <ref> is thus not exact, but  as β_1 is assumed to be small, ϕ_j changes slowly. Therefore, (<ref>) can be used as an effective Hamiltonian for the new Wannier states of separation 1/β. That is, upon rescaling by 1/β, the resulting model is described by the Hamiltonian

    H =  -2J'cosp̂' -(Jλ)' cos (2πβ_1 x̂'),

an Aubry–André model of parameter β_1 with renormalised potential and hopping terms. The same procedure can then be repeated with step sizes 1/β_k to obtain a renormalisation group treatment of the full spectrum.

The terms J', (Jλ)' entering (<ref>) may be estimated numerically from the scaling of bandwidths over a single step of the procedure. In the limit of large n_k, the scaling of J for states sufficiently far from E=0 can be calculated analytically using the WKB approximation (see Appendix <ref>). These calculations show that the renormalisation of the potential-to-hopping ratio λ does not depend on energy (and hence the place of the subband in the spectrum), and is given by (see Appendix <ref>)

    λ' = 2 (λ/2)^1/β.

Iterating this procedure on the emerging quasiperiodic lattices gives the effective amplitude λ^(k) on length scale ≈ N_k as

    λ^(k) = 2 (λ/2)^1/β…β_k-1≃ 2 (λ/2)^N_k.

For λ<2, λ > λ' > λ” >…: the RG procedure tends to λ^(k)=0, that is, the quasiperiodic modulation becomes irrelevant and hence all eigenstates are extended. On the other hand, if λ>2, λ increases upon renormalisation: the system flows to λ^(k)→∞ where hopping is irrelevant, and all eigenstates are localised. The critical exponent of the reduced tuning parameter g = log(λ/2) is ν=1: indeed, according to (<ref>), |g|=ξ^-1.



















§ CRITICAL SCALING NEAR THE SUPERFLUID–INSULATOR TRANSITION





We performed exact diagonalisation on the single-particle Aubry–André Hamiltonian (<ref>) and extrapolated the behaviour of the truly incommensurate model from the sequence of rational approximations M_k/N_k of β, all implemented with periodic boundary conditions. 


The key quantity we considered was the curvature of the lowest band,

    Γ = 1/2Ja_0^2.^2 ε(k)/ k^2|_k=0 = 1/Jħ^2/m_effa_0^2,

where m_eff is the effective mass of particles near the bottom of the band and a_0 is the lattice spacing. The normalisation is chosen such that Γ for an unmodulated tight-binding chain is unity. In an extended phase, the motion of a single particle becomes ballistic beyond a length scale, therefore, its effective mass tends to a finite value in the limit of an infinite system. Bands of a localised model, however, become completely flat, resulting in an infinite effective mass and thus zero curvature. As a consequence, the limit lim_N→∞Γ_N, where N is the period of the lattice,  can be used as an order parameter in a quantum localisation transition. We approximate the second derivative using the energy difference over a finite segment of the lowest band:

    Γ = 1/JE_Θ-E_0/(Θ/N)^2,

where E_Θ and E_0 are the ground state energies of the system in periodic boundary conditions twisted by Θ or without twist, respectively. Equations <ref> and <ref> are equivalent for Θ→0, but in practice, Γ remains essentially unchanged for significant fractions of π. In this paper, Θ=π/20 was normally used.
In interacting many-particle systems, the appropriate generalisation of Γ gives the superfluid fraction or superfluid stiffness, which is widely used to analyse superfluid–insulator transitions <cit.>.





We note that the curvature of a band is related to its width Δ E  and  therefore can be used to extract the scaling properties of the bandwidth; in a homogeneous or crystalline system, this scaling is governed by the dynamical exponent z:

    Δ E ∼ξ^-z.

To elucidate this connection,  the typical band structure in the extended phase is sketched in Fig. <ref>. 
On length scales N_k≫ξ, the effective potential is irrelevant compared to the effective hopping (that is, the renormalised λ^(k)≪1) and the spectrum of any periodic approximation with N≫ξ becomes similar to the spectrum for N≈ξ: the effective lowest band is folded up, largely conserving the continuity of the spectrum. In particular, the small gaps introduced by the remaining weak potential do not affect the curvature at k=0. That is, regardless of the period N_k≫ξ of the lattice, the lowest dynamical band is ∼π/ξ wide in k-space. Approximating its dispersion by

    ε(k) ∼ -Δ E/2 cos(ξ k),

the band curvature follows as

    Γ∼Δ E/(1/ξ)^2∼ξ^2-z∼ |λ-2|^ν(z-2)

by the definition of z and ν: note that ν=1 for all β in the Aubry–André model (cf. Eq. <ref>). We note that the scaling behaviour of the many-particle superfluid fraction is also given by (<ref>) <cit.>, as expected given its relation to Γ.




 §.§ Results for the Aubry–André model







The curvature of the lowest band was calculated for the Aubry–André model (<ref>) near λ=2 for several different incommensurate ratios and plotted in Fig. <ref>. The rational approximations to β were always chosen  such  that the period of the resulting superlattice was much larger than the longest correlation length considered, ξ_max=10^4. In contrast to homogeneous systems, the order parameter never follows a power law, even when ξ is on the order of thousands of lattice sites. This contradicts the conventional notion of a `scaling regime' where the only relevant length scale is the correlation length, resulting in power law behaviour <cit.>. 

The origin of this discrepancy is the emergence of the arbitrarily large `microscopic' length scales N_k discussed in Secs. <ref> and <ref>. Consider a near-critical Hamiltonian with extended eigenstates of correlation length ξ≈ N_k: broadly speaking, its spectrum displays the first k levels of the hierarchical critical spectrum, but further ones are not resolved and thus have no effect on Γ  (cf. Figs. <ref> and <ref>). As a result, its critical scaling  at ξ≈ N_k depends on the kth step of the renormalisation protocol of Sec. <ref> which is in turn governed by β_k. In particular, the slope of the log-log plot in Fig. <ref> is determined by the local dynamical exponent z_k = z(β_k) defined by

    Δ E(N_k+1)/Δ E(N_k) =(N_k+1/N_k)^-z_k≅β_k^z_k.



In Fig. <ref>(a), the continued fraction expansions of all values of β become periodic with identical periods; this implies that the sequence {β_k} and thus the scaling behaviour is identical from a point on. This is manifest in the identical but shifted curves in the plot; the difference in overall scaling stems from the different initial terms in the continued fraction expansion which result in different N_k's corresponding to the same β_k's. 

In Fig. <ref>(b), the values of β are very close to each other, and so their continued fraction expansions  start with the same terms. Since the first few β_k differ by very little, the critical scaling is almost identical for relatively small ξ: this changes noticeably as further terms in the continued fraction expansions become different, giving rise to completely different scalings. This behaviour demonstrates that while the structure of quasiperiodic systems described by only slightly different incommensurate ratios may be very different on sufficiently long length scales, such differences are immaterial in short samples. Such unpredictability of the large-scale behaviour of quasiperiodic systems also plays a key role in quantum complexity theory <cit.>.


We have thus found that the existence of `microscopic' structure on all length scales prevents the formation of a conventional scaling regime where exact power-law scaling relations such as (<ref>) would hold. For β's with periodic continued fraction expansions, however, the sequence  {β_k} itself is periodic and so the scaling behaviour repeats itself on arbitrarily long length scales.  In this case, one can combine all renormalisation steps in one period into a discrete RG protocol where all steps are identical. For such RG schemes, it is common to find a power-law behaviour on average, with log-periodic oscillations around it <cit.>: indeed, we observe such oscillations in Fig. <ref>(a). Nonetheless, a single period of these oscillations may contain  an arbitrarily complex pattern of the RG steps defined in Sec. <ref>, and so a description of the critical behaviour in terms of log-periodic oscillations is not generally practical.


The continued fraction expansions of almost all irrational numbers are, however, not periodic. For these numbers, the RG protocol cannot be described in terms of a single, if complex, step, resulting in a situation more complicated than the log-periodic oscillations discussed above. In particular, there is no way to sensibly define single critical exponents for the Aubry–André model for these values of β. The critical behaviour is only appropriately described by the detailed dependence of observables on the length scale, an example of which is the set of local dynamical exponents (<ref>). Using the analytic RG procedure discussed in Sec. <ref>, z_k can be calculated for β_k≪1 (see Appendix <ref>). To leading order,

    z_k≈ 1.166 β_k^-1/log(β_k^-1),

meaning that z_k→∞ as β_k→0. Therefore, for an incommensurate ratio β=[0;n_1,n_2,…] with n_k_1<n_k_2 for all K<k_1<k_2 for some K, the conventional definition of the dynamical exponent,

    z = lim_ξ→∞logΔ E/logξ,

diverges: we note that these numbers form a dense, uncountable subset of [0,1]. This marks a completely novel critical behaviour, one not even approximated by power laws.








 §.§ Ground state universality of quasiperiodic models






In addition to the Aubry–André model, we investigated a generalised Hamiltonian that also allows for quasiperiodic modulation of the hopping <cit.>:

    H = - ∑_n   [J + Jμcos(2πβ(n+1/2)+ϕ) ]
       ×( a_n^† a_n+1 + H.c.) 
        -Jλ∑_ncos(2πβ n) a_n^† a_n

where μ is the dimensionless modulation amplitude of the hopping.  Remarkably, (<ref>) still has no mobility edges: localisation transitions occur  simultaneously in all eigenstates, similarly to the simple Aubry–André case <cit.>. The boundary between extended and localised phases is given by

    ∑_±√((λ/2)^2 ±λμcosϕ + μ^2) = 2

for ϕ≠0, regardless of the value of β <cit.>. For ϕ=0, the phase diagram consists of an extended (λ<2;μ<1), a localised (λ>2,2μ), and a critical phase (2μ>2,λ) <cit.>. As examples, localisation transitions along the following paths were considered (see Fig. <ref>):



    ϕ   =0,    μ   =1/2,    λ   ≈2;
        
    ϕ   =π/2,     2μ   =λ=Λ/√(2),    Λ   ≈2.


Even though the hopping in these models is no longer uniform, Γ was calculated  using the unchanged definition (<ref>): it is an appropriate order parameter of the localisation transition regardless of normalisation.

Further to this generalised Aubry–André model, we considered the continuum quasiperiodic Hamiltonian

    H = p̂^2/2m + V_1 cos^2(kx)+V_2cos^2(β kx).

Equation <ref> reproduces the Aubry–André model in the limit V_1≫ E_r≫ V_2 where the recoil energy,

    E_r = ħ^2 k^2/2m,

is the typical kinetic energy scale of the system. In addition to this limit, we studied the case of equal absolute lattice depths V_1=V_2=VE_r/2. Periodic approximations to the Hamiltonian were implemented in momentum space and the curvature of the lowest band was calculated by exact diagonalisation using a formula adapted from (<ref>) <cit.>:

    Γ = .π^2/E_rE_Θ-E_0/(Θ/N)^2|_Θ→0 = m/m_eff.

A localisation transition was observed in the ground state of this model for all tested values of the incommensurate ratio β at a β-dependent critical V_c. Unlike the generalised Aubry–André model, however, its spectrum is unbounded, and several mobility edges appear in the spectrum of excited states. Nevertheless, we expect that the structure of the ground state has a hierarchical structure  similar to that discussed in Sec. <ref>.





To test this hypothesis, the curvature of the lowest band was computed for several rational approximations of β=[0;2] near the transition point of all these models. Since this continued fraction expansion is periodic, effective critical exponents ν and z exist and can be determined using a finite-size scaling method <cit.>. For a homogeneous system near a localisation transition, the finite-size scaling hypothesis can be applied to (<ref>) to give

    Γ = L^z-2Φ(L^1/νδ),

where L is the size of the finite system, δ is the distance from the transition point [e.g., (λ-2) for the Aubry–André model] and Φ(x) is a scaling function determined by the universality class <cit.>. In such systems, all sufficiently large length scales are equivalent: taking Γ(δ) for several different system sizes, critical exponents can be found accurately as the ones resulting in the best collapse of the scaled curves on each other <cit.>. For quasiperiodic models, (<ref>) does not hold in general, but for β's with periodic continued fraction expansions, N_k's separated by a full period of the expansion correspond to the same β_k and thus display the same emergent structure. Using these values of N_k as system sizes or period lengths, (<ref>) applies and fitting to it yields the average dynamical exponent discussed in Sec. <ref>.


The result of such a  fit is shown in Fig. <ref>  for the Aubry–André model, the generalised models (<ref>, <ref>) and the continuum model (<ref>) with β=[0;2]. The resulting critical exponents are the same as are the scaling curves apart from overall rescaling. This suggests strongly that both the generalised Aubry–André transitions and the continuum quasicrystal belong to the same ground state universality class as the Aubry–André model.





For general β, however, the Aubry–André phase transition has no well-defined critical exponents, and so the critical behaviour depends qualitatively on the correlation length. Therefore, such a universality class cannot be described in terms of critical exponents and finite-size scaling functions, only through the detailed dependence of observables on the length scale. To illustrate such universality, the curvature of the lowest band in all models was plotted in Fig. <ref> as a function of the distance from the transition point. The curves can be collapsed on top of each other: points mapped onto each other correspond to an  equivalent correlation length. This notion of universality is markedly different from the conventional one based on the existence of a scaling regime in which the only effect of microscopic structure is to set critical exponents.
















§ MULTIFRACTAL ANALYSIS




The ground state dynamical exponent considered in Sec. <ref> is a key quantity in quantum phase transitions, since at zero temperature, only the behaviour of the ground state is relevant. Unlike most quantum phase transitions, however, localisation transitions in the Aubry–André model and its generalisation (<ref>) occur at the same point for all eigenstates <cit.>, resulting in a fully singular continuous spectrum. A probe of the entire spectrum, as opposed to the ground state only, is also more relevant to experiments on Anderson and many-body localisation.

To explore the overall behaviour of the spectrum, we employed a multifractal scaling technique which yields statistics describing differences in the scaling behaviour at different parts of the spectrum. Furthermore, we demonstrate the connection between the structure of the spectrum and the resulting quantum dynamics by analysing the anomalous diffusion dynamics, a key experimental diagnostic, of the same models at criticality.




 §.§ Formulation




Consider a periodic approximation β=M_k/N_k of the incommensurate Hamiltonian. The singularity strength α_i of the ith subband is defined by

    Δ_i ∼ N_k^-1/α_i,

where Δ_i is the width of the subband; by comparison to (<ref>), the ground state dynamical exponent is 1/α for the lowest subband. For an incommensurate ratio with periodic continued fraction expansion, and hence a uniform scaling behaviour over different length scales, it is expected that the subbands of given singularity strength form a fully self-similar structure, the fractal dimension f(α) of which is given by <cit.>

    Ω(α) ∼⟨Δ⟩^-f(α),

where Ω(α)α is the number of subbands with singularity strength between α and α+α and ⟨Δ⟩=N_k^-1/α is a typical bandwidth of singularity strength α. The function f(α) contains complete information about the scaling behaviour of the spectrum and is routinely used to characterise critical spectra of various systems <cit.>. We note that the Hausdorff dimension of the entire spectrum is the maximum value of f(α) <cit.>.


To accurately find f(α) numerically, we considered the scaling exponents τ_q defined through <cit.>

    ∑_i=1^N_kΔ_i^-τ_q∼ N_k^q.

This set of dimensions gives f(α) through the Legendre transform <cit.>



    α   = τ_q/ q;
    
        f(α)    =  qα - τ_q.

 
It is now straightforward to show (see Appendix <ref>) that the slope of a straight line fit to


    ℵ(N_k;τ) = -∑_i μ_i^(τ)logΔ_i

and

    ϕ(N_k;τ) = -∑_i μ_i^(τ)logμ_i^(τ),

respectively, as a function of log N_k, where

    μ_i^(τ) = Δ_i^-τ/∑_j Δ_j^-τ,


gives α^-1 and f/α corresponding to a particular value of τ; from this, the f(α) curve can be obtained parametrically.








 §.§ Results, universal multifractality






Multifractal analysis using the above formalism was carried out for β=[0;n] (n=2,6,10): the resulting f(α) curves for the Aubry–André model are shown in Fig. <ref>. f(α) is only defined on an interval α_min≤α≤α_max and f(α_min) = f(α_max) = 0: α_min,max give the scaling exponents of the smallest and largest bandwidths of the system, respectively, but these represent a vanishing minority of all bands.  In fact, α_min^-1 equals the ground state dynamical exponent (<ref>) in all cases we considered. This suggests that the narrowest bands of the spectrum are near the bottom (and the top) of it and their scaling behaviour is atypical for the spectrum.


Localisation transitions in the generalised Aubry–André model (<ref>) were also observed to occur simultaneously in all eigenstates <cit.>, giving rise to fully critical spectra at transitions. The multifractal dimensions f(α) at the transition points (<ref>) were thus obtained using the same method. The f(α) curves for the simple and generalised Aubry–André models are identical for a given β: this directly shows that the universality observed in the ground state also applies to the entire spectrum. For β=[0;1], the golden mean, and ϕ=0, this behaviour was already known <cit.>. In this particular case, singular continuous spectra appear away from the localisation transition line as well (cf. Fig. <ref>): in accordance with Ref. Han1994, we found that the multifractal structure of these critical spectra is markedly different from the ones on the transition line (not shown). However, the existence of a critical region appears to be a peculiarity of the ϕ=0 phase diagram <cit.>, thus no universal features are expected of it.


It has been conjectured that the peak of the f(α) curve is at α^* = 1/2 for all β, that is, the Hausdorff measure of the spectrum is dominated by bands scaling as Δ∼ N^-2  <cit.>. While this appears to be  the case for β=[0;2] and maybe for [0;6], it is certainly not  for [0;10] where α^*≈0.515 (the numerical error of α is at most ≈ 0.005). Lower quality evidence for [0;n] with large n suggests α increases further with n: the observation of Ref. Tang1986 appears to be a consequence of only using (more easily accessible) β's with small continued fraction terms. 











 §.§ Expansion of a wave packet







The multifractal dimensions f(α) contain full information on the scaling behaviour of the spectrum, and since the dynamics of a quantum system depends on differences between its energy levels, they capture the dynamical behaviour of the critical system. A straightforward example is the diffusion dynamics of an initially site-localised particle after a sudden quench onto the Aubry–André Hamiltonian (<ref>). This expansion can be characterised through the evolution of the pth moment of the resulting quantum state:

    μ_p    = ⟨ |x-x_0|^p ⟩;    ℓ_p    = μ_p^1/p,

where x_0 is the position where the wave function is initially localised and p is an arbitrary positive real number. In a conventional critical system, ℓ_p∼ t^1/z because t is a characteristic time scale corresponding to the length scale ℓ_p(t) <cit.>. In this context, σ=1/z is commonly referred to as the anomalous diffusion exponent.


Using exact diagonalisation, the time evolution of the initial state can be obtained directly from

    |ψ(t)⟩ = ∑_n |n⟩ e^-iE_nt⟨ n | ψ(0)⟩,

where |n⟩ are the eigenstates of the Hamiltonian with energy E_n: given |ψ(t)⟩, μ_p can be calculated straightforwardly. 
As the details of the expansion dynamics will depend on the choice of initial state <cit.>, we show in Fig. <ref>(a) the evolution of the rms width ℓ_2 averaged over all initial sites x_0 for periodic approximations of β=[0;n] (n=2,6,10) in the Aubry–André model. 
Apart from finite size effects, each expansion follows an approximate power law: fitting a power law to each plot resulted in a diffusion exponent σ≈0.5 within the error of the fit. Similar behaviour has previously been found for other values of β as well <cit.>. On the other hand, σ for a fixed value of β does depend on p, as shown in Fig. <ref>(c) for β=[0;6] and p=1,2,4. This unusual behaviour is readily accessible by measuring higher moments of the diffused density distribution in typical sudden expansion experiments <cit.>.


In addition to the Aubry–André model, ℓ_p(t) was calculated by the same method for the critical point λ=2, μ=1/2, ϕ=0 of the generalised Aubry–André Hamiltonian: ℓ_2 for β=[0;10] was plotted in the inset of Fig. <ref>(a) together with ℓ_2 for the simple Aubry–André model. The exponents of the approximate power laws were found to match, together with the structure of oscillations around it:

    ℓ_2^GAA(t) = ℓ_2^AA(0.867t)

holds accurately for all but the shortest time scales.








 §.§ Connection between expansion dynamics and spectrum multifractality


In order to connect the expansion dynamics in a critical tight-binding model to the multifractal properties of the spectrum, consider the Aubry–André model with an arbitrary value of β with periodic continued fraction expansion. Since the only natural length and time scales of the problem are the lattice spacing a_0 and the `hopping time' ħ/J, μ_p depends on these scales as

    μ_p(t;x_0,β,λ) = a_0^p m_p(Jt;x_0,β,λ),

where m_p is now a dimensionless function of dimensionless variables. To get an overall description of the critical dynamics, we set λ=2 and average over the position of the initial site:

    m_p(Jt;β) = lim_N→∞1/2N+1∑_x_0=-N^N m_p(Jt;x_0,β,λ=2).

Consider now the kth step of the renormalisation process outlined in Sec. <ref>: the spectrum consists of N_k critical subbands with incommensurate ratio β_k; let the effective hopping term in each be J_i (1≤ i≤ N_k). Provided the time t is longer than the time scales corresponding to typical band gaps, interference between bands averages out, leaving

    m_p(Jt;x_0,β) ≃∑_i=1^N_k|⟨ x_0^(i)| x_0⟩|^2 N_k^p m_p(J_it;x_0^(i),β_k),

where |x_0^(i)⟩ is the Wannier state of the ith subband living (among others) on site x_0;  the factor N_k^p is due to the renormalisation of the lattice spacing. To average (<ref>) over lattice sites, we note that each renormalised band has one Wannier state per N_k lattice sites and the sum of the overlap integrals |⟨ x_0^(i) | x_0⟩|^2 over all x_0 is 1 since the |x_0⟩ form a basis. As a result, the overlap integrals average to 1/N_k for all lattice sites, and hence

    m_p (Jt;β) ≃ N_k^p-1∑_i=1^N_km_p(J_it;β_k).



Now, consider those k that correspond to full periods of the continued fraction expansion, that is, β_k=β. Assuming that the expansion is governed by a power law at long times,

    m_p(Jt;β)∝ (Jt)^pσ_p   (Jt→∞),

Eq. <ref> gives

    N_k^p-1∑_i=1^N_k (J_it)^pσ_p   ≃ (Jt)^pσ_p
    ∑_i=1^N_k J_i^pσ_p   ∝∑_i=1^N_kΔ_i^pσ_p∝ N_k^1-p

where Δ_i is the width of the ith subband for β=M_k/N_k, 4J_i in the unmodulated tight-binding approximation. In terms of the multifractal dimensions introduced in Sec. <ref>, the anomalous diffusion exponents σ_p are given by

    σ_p =-τ_1-p/p.
 
In contrast to conventional diffusion dynamics, σ_p now depends on p and is not equal to the inverse of the ground state critical exponent. The only crucial assumption in deriving (<ref>) is the self-similarity of the spectrum, therefore, we expect it to hold for the dynamics of other singular continuous spectra, e.g., the  Fibonacci quasicrystal <cit.>. In particular, as the spectra of all generalised Aubry–André transition points are described by the same multifractal exponents, the σ_p are universal too. The differences seen at very short times can be attributed to initial renormalisation steps required to attain a fixed point. 


An interesting special case is that of p=2. There is strong numerical and analytical evidence <cit.> suggesting that for the Aubry–André Hamiltonian with rational β=M/N, the sum of bandwidths scales as

    lim_N→∞ N∑_i=1^N Δ_i ≈ 9.3299

regardless of M. This implies that τ_-1=-1 for any β: comparing with (<ref>), we find that σ_2 = 1/2, as seen numerically in Fig. <ref>(a). Unlike diffusive systems, however, σ=1/2 here cannot be regarded as the consequence of a random walk between scatterers since σ_p≠1/2 in general.



In Fig. <ref>(b), we note that oscillations around the approximate power law scaling of ℓ_p(t) decrease with time and become unnoticeable for sufficiently long times. The origin of this behaviour is clear from (<ref>): for Jt≫1, the expansion dynamics can be regarded as a superposition of the same dynamics at earlier times J_it. Since these J_i range over several orders of magnitude for sufficiently large N_k, m(J_it) probes any short-time oscillations over several periods, thus averaging them out. That is, expansion length scales in different subbands can be very different, of which ℓ_p is only an average. This distinction becomes manifest in the expansion dynamics for β's with aperiodic continued fraction expansions: while expansion dynamics at different length scales is different, at any particular time, these are averaged out, preventing the formation of clean crossovers similar to those seen in Fig. <ref> for a single sequence of subbands (namely, the ground state).



















§ CONCLUSION




We have investigated the critical behaviour of the Aubry–André model and other one-dimensional quasiperiodic systems near their localisation transitions. 
In particular, we considered the dependence of energy scales near the ground state, Δ E, on the correlation length ξ. While the standard theory of phase transitions dictates that for large ξ, the system attains a scaling regime in which Δ E∝ξ^-z, we found that the critical behaviour is not described accurately by a power law on arbitrarily large length scales. 


This is caused by the hierarchical structure of the critical spectrum of quasiperiodic models, captured by the continued fraction expansion of the irrational number β=[0;n_1,n_2,…] describing their incommensurability. Each continued fraction term n_k has associated with it a length scale N_k: scaling properties of the critical spectrum near this length scale were found to be fully determined by n_k.
Since the spectrum of a system near a phase transition  is sensitive to  spatial features on length scales up to the correlation length ξ, the critical behaviour of quasiperiodic models at ξ≈ N_k will also be governed by n_k. As the sequence of these n_k can be arbitrary and is controlled by the precise value of β, the dynamical exponent z can typically not be defined for quasiperiodic models.  As an example, we found that for a wide class of β's, Δ E tends to zero faster than any power of ξ, heralding critical behaviour qualitatively different from any conventional system.
Furthermore, the dependence of the critical behaviour on the incommensurate ratio is unusual: arbitrarily close values of β can result in qualitatively different asymptotic behaviours very near the transition, as their continued fraction expansions eventually start to deviate.

Even though the localisation transition of one-dimensional quasiperiodic models cannot be described by power law relations, we find numerically that transitions in different models sharing the same value of β display universal features. Instead of critical exponents, such universality classes are described  by the detailed dependence of observables such as Δ E on the correlation length. For models belonging to the same universality class, such functions can be scaled onto each other, similarly to finite-size scaling techniques for conventional phase transitions. The origin of such universality remains the identical behaviour under the renormalisation of length scales; the key difference is that quasiperiodic systems only admit a single sequence of discrete renormalisation steps that themselves depend on the length scale.


To complement studies of the ground state, we considered scaling properties of the entire spectrum on different length scales. For β's with a periodic continued fraction expansion, the spectrum is expected to be self-similar at the Aubry–André critical point: its structure was found to be a multifractal, and multifractal dimensions were calculated for several values of β. We also investigated the expansion dynamics of a localised wave packet and found that the evolution of the spread ⟨ r^p⟩^1/p of the wave function is described by a power law the exponent of which depends on p and β. This is at odds with the behaviour of diffusive systems, where this exponent is 1/2 for all p. Similarly to ground state properties, we again found universality between transition points of different quasiperiodic models in both their multifractal spectrum and expansion dynamics.


For the Aubry–André model, we used a discrete renormalisation group protocol <cit.> to construct the critical spectrum and thus explicitly calculate the scaling of Δ E with ξ; non-power-law universality classes could be understood through the renormalisation behaviour of  other types of quasiperiodic models near phase transitions. 

Quasiperiodicity in higher dimensions leads to the emergence of arbitrarily large `microscopic' length scales the same way as in one dimension: this discrete large-scale structure is manifest in sharp diffraction peaks at progressively smaller momenta <cit.>. Therefore, it is reasonable to expect that phase transitions in such systems (including material quasicrystals) also display non-power-law behaviour. 
In general, quasiperiodic systems open the door to more complex large-scale behaviours, especially with interactions, which can show up, for instance, in increased quantum complexities <cit.>, as novel universality classes for the many-body localisation transition  <cit.>, and in conjunction with their inherited topological features <cit.>.
















§ ACKNOWLEDGEMENTS

We are grateful to Ehud Altman, Bartholomew Andrews, David Huse, and Austen Lamacraft for stimulating discussions and insights. This work was partly funded by the European Commision ERC starting grant QUASICRYSTAL and the EPSRC Programme Grant DesOEQ (EP/P009565/1). 













§ WKB THEORY OF TIGHT-BINDING MODELS




In this appendix, we develop a semiclassical theory of tight-binding lattices with potentials slowly varying compared to the lattice spacing. The derivations presented here follow closely the standard derivations of WKB theory for an ordinary, quadratic dispersion relation <cit.>. Since the period of the incommensurate modulation, β^-1 is large, this theory is applicable to the Aubry–André model for the class of β's considered, and can be used to accurately estimate the renormalised hopping and thus the critical exponents ν and z <cit.>.








 §.§ Construction of the wave function



We assume that the period of the modulating potential is very much larger than the lattice spacing. In this case, the discreteness of the wave function becomes irrelevant, and the Hamiltonian can be written as (the unit of length is the lattice spacing, ħ=1)

    H = -2Jcosp̂ + V(x̂).

where the nonquadratic dependence on p̂ follows from the tight-binding dispersion relation. Due to this nonquadratic dispersion relation, the quasiclassical wave numbers depend differently on energy:

    k(x)   = arccos(V(x)-E/2J) ;
        
    κ(x) = ik(x)    = arcosh| E-V(x)/2J|,

Using k(x), the Schrödinger's equation (<ref>) and the WKB ansatz can be written as

    0    = -(cosp̂)ψ(x) + cos k(x) ψ(x)
        
    ψ(x)    ≈ A(x)  exp(± i ∫^x k(x') x') = A(x)ϕ(x).

where both A(x) and k(x) are assumed to vary slowly. Due to this slow variation, considering terms with a different number of derivatives amounts to separation of scales: in first order WKB approximation, only terms with zero or one derivatives are retained. The nth derivative of ψ(x) is given by

    ψ^(n)(x)    = Aϕ^(n) + nA'ϕ^(n-1) + O(A”)  where
    ϕ'(x)    = ± ik ϕ
    ϕ”(x)    = (± ik)^2 ϕ± ik'ϕ
    ϕ”'(x)    = (± ik)^3ϕ + 3(± i)^2 kk'ϕ± ik”ϕ
    ∴ϕ^(n)(x)    = (± ik)^n ϕ + n2 (± i)^n-1 k^n-2k'ϕ +O(k”);

Eq. <ref> can be proved by induction. Combining (<ref>) and (<ref>) gives ψ^(n) and p̂^n ψ as 

    ψ^(n)   = (± ik)^n Aϕ + n2 (± ik)^n-2 (± ik')Aϕ
        +n (± ik)^n-1 A'ϕ;
    p̂^n ψ   = (± k)^n Aϕ - i (±1)^n-1n2 k^n-2 k'Aϕ
        -i n (± k)^n-1 A'ϕ.

Writing cosp̂ as a Taylor series, we finally obtain

    (cosp̂)ψ   = ∑_n=0^∞ (-1)^n p̂^2nψ/(2n)!
         = ∑_n=0^∞{ (-1)^n k^2n/(2n)!ψ∓
          i [(-1)^n/2k^2n-2/(2n-2)!k'Aϕ + (-1)^n k^2n-1/(2n-1)!A'ϕ] }
       = cos k×ψ± i (12cos (k) k'A -sin (k) A')ϕ.

Writing this into (<ref>) yields

    1/2cos(k) k'A -sin (k) A'    = 0
    
        A(x)    ∝1/√(sin k(x)).

Noting that the velocity of a classical particle moving under this Hamiltonian would be

    v = ẋ = ∂ H/∂ p = 2Jsin p,

A(x) can be interpreted as reproducing the classical probability of the particle being found at x, similarly to the amplitude in standard WKB theory <cit.>.


The derivation above does not depend on k(x) being real. At points with too large potentials, k = iκ with κ defined in (<ref>), and the wave function (<ref>) becomes

    ψ(x) ∝1/√(sinhκ(x))exp(-∫_x_0^x κ(x') x'),

where the classical turning point x_0 is given by E = V(x_0)-2. At this turning point, k=0, and so the cosine dispersion may be replaced with a quadratic one: as a result, the Schrödinger's equation near the turning point reduces to the Airy equation. Solving this equation gives connection formulae equivalent to those in standard WKB theory:

    C/2√(sinhκ)   exp(-∫_x_0^x κ x') 
       ⟷C/√(sin k)cos(∫_x_0^x k x -π/4).

The similarity of the connection formulae to standard WKB also means that the Bohr–Sommerfeld quantisation condition holds for this dispersion relation too:

    ∫_x_0^x_1 k(x) x = (n+1/2)π.

We note that the region E>V(x)+2 is also inaccessible classically. There, k(x)=iκ(x)+π, corresponding to an exponentially decaying wave function changing signs at every lattice site. Eqs. <ref>, <ref>, and <ref> generalise straightforwardly; we shall not discuss them in detail as they only become relevant near the top of the Aubry–André spectrum.


Finally, we find the normalisation constant C for a wave function living in a single potential minimum. Ignoring the exponentially decaying part, the normalisation requirement is

    1    = ∫_x_0^x_1C^2/sin pcos^2 ϕ(x)  x ≈∫_x_0^x_1C^2/sin p1/2
       = C^2/2∫_x_0^x_12J/v x = JC^2  ∫_x=x_0^x_1 t = JC^2 T/2= C^2 Jπ/ω


    C = √(ω/Jπ) ,

where ω is the frequency of classical oscillations in the well.









 §.§ Hopping between neighbouring wells



Consider a potential consisting of identical, centrosymmetric wells centred on x=na,  n∈ℤ. If a is large compared to the classically allowed region near the minimum of the potential, there is only appreciable hopping between neighbouring minima, and its value can accurately be estimated using WKB approximation. This calculation follows that of Ref. landau9 (55, Problem 3) which solves the same problem for a quadratic dispersion.


Assuming that the overlap between wave functions Ψ(x-an) living in neighbouring wells is small, each one can be treated as a Wannier function, that is, Bloch states are of the form

    ψ_k(x) = C∑_n=-∞^∞ e^ikanΨ(x-an).

The Schrödinger's equation for a single well and for the Bloch state are then 


    -Jψ_k(x-1) -Jψ_k(x+1) + [V(x)-ε_k]ψ_k(x)    = 0
        
    
        -JΨ(x-1) -JΨ(x+1) + [V(x)-ε]Ψ(x)    = 0,


where ε is the energy of a well state in isolation and ε_k is the dispersion of the resulting band. Multiplying (<ref>) by Ψ(x), (<ref>) by ψ_k(x), subtracting and integrating from x=-a/2 to a/2 gives


    (ε_k-ε)ψ_k(x)Ψ(x)    +J[ψ_k(x-1)Ψ(x) +ψ_k(x+1)Ψ(x) -ψ_k(x)Ψ(x-1) -ψ_k(x)Ψ(x+1)] = 0
     (ε_k-ε) C    = -J [ (∫_a/2-1^a/2 - ∫_-a/2-1^-a/2) Ψ(x)ψ_k(x+1) x + (∫_-a/2^-a/2+1-∫_a/2^a/2+1)Ψ(x)ψ_k(x-1) x]
       = -J[∫_α-1/2^α+1/2{Ψ(x-1/2)ψ_k(x+1/2) - Ψ(x+1/2)ψ_k(x-1/2)}]^a/2_α=-a/2.


Consider the integral in brackets. At α=a/2, only the n=0 and n=1 terms are relevant in (<ref>):

    ψ_k(x) = C(Ψ(x)±Ψ(a-x)e^ika)

where the sign depends on whether the well eigenstate in question is even or odd. Substituting this form in the integral of (<ref>) gives

    ∫_(a-1)/2^(a+1)/2⋯  x    = C ∫_(a-1)/2^(a+1)/2[Ψ(x-1/2)Ψ(x+1/2) ±Ψ(x-1/2)Ψ(a-1/2-x)e^ika.
       . -Ψ(x+1/2)Ψ(x-1/2) ∓Ψ(x+1/2)Ψ(a+1/2-x)e^ika] x
       = ∓ Ce^ika∫_-1/2^1/2[Ψ(a+1/2+x)Ψ(a+1/2-x) -Ψ(a-1/2+x)Ψ(a-1/2-x) ] x.

Similarly, for α=-a/2, the n=0 and n=-1 terms yield

    ∫_-(a+1)/2^-(a-1)/2⋯  x = ± Ce^-ika∫_-1/2^1/2[Ψ(a+1/2+x)Ψ(a+1/2-x) -Ψ(a-1/2+x)Ψ(a-1/2-x) ] x,

and hence by (<ref>), 

    ε_k-ε = ± 2Jcos ka ∫_-1/2^1/2[Ψ(a+1/2+x)Ψ(a+1/2-x) -Ψ(a-1/2+x)Ψ(a-1/2-x) ] x.

To evaluate each term of the integral in (<ref>), we employ a saddle point approximation to (<ref>): writing

    lnsinhκ(α+x) = lnsinhκ(α) + rx + sx^2 + O(x^3),

we obtain

    ∫_-1/2^1/2Ψ(α+x)Ψ(α-x) x    ≈Ψ(α)^2 ∫_-1/2^1/2exp(-κ x -κ'/2x^2 +r/2 x+s/2x^2) exp(+κ x -κ'/2x^2 -r/2 x+s/2x^2) x 
       =Ψ(α)^2 ∫_-1/2^1/2 e^-(κ'-s)x^2 x.

The Gaussian integral is only significantly different from 1 if the factor multiplying x^2 is O(1); however, under the WKB approximation, κ changes very slowly and so κ', s∝κ”≪ 1. That is, (<ref>) can be written as 

    ε_k-ε   ≈± 2Jcos ka [ Ψ^2(a+1/2) - Ψ^2(a-1/2)]
        ≈± 2Jcos ka·Ψ^2(a/2) ( e^-κ - e^κ) 
        = ∓ 4Jcos ka ·Ψ^2(a/2) sinhκ,

where only the dominant variation in Ψ(x) due to exponential decay was retained. Finally, substituting the wave function (<ref>, <ref>, <ref>) yields

    ε_k-ε   ≈∓ω/πexp(-2∫_x_0^a/2κ(x) x)cos ka.

That is, each well eigenstate broadens into a tight-binding type band with effective hopping term

    J' = ±ω/2πexp(-2∫_x_0^a/2κ(x) x).




















§ RENORMALISATION OF Λ




From (<ref>), the renormalised hopping is given by

    λ' = λJ'(4/λ,2E_0/λ;β)/J'(λ,E_0;β)

where J'(λ,E;β) is the hopping term (<ref>) for a band at energy E in an Aubry–André model with parameters λ and β. Substituting (<ref>) gives

    λ' =  λT/T'exp(2∫_x_0^1/(2β)κ(x) x - 2∫^1/(2β)_x_0'κ'(x) x)

where T is the classical period of oscillation around the minimum and κ(x) is the imaginary wave vector (<ref>); primes denote quantities of the dual model. For brevity, we write ε = -E/J.








 §.§ Relation of T and T'



From the classical velocity–momentum relation (<ref>),

    T    = ∮ x/ẋ = 4∫_0^x_0 x/2Jsin k(x)
       = 4/J∫_0^x_0 x/√(4 - (ε-λcos(2πβ x))^2)
       = 4/2Jπβ∫_0^α_0α/√(4-(ε-λcosα)^2)
        = 4/2Jπβ∫_ε-2^λ y/√(λ^2-y^2)√(4-(ε-y)^2),

where α=2πβ x is the phase of the modulating potential and y = λcosα. Very similarly, the classical period of the dual is 

    T'    = 4/2Jπβ∫_0^α_0'α/√(4-(2/λε -4/λcosα)^2)
       =2λ/2Jπβ∫_ε-λ^2  y/√(4-y^2)√(λ^2-(ε-y)^2),

where now y = 2cosα. The two integrals can be turned into each other by changing y into y'=ε-y, therefore they are equal: (<ref>) becomes

    λ' = 2  exp(2∫_x_0^1/2βκ(x) x - 2∫^1/2β_x_0'κ'(x) x).










 §.§ Evaluating the integrals ∫κ  x






To evaluate (<ref>), we first rewrite the integrals in terms of the phase α=2πβ x:

    λ'    = 2exp(2/2πβ∫_α_0^πκ(α)α - 2/2πβ∫_α_0'^πκ'(α)α) 
       = 2exp(I-I'/πβ),

where I,I' are integrals independent of β, defined as



    I    = ∫_α_0^πarcosh(ε-λcosα/2) α   α_0    = arccos(ε-2/λ)
        
    
        I'   = ∫_α_0'^πarcosh(ε-2cosα/λ) α   α_0'     =arccos(ε-λ/2).


These integrals can be thought of as the area in (α,κ) space bounded by 2coshκ+λcosα=ε and λcoshκ+2cosα=ε, respectively. Introducing the variables a = cosα, k=coshκ, the area integrals can be rewritten as

    I,I' = ∬ a  k/√((1-a^2)(k^2-1));

the integration areas are bounded by the lines a=-1, k=1, and 2k+λ a = ε (for I) or λ k + 2a =ε (for I'; see Fig. <ref>). It follows that I-I' entering (<ref>) is the integral of the same integrand over the difference of the two domains. Since for λ<2, 

    A' = ε-λ/2   > ε-2/λ= A;    K' = ε+2/λ   > ε+λ/2 = K

(and vice versa for λ < 2), this area difference is a quadrilateral bounded by all four lines bounding the triangles (see Fig. <ref>).


For simplicity, we assume that λ is infinitesimally close to 2: λ=2+η, |η|≪1. In this case, the difference quadrilateral is infinitesimally thin: slicing it along lines of constant a gives the integral

    I-I' ≃∫_-1^A_0 a [k(a)-k'(a)]/√(1-a^2)√(k_0(a)^2-1).

In writing (<ref>), we have ignored the variation of k across one slice in the denominator, and replaced it with k_0(a) corresponding to λ=2: this introduces first order corrections to the denominator which, since k-k' is first order in η, can be ignored. Now,

    k(a)-k(a') = ε-λ a/2 - ε -2a/λ≃(ε/4-a)η


    ∴ I-I'   ≃η∫_-1^(ε-2)/2 a(ε/4-a)/√(1-a^2)√((ε/2-a)^2-1)
       = η∫_-e+1^e+1x  x/√((1-(e-x)^2)((e+x)^2-1))
       = η/2∫_(1-e)^2^(1+e)^2(x^2)/√(4e^2-(x^2-(1+e^2))^2)
    
        I-I'   ≃πη/2.

Writing this into (<ref>) gives

    λ' ≃ 2e^η/(2β)≃ 2+η/β;

that is, the reduced tuning parameter λ-2 increases by a factor of β^-1 on rescaling. 


It is possible to evaluate I-I' for an arbitrary value of λ. We omit the derivation due to its length and report that

    I-I' = πlogλ/2λ' = 2(λ/2)^1/β

as stated in Sec. <ref>.

















§ RENORMALISATION OF THE HOPPING, THE DYNAMICAL EXPONENT







To estimate the ground state dynamical exponent corresponding to a particular length scale, we consider the definition (<ref>),

    z_k = log(Δ E(N_k+1)/Δ E(N_k))/logβ_k≈log(J_k+1/J_k)/logβ_k,

at λ=2; since the renormalisation of J in one RG step only depends on β in that step, we anticipate that z_k only depends on β_k.
We first consider the Bohr–Sommerfeld quantisation condition (<ref>) for β≪1: in terms of the phase α=2πβ x,

    ∮ k(α)α = 2(2n+1)π^2β,

where the integrand is given by 2cos k +2cosα = -E/J =ε. For small values of k and α, both cosines can be approximated as quadratics: the contour of the area integral becomes approximately a circle, and thus

    ∮ k(α)α   ≃πα_0^2 ≃π(4-ε)
    ε_n    ≃ 4-2(2n+1)πβ.

In particular, n=0 in the ground state, and so ε_0 ≃ 4-2πβ. The most important consequence of this is that the ground state energy in the limit β≪1 is close to -4J and thus most of the distance between two neighbouring minima is classically unaccessible.


Consider now the expression (<ref>) of the renormalised hopping. By the quadratic approximation introduced above, the classical motion around a minimum can be treated as harmonic; the frequency follows from the coefficients of p^2 and x^2 as ω≃ 4Jπβ. Similarly to (<ref>), J' can now be written as

    J'    ≃ 2Jβexp(-I/πβ)
    ∴ z(β)    ≃ 1 + I/πβ|logβ|

where I is given by (<ref>); it also depends on β through the ground state energy. Since ε≈4 for any small β, the leading order term in z(β) can be obtained by assuming ε=4 and thus α_0=0:

    I    ≃∫_0^πarcosh(2-cosα)α≈ 3.6639
    
        z(β)    ≃ 1.1662 β^-1/log(β^-1).

That is, the ground state dynamical exponent diverges as β_k→0, as discussed in Sec. <ref>. More accurate estimates can be obtained by numerically solving (<ref>) for ε and evaluating (<ref>) directly. 


To provide a numerical check on this result, the ground state dynamical exponent was obtained by the finite-size scaling method outlined in Sec. <ref> for β_n=[0;n], 2≤ n≤ 10. For these numbers, β_k=β for all k, and so the average dynamical exponent yielded by the finite-size scaling procedure equals z(β). In addition, z(1/n) was estimated by calculating the lowest bandwidth for β=1/n and equating it to 4J' in the first and only step of the RG procedure. The resulting critical exponents are plotted against β^-1 in Fig. <ref> together with the z(β) curve predicted by WKB theory. The correspondence between numerical and analytic results improves with decreasing β_n, as expected from the underlying assumptions of the analytic theory.















§ NUMERICAL COMPUTATION OF F(Α)




Due to its definition (<ref>), it is more straightforward to obtain q for a given value of τ than the other way around. Therefore, we consider the alternative Legendre transform



    α^-1   =  q/τ;
    
        f' = f(α)/α   = q - τα^-1.


In principle, this Legendre transform  could now be obtained from q(τ), given by power law fitting to (<ref>), numerically: however, taking derivatives numerically tends to introduce significant noise. To mitigate this, we perform the Legendre transform before power law fitting, as suggested by Ref. Chhabra1989. Equation <ref> is equivalent to

    q = lim_N_k→∞log(∑_i Δ_i^-τ)/log N_k;

writing this into (<ref>) gives

    α^-1(τ)    = /τlim_N_k→∞log(∑_i Δ_i^-τ)/log N_k
       = lim_N_k→∞1/log N_k/τlog(∑_i Δ_i^-τ)
       = -lim_N_k→∞1/log N_k∑_i Δ_i^-τlogΔ_i/∑_j Δ_j^-τ;
    
        f'(τ)    = q - τα^-1
       = - lim_N_k→∞1/log N_k[∑_i Δ_i^-τlogΔ_i^-τ/∑_j Δ_j^-τ - log(∑_i Δ_i^-τ)]
       = - lim_N_k→∞1/log N_k∑_i Δ_i^-τ/∑_j Δ_j^-τlogΔ_i^-τ/∑_j Δ_j^-τ.

That is, α^-1(τ) and f'(τ) are given by fitting a straight line to


    ℵ(N_k;τ)    = -∑_i μ_i^(τ)logΔ_i,
    ϕ(N_k;τ)    = -∑_i μ_i^(τ)logμ_i^(τ),

respectively as a function of log N_k, where

    μ_i^(τ) = Δ_i^-τ/∑_j Δ_j^-τ,


as stated in Sec. <ref>.















