Universality and quantum criticality in quasiperiodic spin chains

Quasiperiodic systems are aperiodic but deterministic, so their critical behavior differs from that of clean systems and disordered ones as well. Quasiperiodic criticality was previously understood only in the special limit where the couplings follow discrete quasiperiodic sequences. Here we consider generic quasiperiodic modulations; we find, remarkably, that for a wide class of spin chains, generic quasiperiodic modulations flow to discrete sequences under a real-space renormalization-group transformation. These discrete sequences are therefore fixed points of a functional renormalization group. This observation allows for an asymptotically exact treatment of the critical points. We use this approach to analyze the quasiperiodic Heisenberg, Ising, and Potts spin chains, as well as a phenomenological model for the quasiperiodic many-body localization transition.


Supplementary Note 1.
This section elaborates on various aspects of the ground-state physics of quasiperiodic spin chains. The key result of this section is an explanation of why generic quasiperiodic patterns flow to sequences; we also discuss other aspects of quantum criticality in these spin chains.
A. Pattern of the local minima in the initial potential We noted in the main text that the local minima of the initial potential 1 + cos(2πϕi + θ)-i.e., sites at which is smaller than at either of the neighbors-follow a sequence structure that sharpens into a discrete sequence under renormalization. We here provide a proof of this statement for this specific potential. It is helpful, for this section, to define the fractional part of a real number x, denoted {x}, as the difference between x and the largest integer less than x, i.e., {x} ≡ x − x .
Thus we have the following chain of results: letter 'B' is at n th position in the Fibonacci word sequence shifted from the standard Fibonacci sequence by distance n 0 such that proving the claim that local minima of the initial potential indeed follow the pattern of the letter 'B' in the Fibonacci sequence.
Thanks to the above result, we can decompose the cosine coupling into A-couplings and B-couplings arranged in a Fibonacci sequence. To emphasize this decomposition, we shall relabel the couplings n as follows: Note that A 0 (n) and B 0 (n) are not defined for all values of n: for a given value of n, either A 0 (n) or B 0 (n) is defined.

