Enhanced Superconductivity in Quasi-periodic Crystals

We study superconductivity in a family of one dimensional incommensurate system with $s$-wave pairing interaction. The incommensurate potential can alter the spatial characteristics of electrons in the normal state, leading to either extended, critical, or localized wave functions. We find that superconductivity is significantly enhanced when the electronic wave function exhibits a critical multifractal structure. This criticality also manifests itself in the power-law dependence of superconducting temperature on the pairing strength. As a consequence, an extended superconducting domain is expected to exist around the localization-delocalization transition, which can be induced by either tuning the amplitude of the incommensurate potential, or by varying the chemical potential across a mobility edge. Our results thus suggest a novel approach to enhance superconducting transition temperature through engineering of incommensurate potential.


I. INTRODUCTION
Electronically incommensurate potential appears in many condensed matter systems. Prominent examples include quasicrystals, borken symmetry with incommensurate order parameters, and the Moiré superlattice in twisted van der Waals heteroustructures. Because of the incommensurability between the emergent superstructure and the underlying lattice, the crystal momentum is no longer a good quantum number, which invalidates the conventional band-structure description of electronic states. More importantly, incommensurability can have significant effects on the electron eigenstates. For instance, incommensurate potential can render the electronic wave functions localized or critical [1]. This has been demonstrated in the Aubry-André model, a canonical system for studying the incommensurability-induced electron localization, and its variants [2][3][4] and quasicrystal systems [5][6][7].
Collective electron behaviors are also expected to be modified by the presence of incommensurability, due to the altered nature of single-particle wavefunction. Indeed, superconductivity [8] and unusual quantum critical state [9] have been observed in quasicrystals. The recent experimental observation of superconductivity and correlated insulating states in the twisted bilayer graphene (TBG) with incommensurate structure is another example [10][11][12][13]. Although the single particle physics in TBG can be satisfactorily described by a continuum model neglecting the incommesurability of the Moiré pattern [14,15], the role of incommesurability on many body states remains unexplored [16]. Motivated by these recent experimental progress, we study the superconductivity in a family of quasi-periodic systems both in one and two dimension.
In quasi-periodic systems, the electronic states can be categorized into extend, localized, and critical states depending on the spatial characteristics of the wave functions. In the extended state, the wave function spreads extensively over the whole system even in the thermodynamic limit, and are analogous to the Bloch states in crystals. The localized state, on the other hand, exhibits a wave function that is confined to only a finite number of lattice sites. Most interestingly, a multifractal, self-similar structure emerges in the wave function of the critical state [17]. The different nature of these elec-tron eigenstates also highlights a trade-off between the pairing strength and phase coherence of superconductivity. On one hand, although superconducting pairing can be maximized locally through confinement of electrons, superconductivity is disrupted due to the localized condensates. On the other hand, while a better phase coherence can be maintained by an extended wave function, delocalized electrons in such a state do not take full advantage of the short-range pairing interaction. As a consequence, the superconducting transition temperature T c is exponentially weak according to the BCS theory. This implies that T c may be enhanced in the case of critical states by optimizing the local pairing interaction while maintaining the long range phase coherence. This is indeed the case as will be revealed below.

II. MODEL
We study a one dimensional s-wave superconductor with an incommensurate potential, described by a Hamiltonian H = H 0 + H sc , with Here c † iσ (c iσ ) is creation (annihilation) operator of electron with spin σ on the i-th site of a periodic chain, µ is the chemical potential, U i is the on-site potential, and t is the nearestneighbor hopping constant, which is set to t = 1 for convenience in the following discussions. The Aubry-André (AA) model corresponds to an incommensurate U i = J cos(2πQx i ), where Q is an irrational number and x i is position of the i-th site, so that the local potential becomes incommensurate with the underlying lattice. In this study, Q is set to be the golden ratio, ( √ 5 − 1)/2, and is approximated by the Fibonacci sequence Q ≈ F n−1 /F n , where F n is the n-th Fibonacci number. We consider half filling by tuning µ. H sc describes the s-wave superconducting coupling and V is the pairing strength. The AA model, described by H 0 , exhibits a self-duality and a sharp localization-delocalization transition driven by J. It displays a spectrum consisting entirely of extended states for J < 2, and of localized states for J > 2. The quantum critical point J = 2 is characterized by a self-similar spectrum with all eigenstates becoming critical [18].
Standard Bogoliubov-de Gennes (BdG) method is used to solve this system [19]. The BdG Hamiltonian is where ∆ i = V c i↑ c i↓ is the local pairing amplitude. Here H eff can be diagonalized by Bogoliubov transformation, where γ † n and γ n are the creation and annihilation operators for Bogoliubov quasiparticle at state n and the prime sign means the sum is over all positive quasiparticle state E n > 0. The u and v coefficients are obtained from the BdG equations, where

III. BDG RESULTS
First we show the Bogoliubov-de Gennes (BdG) calculation results for the 1D s-wave superconductor under incommensurate potential U i = J cos(2πQx i ). The results of local order parameter, probability distribution of local order parameter and density of states, are shown in Fig. 1. It can be seen that for the extended states with J = 0, the system is a standard homogeneous s-wave superconductor. When the system is critical at J = 2, superconducting order parameter oscillates in space as evidenced by double peaks in the distribution P(∆). In the localized region with J = 4, there are superconducting islands with locally enhanced superconductivity separated by weak superconducting regions. In all cases, the spectrum is gapped around the chemical potential E = 0.
Since the emergence of superconductivity requires the phase rigidity of the Cooper-pair condensates in an inhomogeneous state, here we use the superfluid stiffness D s to characterize the long-range phase coherence. It is given by [19][20][21] where K x is the averaged kinetic energy and is the retarded correlation function of the particle current operator, The details on the calculation of D s based on the BdG method is presented in Appendix A. In the well localized phase, the global phase coherence is established by a weak Josephson type coupling between strong superconducting islands separated by weak superconducting regions [22]. The energy of the superconducting condensate can be approximated as E ∝ −D s i, j cos(θ i − θ j ), where θ i is the phase of the superconducting order parameter at i-th strong superconducting island [21]. D s increases with the amplitude of the superconducting order parameter in the strong superconducting islands and the overlap of the order parameter between these islands. Although strictly speaking, there is no long-range superconductivity order in 1D, our mean-field approach to the superconducting AA model should be viewed as a quasi-1D approximation to either 2D or 3D incommensurate superconductivity. (A true 2D model calculation will be presented in Sec. VI) With this understanding, superconductivity is destroyed by suppressing either the amplitude of the order parameter or the phase coherence. In the extended state as in the case of conventional BCS theory, T c is limited by the averaged amplitude of the order parameter over the whole system. In the localized state, T c is limited by phase fluctuation and is proportional to the zero temperature D s , which measures the coupling between different superconducting islands. Here we define two temperature scales: T c1 is the temperature when ∆ vanishes throughout the system, and T c2 = D s (T = 0) represents the energy scale of phase coherence. And we estimate the transition temperature T c of our system by min(T c1 , T c2 ). As J increases, the system becomes more spatially localized. T c1 increases because the system can take the advantage of local pairing interaction. Meanwhile T c2 diminishes as the superconductivity becomes more spatially localized. Therefore the superconductivity is limited by T c2 for a strong incommensurate potential. BdG calculation results of T c1 and T c2 of a system with size L = 233 are presented in Fig. 2. In Appendix D, We check that the finite size effect is negligible by comparing to the results with a larger L. T c1 increases monotonically with J. In the localized state, T c1 corresponds to the highest transition temperature among all superconducting islands. As the electronic wave functions become more localized by increasing J, the electrons can take full advantage of the local pairing interaction, and as a consequence, local superconductivity is enhanced. On the other hand, T c2 is first enhanced with increasing J until a critical J c , whose origin is unclear. As J is further increased, T c2 starts to decrease due to the loss of phase coherence between spatially localized superconducting islands. The dependence of T c1 and T c2 on J indicates the existence of a superconducting dome near the localization transition at J = 2, as schematically depicted in the inset of Fig. 2(a). It is also worth pointing out the different nature of superconducting transition on the two sides of the dome. In the extended regime corresponding to small J, the system undergoes a superconductivity to metal transition due to a vanishing amplitude of Cooper pairing upon increasing temperature. In the localized phase, there is a temperature-driven superconductorto-insulator transition caused by the phase fluctuations of the superconducting order parameter. At finite temperature, there is no sharp distinction between metals and insulators and we expect a smooth crossover tween metallic and insulating state for temperature above T c .
The dependence of T c1 and T c2 on V are plotted in Fig. 2(b) for the three different types of electron eigenstates in AA model. For extended wave function, T c1 has a standard BCS exponential relation with 1/V. In the localized state, the relation between T c1 and V is almost linear. Interestingly, at the critical point J = 2, T c1 increases with V according to a power law: T c1 ∝ V 1.6 . On the other hand, T c2 ∝ exp(−1/aV) in the localized state. Therefore in the weak coupling limit V t, T c is exponentially weak in 1/V both in the localized and extended regions, while T c is enhanced significantly near the localization transition as it depends on V by a power law.
In the AA model, the electronic spectrum form bands for the extended states. Both for the critical and localized states, instead of form bands, the spectrum is point-like [23]. Therefore, it is likely that the chemical potential locates in the gap of the single-particle spectrum, and therefore a threshold V is required to trigger superconductivity.