B. Sharpening of the sequence structure under RG
For the XXX spin chain, the RG rules are given by with c = ln 2 where n is the smallest . We remark that the whole Fibonacci sequence is built up of word patterns 'ABABA' and 'ABA'. This means that we can apply the decimate all the local minima B in one go, arriving at the following renormalized couplings: where m labels the number of such Fibonacci RG steps. In order to get A m+1 we decimate two B m 's, and to get B m+1 we decimate one B m . If A m and B m form a Fibonacci sequence then A m+1 and B m+1 also follow the Fibonacci pattern. This follows from the inflation rule: A → ABABA and B → ABA. Thus starting from a Fibonacci word sequence, if we replace 'ABABA' by 'A' and 'ABA' by 'B', we again get a new Fibonacci sequence.
Note: We are considering periodic boundary conditions, so our initial system size should be an even number. Also we will consider system sizes given by Fibonacci numbers, F n , as the fractional part of F n ϕ is ϕ −n , i.e. F n ϕ ≈ F n+1 . Since all even Fibonacci numbers can be written as F 3l , we consider initial system sizes of the form N = F 3l . In a Fibonacci sequence of size F 3l , we have F 3l−2 B letters and F 3l−1 A letters. Since this is a word sequence we can decimate all the F 3l−2 'B' couplings in one go. Each decimation decreases the number of couplings (bonds) by 2. So after F 3l−2 decimations, the number of remaining bonds is F 3l − 2F 3l−2 = F 3(l−1) . This corresponds to what we call a Fibonacci step.
(1 + cos(1 + 2πϕ(n + i))) − (1 + cos(2πϕ(n)) + c, = 1 + c − cos(2πϕn)(1 − 2cos(2πϕ)), = 1 + c + λ 2 cos(2πϕn), where (1 + 2(cos(4πϕ) − cos(2πϕ)) ≡ λ 1 ≈ 2.64 and −(1 − 2cos(2πϕ)) ≡ λ 2 ≈ −2.47. The label n in A 1 (n) and B 1 (n) only takes certain integer values, corresponding respectively to the positions-on the original lattice-of the central bonds on the ABABA and ABA clusters. Thus for instance ( Fig. 1) if A 1 (0) is defined, at the center of the A cluster, then B 1 (n) can only be defined for n = 4. One can straightforwardly check that at this level all B 1 (n) that are defined are local minima, i.e., they are smaller than the neighboring A 1 couplings. By considering the range of possible values of n [i.e., the on-site phase] for which an ABABA sequence is possible, we find the following constraint: Thus the bandwidth for A-type couplings on the new, decimated chain is ≈ 0.1. A similar analysis for the B couplings-now considering patterns of the form AABAA-yields the constraint cos(2πϕn) ∈ (−1, −1 + δB 1 ) where δB 1 = 1 − cos 2π 5ϕ−8 2 (≈ 0.04). See Fig. 1. To summarize, after one step of decimation, all the couplings within each "type" (A or B) become similar in value, but are separated from couplings of the other type by ≈ c. This leads to the formation of well-separated bands, which then become increasingly well-separated under the action of the RG.
To form A and B couplings in the subsequent Fibonacci steps, the restriction on the initial value of the middle couplings becomes tighter. For example, to get an A 2 coupling we need to have 'ABABAABAABABAABAABABA' Supplementary Figure 1. Part of the initial coupling distribution. The green 'W' shape is the pattern 'ABABA' in the Fibonacci sequence; this will give rise to A block in the next generation of Fibonacci sequence. The orange 'V' shape is the pattern 'ABA' which will lead to a B block in the next generation. Any block with length in the region above the brown dashed line will be the middle 'A' block in the pattern 'ABABA', while blocks in the region between the red and black dashed line will constitute the 'B' blocks in 'ABABA'. Blocks with length below the blue dashed line will be a 'B' block in the pattern 'ABA'.
patterns in the initial potential, but to have such a long pattern the middle 'A' coupling needs to be closer to cos 2πϕn = 1 than it was required to form A 1 . Thus we expect the fluctuations δA m and δB m to go to zero as m → ∞. Empirically these fluctuations decay exponentially with m. This can be iterated to get (we are using C(n) ≡ cos 2πϕn and S(n) ≡ sin 2πϕn for brevity), Similarly we have For large m, A m (n) ≈ 1 + m(m + 1)c + 1 cos πϕ and B m (n) ≈ 1 + m 2 c + 1 cos πϕ . A m − B m m→∞ − −−− → mc. If c = 0 (corresponding physically to Ising and XX chains that can be mapped onto free fermions), then the sequence structure does not survive under RG.

C. Generic potentials
We showed above that for couplings of the form − log(J n ) = a + cos(2πϕn + θ), the initial distribution flows to a sequence under RG. But this result is quite general and almost all initial distribution of the couplings flow to a sequence. If f is any monotonic function in the range (−1, 1), is bounded and n = − log J n = f (cos(2πϕn + θ)), then the minima can be identified exactly as in the analysis above. Therefore, under RG, the couplings flow to the Fibonacci sequence. Natural functions like logarithm, exponential, etc. all satisfy the above criteria. For example, couplings of the form n = − log(1+ +cos(2πϕn+θ)) (corresponding to the natural choice J n = 1+ +cos(2πϕn+θ)) will flow to a sequence -and we verified this numerically -where is a positive constant which prevents the couplings from hitting the singularity. We have also checked that adding higher harmonics to J n such as cos(4πϕn + δ) with a different phase δ does not affect the flow to sequences.
Note that f need not be a monotonic function. E.g for f (x) = |0.5 + x| and f (x) = sin(2x), we get Fibonacci sequences under RG flow even though the above functions are non-monotonic in the range (-1,1). But if the function f has too many extrema, we find that the Fibonacci sequence structure gets destroyed, e.g f (x) = sin(4x).

D. Correlation length exponent
The fixed points of the above RG is unstable against dimerization of even and odd couplings. To find the critical exponent ν associated with the fixed point, we study the behavior of the sequence under an asymmetric perturbation. We dimerize the perfect sequence as Under RG, after m Fibonacci steps the sequence flows to, The asymmetry between even and odd blocks keeps increasing with m and eventually will become positive, destroying the Fibonacci pattern and driving the system to a phase.
In a single Fibonacci step, the system size is scaled by ϕ 3 implying that b = ϕ 3 . This gives us ν ≡ 1/y δ = 1.

E. Spin exponent for the quasiperiodic quantum Potts chain
At the critical point, the 2-point correlation function of the Potts order parameter scales as, Assuming for now that the critical point of the quantum Potts model is given by the Fibonacci fixed point discussed above, we can calculate this 2-point function analytically. Decimating the J i couplings in the RG leads to the formation of spin clusters which are locally ferromagnetic, while the decimations of a transverse field term h i freezes the corresponding spin (or effective spin cluster) in a given configuration with zero magnetization. We call the spin clusters that are not yet decimated active. The 2-point function is equal to the probability of the spins at 0 and r to be part the same active cluster, which is true if and only if all the spins between them are either part of the same cluster to which the spins at 0 and r will eventually belong to, or are decimated. This implies that the probability for the two spins to be in same cluster is equal to the probability for them not to be decimated at the length scale where there is no other active cluster in between them. If we denote by P (r) the probability for a given spin not be decimated under RG over a distance r, we thus have C(r) ∼ P (r) 2 . We then compute P (r) ∼ r −∆σ as follows: we start with system size r = F 3m with h i and J i taken from the same Fibonacci sequence, with h corresponding to odd links and J even ones. Let h 0 denote the magnetic field associated with the cluster under consideration. For a Fibonacci sequence, we know that all 'B' couplings are decimated in a single Fibonacci step. Thus we have following relation, P 0 represents the probability that the relevant field h 0 is a type 'A' coupling in the initial distribution, and is given by P 0 = ϕ/(1 + ϕ). (P 0 is the probability of getting a 'A' coupling at some particular site for a randomly chosen global phase.) P 1 is the probability that the h 0 coupling is a type 'A' link after 1 Fibonacci step. P 1 thus denotes the probability that the bond in question is part of a sub-pattern in previous Fibonacci step which gives birth to a 'A' bond in this step. For the standard Fibonacci sequence we know this sub-pattern to be 'ABABA'. Thus P 1 = probability of the link to be in the sub-pattern 'ABABA' in the previous Fibonacci step = 3ϕ/(3ϕ + 2). In a similar fashion we can show that all remaining P i = 3ϕ/(3ϕ + 2). This leads to, For a quasiperiodic Potts chain, it is much more natural to take h and J fields from independent distinct quasiperiodic potentials. The above analysis for the Fibonacci fixed point will not apply for such initial distribution. The system, however, still flows to sequences, but the fixed point oscillates between multiple sequences. Though the details of these sequences differ, we observe numerically that they have same universal features which leads to same spin and correlation length exponent as in above calculations (see main text). We remark that our argument only used the probabilities/ratio of various sub-patterns and letters, which could be same for a wide class of sequences. The details of the sequence is unimportant. As an example, for many choices of global phase we get a 3 letter sequence defined by the inflation rules A 1 → A 1 BA 2 BA 1 , A 2 → A 2 BA 1 BA 2 , B → A 1 BA 2 . This sequence has same ratios and probabilities of relevant sub-patterns as the Fibonacci sequence and hence the above calculations carry forward.

Flow to Fibonacci Sequence
The arguments above can be readily generalized to the symmetric MBL RG (see main text). The RG assumes alternating thermal and insulating blocks parametrized by "lengths", T /I . The RG flow is dictated by the decimation of the lowest length to get a new length given by the rule, where β I = 1/β T ≥ 1. We take the initial distribution of lengths to be a quasi-periodic potential given by T /I n = W T /I (1 + cos (2πϕn + θ)). The initial length distribution can be thought of as a sequence like in the case of XXX Supplementary Figure 2. Critical self-similar sequence for the MBL transition RG showing the lengths n after 3 Fibonacci step for different values of β I , and δW = 0. The red (blue) dots represent the thermal (insulating) blocks. We see that on increasing β I slightly, the length of thermal blocks become greater but the the sequence structure is not broken. However, when β I is increased further, low thermal blocks become greater than high insulating blocks (see β I = 1.68).
spin chain, and for the symmetric case of β I = 1 we get a Fibonacci sequence under RG with the following values, (1 + cos(2πϕ(n − i))) (7) since sin(F 3m+3 πϕ) = 0 as F 3m+3 ϕ ≈ F 3m+4 . Similarly one can show that, , implying that under RG the initial length distribution flows to a fixed point characterized by infinite difference in length. Thus the assumption that n n±1 gets better with the RG flow.
As we increase β I slightly, the thermal blocks are "favored" by the RG rules. This changes (7) and (8) slightly to give (ignoring fluctuations), where δ T,I 1,2 > 0, and represents the fact that the size of insulating blocks are reduced due to β T < 1 and the size of Supplementary Figure 3. Left: Probability of getting an insulating block at the end of RG is plotted against δW ≡ W I − W T for N initial blocks. The '+' markers correspond to β T = β I = 1 (symmetric case), while the 'o' markers correspond to β I = 10/7. The data for both β I overlap, as expected from the arguments given in the text. Right: Finite size collapse for β I = 1 with ν = 1. The error bars represent standard error.
thermal blocks are increased due to β I > 1. Under RG we will have, The first term in the above expression increases by a factor of ϕ 3 and the correction due to the asymmetry between thermal and insulating blocks also increases by an almost same factor. Thus if initially the corrections are small enough we see that the sequence structure is not broken under the RG flow. But for a large asymmetry, δ 1,2 can be large enough to destroy the sequence pattern (as we get B T > A I ). In fact studying the dependence of δ T,I 1,2 on β I , we find that for 1 < β I < 1.6 we have the same Fibonacci self-similar sequence as in symmetric case. This can be confirmed numerically by running the RG for different values of β I . In the left pannel of Fig. 3 we compare transition data for β I = 1 and β I = 10/7. They are overlapping perfectly in agreement with our argument.

Correlation length exponent
We now compute the critical exponent ν in the symmetric case (β I = 1). The mechanism responsible for driving the system to a phase (either thermal or MBL) is by introducing an asymmetry in the initial distribution via W I = W T . This leads to formation of "defects", see Fig. 4.
Let us assume that W T = 1 and δ ≡ (W I − 1) > 0 is small, so that the change in insulating blocks close to the red dash line in Fig. 4 can be written as, W I (1 + cos (2πϕn + θ)) ≡ 1 + cos (2π(ϕn + ∆) + θ), where ∆ is a function of δW with ∆ ∝ δW for small δW (provided the derivative of the initial potential is well defined at that point), see Fig.  4 The phase plays no important role so we take θ = 0 for simplicity. At criticality (δW = 0) and at the boundary of the minima region, i.e on the dashed red line in Fig. 4, two adjacent insulating and thermal blocks have the same length, that is, if n lies on the red line then {nϕ} = 1 − {ϕ}/2 and {(n + 1)ϕ} = {ϕ}/2. Now, imagine a local minimum insulating block slightly below the dashed line, i.e {nϕ} = 1 − {ϕ}/2 − γ, then the adjacent thermal block will be above the red line by same amount with {(n + 1)ϕ} = {ϕ}/2 − γ. As we move away from criticality δW = 0, we have If −γ + ∆ > γ then the insulating block will become larger than the neighboring thermal block. This implies that any insulating block for which Figure 4. Formation of a defect in the sequence. The green dots represent thermal blocks and red dots are insulating. The red dashed line represents the boundary for a block to be a local minimum: all the local minima are below that line for W T = W I = 1. Here 4 consecutive blocks are shown. As we increase W I , the thermal (red) dots go up. The initial potential (dashed black line) gets transformed into black line, showing the shift of the minima from the red dot to the green dot.
leads to a defect.
To compute the correlation length, suppose that the block at n 0 is arbitrarily close to the red line and satisfies the above criterion for creating a defect. We are now interested in the position of the next defect in the sequence pattern. Let n 0 + k be the position of the next defect. We want k to satisfy −∆/2 < kϕ − l < 0 for some integer l. Since n 0 + k is the first defect after n 0 , we should also impose the constraint that 0 < pϕ − q < 1/2 or −1/2 < pϕ − q < −∆/2, ∀p < k and for some q ∈ Z. We know by Diophantine approximation that the quantity |aϕ−b|, where a, b are integers, is lowest among all the integers c < a if and only if a is a Fibonacci number. Thus k has to be a Fibonacci number, k = F i0 with kϕ − F i0+1 = (−ϕ) −i0 and the constraint that ϕ −i0 < ∆/2 < ϕ −(i0−1) . (If kϕ − F i0+1 > 0, we can simply take the next Fibonacci number with k = F i0+1 .) The distance between defects sets the correlation length ξ = F i ≈ ϕ i ∼ ∆ −1 ∼ δW −1 , which implies that ν = 1. We checked numerically that ν = 1 leads to a perfect finite size collapse of the probability for the RG to end with an insulating block (Fig. 3).

B. Singular example:
The above argument for ν = 1 suggests that more singular potentials can lead different exponents ν > 1. When we wrote ∆ ∝ δW , we implicitly assumed that the derivative of the potential is not singular at the point where {ϕx + θ/(2π)} = {ϕ} 2 . However there exist some functions for which this is not true. For example, consider n = f (cos(2πϕn + θ)) with f (x) = 2 + sgn cos(2π {ϕ} 2 ) − x | cos(2π {ϕ} 2 ) − x|. f (x) is a bounded monotonic function in the range (−1, 1) and hence flows to sequence under RG, but it has a singular derivative at x = cos 2π {ϕ} 2 . This implies that δW ∼ √ ∆, so we expect ν = 2 in this case, as confirmed by numerical simulations of the RG (Fig. 5).
C. β I > 1 As mentioned above, slightly increasing β I is an irrelevant perturbation to the sequence, and does not alter the transition values. But if β I is increased beyond a certain threshold, the asymmetric perturbation can dominate the sequence values and break down the sequence pattern (see Fig. 2). We get different self-similar sequences as β I is increased. The next self-similar sequence is for 3 < β I < 12.75, see Fig. 6. This sequence is self-similar under one Fibonacci step. More concretely if the sequence size is F 3(m+1) then under RG, the sequence pattern will re-emerge when the system size becomes F 3m . To formalize this, we define the so-called substitution matrix of the sequence. If we define w 1 = A T , w 2 = A I , w 3 = B I , w 4 = B T , then the substitution matrix M is given by M ij = #w i in the inflation rule for w j .
For the Fibonacci sequence discussed above for the case of β I close to 1, the inflation rules were A T → A T B I A T B I A T , B T → A T B I A T , and so on. This corresponds to the substitution matrix, Let w = [#w 1 , #w 2 , #w 3 , #w 4 ] be the vector representing the number of various letters in the initial word sequence. Let λ i be the eigenvalues of M (with λ 1 the largest value) and v i being the corresponding eigenvectors. We assume that we can write w = i a i v i with some real coefficients a i . The vector M w denotes the number of letters after a single application of inflation rules. Thus after n inflations, the total number of letters is given by, implying that the number of letters changes by a factor of λ 1 under a single application of the inflation rule. For (9) λ 1 = ϕ 3 , implying that under one Fibonacci step, the number of blocks is decreased by a factor of ϕ 3 .
Moving back to the sequence in Fig. 6 for 3 < β I < 12.75, the substitution matrix is given by Supplementary Figure 6. Critical self-similar sequence after 3 Fibonacci step for different values of β I at the critical value of δW ≈ 1.97. Red dot represent the thermal blocks and blue represents insulating. We see that on increasing β I , the length of thermal blocks become greater and length of insulating blocks smaller (relative to one another). We checked that the sequence structure is preserved at all Fibonacci steps. Eventually on increasing β I further the sequence pattern will get destroyed, i.e the smaller thermal blocks will be larger than the largest insulating block.
with the largest eigenvalue given by ϕ 3 . This shows that the sequence repeats itself after a single Fibonacci step. We remark that the range of β I over which the above sequence is well defined can be identified semi-analytically. We can obtain the lengths of the block after the 1st Fibonacci step as a function of β I by summing over the cosines (like we did for symmetric case). Then the lengths at further steps can be calculated as follows, Eq (11) can be calculated numerically and we can check that the sequence pattern, A T m > A I m , B T m , A I m and A I m > B I m , is true at all Fibonacci steps for 3 < β I < 12.75.
As we increase β I further, we get different sequences with different periods of self-similarity. For 62 < β I < 78, we have another sequence whose substitution matrix is given by, The largest eigenvalue of M 3 is ϕ 6 . This implies that the sequence repeats itself after two Fibonacci steps, i.e if the system size is N = F 3m then under RG flow the sequence re-emerges when the system size is F 3(m−2) . This can be checked numerically (see Fig. 7). A nice finite size collapse for the transition is observed only for system sizes of the form F 3×(2m+1) , in contrast to the previous two transitions, for which collapses were obtained for system sizes of the form F 3m . We have also observed more complicated sequences for larger values of β I , with a period of 3 Fibonacci steps. Studying such sequences numerically is however very challenging as extremely large systems are needed. We expect that as one increases β I , the transition is controlled by larger and larger sequences with complicated inflation rules and scaling of the form ϕ 3k . Our defect argument for ν = 1 is however quite general, and applies generally to all these self-similar sequences.