IV. WEAK COUPLING THEORY
The BdG method is restricted to a large V because the superconducting coherence length ξ increases exponentially fast for a weak V [ξ ∼ exp(1/N 0 V) with N 0 the density of state]. This would require large system size L ξ, which is practically impossible. To reach the weak coupling limit V/t 1 and also to understand the dependence of T c1 on V, we provide analytical description based on Anderson's idea of pairing the time-revered eigenstates of the non-interacting system [21,24]. The non-interacting time-reversal symmetric Hamiltonian H 0 is bilinear and can be exactly diagonalized: H 0 |ψ ασ = α |ψ ασ , where α labels the exact eigenstates of H 0 . We rewrite H in this basis, d † ασ |0 = |ψ ασ and only consider the pairing interaction between time-reversed states, |ψ α↑ and |ψᾱ ↓ : where Here ψᾱ is the time reversal partner of ψ α . The linearized gap The characteristic of the normal state electronic wave function is contained in the M matrix.
For the extended states J < 2, wave functions extend over the entire lattice and the amplitudes of wave functions scale as the inverse square root of the system size L, |ψ α | ∝ 1/ √ L. Thus, M αβ is independent of α, β and scales as 1/L, which leads to a gap equation with the standard BCS form and therefore T c1 ∝ e −1/aV . (see Appendix B for detailed calculations) For the localized states, the wave functions are confined in small regions characterized by a localization length ξ l and they scale as The wave function has negligible overlap with wave functions of other states. As a result, only the diagonal terms of M αβ are important, M αβ = δ αβ /ξ l , which results in a linear dependence of T c1 on V, T c1 ∝ aV. The results of T c1 vs V at J = 4 is shown in Fig. 2 (b). The dependence of T c1 on V deviates slightly from a linear behavior because of the nonzero overlap of wave functions at different energies when J is not large.
In critical state J = 2, the spectrum is self-similar, which is characterized by a multifractal exponent α M and its distribution f M (α M ) [25]. Therefore one would also expect M αβ to be self-similar, i.e. M αβ ≡ M( α , β ) = b −η M( α /b, β /b). For simplicity, we have assumed that M is characterized by a single exponent η. Because the spectrum is discrete, this scaling transformation is valid only for discrete value of b, as shown in Fig. 3. The scaling property of M immediately leads to a power-law relation between T c and V. If the superconducting coupling strength V is scaled by a factor α, V → αV, the energy level must also be scaled by a factor b, → b , to maintain the form of the gap equation unchanged. The term in tanh is dimensionless, thus the temperature T must also be scaled by b. From this scaling argument, we can obtain the power-law dependence T c1 ∝ V 1/(1+η) .
The spectrum of AA model is self-similar at J = 2. At J = 2, the spectrum does not have any continuous bands and for a finite system there is no one-to-one correspondence between the original M matrix and the zoomed one. Therefore, the value of η is estimated by calculating the self-similar scaling of largest elements in each small blocks of M αβ . The scaling exponent of each block is then averaged to obtain the final η. For L = 1597 and L = 6765, the matrix can be rescaled by a scaling factor b ≈ 13.8 and the estimated exponent is η ≈ −0.19 which gives T c1 ∝ V 1.23 . The scaling analysis agrees reasonably well with the numerical fitting result T c1 ∝ V 1.6 . The slight deviation in the exponent could be caused by the finite size effect because extremely large system size is required to capture the self similarity behavior of M αβ with a rescaling factor b ≈ 13.8. The single exponent approximation to the scaling relation for the M αβ can also cause deviation. Now it becomes clear from Eqs. (12) and (14) that the fractal nature of normal state wave function renders the effective pairing interaction being fractal. As a consequence, T c depends on the bare pairing interaction by a power law function, and superconductivity is enhanced.

V. OTHER 1D MODELS
The power law dependence of T c1 on V is due to the self-similarity of M αβ , which can be demonstrated in the Fibonacci model with an onsite incommensurate potential given by m < x < m + 1 − Q, where m is an arbitrary integer. The electronic spectrum is always critical regardless the strength of the potential [17,23,26,27], see also Appendix C. It can be seen from Fig. 4 that in this model, T c and V always has a power law relation provided that the Fermi level is not in a gap of the non-interacting spectrum. In the AA model, all eigenstates have the same spatial characteristics, and the localization-delocalization transition is controlled by the strength J of the incommensurability. It is also possible in certain class of incommensurate models that the localized and extended states coexists in the spectrum and are separated by a mobility edge, at which the wave functions become critical. One can thus change the wave function characteristics by tuning the chemical potential µ. Therefore there can exist a superconducting dome as a function of µ near the mobility edge. We demonstrate this scenario explicitly using the generalized Harper model with a modulated incommensurate potential U i = J cos(2πQx ν i ). Without superconducting coupling, this model exhibit a continuous spectrum with mobility edges at ±(2 − J) for 0 < ν < 1 and J < 2 [28], see also Appendix C. The states at the band edge are localized and the states in the middle of the band are extended. The BdG results of a system of L = 987 show a enhancement of superconductivity near the mobility edge forming a superconducting dome, contributed from the states near the mobility edge which are critical. Hence, there is a power law relation between T c and V at µ = −1.0. Superconductivity is suppressed when the chemical potential is tuned to localized region. Our BdG results for µ = −2.0 show that T c1 depends on V by a power law with a smaller exponent. This power law dependence is due to the relatively large V required by the BdG calculations, where both the localized state and critical states contribute to superconductivity. In the weak coupling limit when only localized states participate in the pairing, T c1 scales with V linearly according to Eq. (14). We next consider a model where the localized and extended states are separated by an energy gap in the spectrum. This is realized using the double cosine potential U i = J 1 cos(2πQx i ) + J 3 cos(6πQx i ) [17]. In this case, no critical state exists, and the dependence of T c1 and V follows either that for the extended states or for localized states as shown in Fig. 6.

VI. 2D AA MODEL
To further support our conclusion in 1D, we perform additional calculations in 2D. The 2D system also allows one to study the superconducting transition in the presence of thermal fluctuations. We consider s-wave superconductivity in a generalized Aubry-André model in 2D with an incommensu- rate potential Here, x i and y i are 2D coordinates of the i-th site. The localization transition occurs at J = 2 when superconductivity is absent [29]. We calculate T c using the weak coupling theory described in Eqs. (12) and (14). As shown in Fig. 7, at the localization transition, T c ∝ V 1.2 . The power law dependence of T c on V means an enhancement of T c in comparison to the standard weak coupling BCS theory results for a uniform system. We note that T c1 drops to zero quickly for a very small V. This is because of the finite size effect, when V is comparable to the discrete single particle spectrum gap. The results for a large V deviate from the power law behavior because the weak coupling approximation in Eqs. (12) and (14) breaks down.

VII. EFFECT OF COULOMB INTERACTION
In the critical or localized states, the effect of Coulomb interaction is also enhanced [30,31], similar to the pairing interaction. The effect of Coulomb interaction can be introduced by an energy dependent pairing strength where ω D is the Debye frequency and ω c is the frequency associated with the Coulomb interaction. V p < 0 is the attractive interaction and V c > 0 is the repulsive Coulomb interaction. We take ω D = 0.3, ω c = 0.5. The results of T c1 vs J of superconducting AA model with Coulomb interaction are shown in Fig. 8. T c1 for a given J is suppressed by the Coulomb interaction, and it is enhanced when the system is tuned to the more localized side at a given V c by increasing J. Note that T c1 is determined from the amplitude of the order parameter. When the system enters the localized region J > 2, the superconductivity is limited by superfluid stiffness, which is suppressed due to localization. Thus, a superconducting dome around the localization transition is expected even in the presence of Coulomb interaction.

VIII. DISCUSSION AND CONCLUSIONS
As shown in Appendix E, the normal state density of state (DOS) at Fermi energy is modified by the incommensurate potential, and hence affects T c . However, the change of DOS cannot explain the power law dependence of T c on the pairing interaction. So far, we have mainly focused on the effect of an incommensurate potential on T c , which is determined by normal state wave functions and spectrum. When superconductivity is fully developed far below T c , superconductivity can affect the localization transition by gapping the quasiparticle energy spectrum. It is possible that the localization transition at T = 0 is completely masked by superconductivity.
Let us discuss the relation of our work to others. Similar phenomenology has been discussed in superconductors with random disorders, where a power law dependence of T c on V is found based on scaling analysis near the localization transition [33][34][35]. The enhancement of T c by disorders due to the multifractal electronic state was studied theoretically [36] and observed in experiments [37]. The effect of random disorder on superconductivity was studied by solving the BdG equation numerically in 2D. No enhancement of superconductivity was found because of the absence of localization transition in the standard Anderson model in 2D. [21,38] A different mechanism for the enhancement of T c due to the enhancement of density of state by disorders was studied in Ref. [39]. In Ref. [40], it is argued that impurities can cause spatial modulation in the pairing potential, which enhances T c . In these cases, the single particle spectrum is continuous, which is different from that in the AA model studied here. We remark that the quasi-periodic potential has weaker effect on the localization of electronic wave function than the random disorders. This allows us to study the enhancement of the superconductivity near the localization transition or mobility edge in 1D models, which is not possible for random disorders, see Appendix F for more detailed discussions.
To summarize, we study the effect of incommensurate potential on 1D s-wave superconductors and found an enhancement of superconductivity near the localization transition in a class of quasi-periodic crystals. At the localization transition, T c depends V by a power law, which gives rise to a superconducting dome near the localization critical point. In the region with extended states, superconductivity is destroyed by the suppression of the amplitude of the superconducting order parameter; while in the localized states, superconductivity is killed by the fluctuations of the phase of the superconducting order parameter. Our results suggest a promising routine to enhance T c of superconductors by incommensurate potentials. As explained in the main text, the superfluid stiffness of the system is necessary for the estimation of critical temperature T c2 in the localized region. Here shows the derivation of superfluid stiffness expressed in the BdG framework. The derivation based on the method used in Ref. [19]. Consider a general Hamiltonian, We consider short range hopping. The particle current and the local kinetic energy associated with the x-oriented hopping can be written as, Here a = 1 is the lattice constant. The local conductivity can written in terms of these two operators, where · · · is the expectation value of the operator. Average over the spatial variable r r r i , where K x = 1 N i K x (r r r i ) . The correlation function is only a function of the time difference t − t , which allows a Fourier transform to frequency domain, where Π xx (q q q, t) = − i N θ(t) J P x (q q q, t), J P x (−q q q, 0) and Π xx (q q q, ω) = ∞ −∞ e iωt Π xx (q q q, t)dt. The superfluid stiffness D s is given by, Here e is set to 1 and dropped in the final expression. In the absence of spin-orbit coupling and other spin-flip scattering terms, the dimension of the BdG equation can be reduced from 4N to 2N. In this case, the BdG transformations are The prime sign above the summation indicates that only states with positive energy are included. Note that the reduction of Hamiltonian also divides the eigenvalues into two groups and the subscript ofñ 1 andñ 2 means that they correspond to different set of eigenvalues Eñ 1 and Eñ 2 . Thus there is a set of anti-commutation relations: γ † n 1 , γm 1 = δñ 1m1 , γ † n 2 , γm 2 = δñ 2m2 , γ † n 1 , γm 2 = γ † n 2 , γm 1 = 0. The kinetic terms can be written in terms of u and v as, Here f (E) is the Fermi function. The current-current correlation function can be written as where A n 1 n 2 σ (q q q) = i e −iq q q·r r r i j>i x j − x i t i j u n 1 * jσ u n 2 iσ − u n 1 * iσ u n 2 jσ , (A12) D n 1 n 2 σ (q q q) = i e −iq q q·r r r i j>i x j − x i t i j v n 1 jσ v n 2 * iσ − v n 1 iσ v n 2 * jσ .
For nearest hopping, j = i + 1 and the summation j>i x j − x i t i j is reduced to t. Using the above equations, one can calculate the superfluid stiffness D s by solving the BdG equation.
ψᾱ is the time reversal state of ψ α . Using the BCS mean-field approximation, ∆ β = V c β↓ c β↑ , the self-consistent equation can be therefore written using the M matrix where E β = ∆ 2 β + 2 β . This formulation allows us to study the relation between the critical temperature T c1 and the superconducting coupling strength V in the weak coupling regime which is inaccessible for the numerical BdG calculations.
For extended state, all wavefunctions extend over the entire lattice and the amplitude of wavefunctions scales as the inverse square root of the lattice size, |ψ α | ∝ 1 √ L . Thus, M αβ becomes independent of α and β and scale as M αβ ∝ 1 /L . The order parameters satisfy ∆ α = ∆ β = ∆, and we obtain We can transform the discrete summation over states to an integration of energy by introducing the density of state (DOS), In the weak coupling limit, we can approximate N( ) by the density of state at Fermi surface N 0 and rewrite equation (B7) as where E = √ ∆ 2 + 2 . When the temperature approaches the critical temperature from below T → T − c , the order parameter goes to zero from above ∆ → 0 + and E ≈ . Introducing an integration cutoff ω c T c , we have where γ is a constant number. Therefore, the critical temperature is given by In the localized state, the wavefunctions are confined in small regions characterized by a localization length ξ l and they scale as |ψ α | ∝ 1 √ ξ l . The wavefunction has almost no overlap with wavefunctions of other states, as a result, only diagonal terms of M αβ are important. Thus the M matrix has the form, M αβ = δ αβ /ξ l . The gap equation (B6) becomes, where E = √ ∆ 2 + 2 . When ∆, this equation only has zero solutions ∆ = 0. When ≈ 0, this equation becomes Near the critical temperature, the order parameter ∆ is small and we can obtain Therefore, the critical temperature is linearly proportional to V In the localized state, T c increases with J because ξ l decreases with J, which is consistent with the BdG results. In general cases, the M matrix can be calculated numerically using the wavefunctions obtained from diagonalizing the Hamiltonian. T c corresponds to the largest eigenvalue of the linearized gap equation. In practice, we find the corresponding V for a given T c . This method gives accurate T c same as the result obtained from BdG calculation. Results for T c vs V obtained by M matrix for the AA model in the extended, critical and localized phases are shown in Fig. 9. I n = x i |ψ n (x i )| 4 x i |ψ n (x i )| 2 −2 , where ψ n (x i ) is the n-th eigenfunction of H 0 . I n is finite for a localized state but vanishes as 1/L d for an extended state. Here L is the linear system size and d the spatial dimension. The states near the mobility edges are critical which leads to a power-law dependence between T c and V when µ ≈ ±1 as shown in the main text.

Fibonacci model
The Fibonacci model has an incommensurate potential: U i = U(Qx i ) and U(x) = −J for m − Q x m, J for m < x < m + 1 − Q, where m is an arbitrary integer and Q = √ 5 − 1 /2. The key characteristics of this model is that it is always critical regardless the strength of J. The spectrum exhibits self-similarity, see Fig. 11. Here we check the finite size effect. We calculate T c1 (temperature when the amplitude of superconducting order parameter vanishes) and T c2 (temperature when the superfluid stiffness vanishes) for different system sizes, and the results are shown in Fig. 12. It can be seen that T c1 and T c2 converge to a fixed value very quickly as one increases L. Therefore, the system size with L = 233 we used in the main text has negligible finite size effect for the V we used. When V is reduced, larger system size is required because the superconducting coherence length increases when V is reduced.  Appendix E: Density of state in the normal state Superconducting transition temperature T c depends on the normal state density of state (DOS) at the Fermi energy. Here we calculate the normal state DOS in the presence of an incommensurate potential using the Aubry-André model in Eq. (1). As displayed in Fig. 13, the DOS at Fermi energy is increased in the presence of an incommensurate potential. However, this increase of DOS cannot explain the power law dependence of T c on the pairing interaction.