Supercurrent Distribution in Real-Space and Anomalous Paramagnetic Response in a Superconducting Quasicrystal

We theoretically study the real-space distribution of the supercurrent that flows under a uniform vector potential in a two-dimensional quasiperiodic structure. This is done by considering the attractive Hubbard model on the quasiperiodic Ammann-Beenker structure and studying the superconducting phase within the Bogoliubov-de Gennes mean-field theory. Decomposing the local supercurrent into the paramagnetic and diamagnetic components, we numerically investigate their dependencies on average electron density, temperature, and the angle of the applied vector potential. We find that the diamagnetic current locally violates the current conservation law, necessitating compensation from the paramagnetic current, even at zero temperature. The paramagnetic current shows exotic behaviors in the quasiperiodic structure, such as local currents which are oriented transversally or reversely to that of the applied vector potential.


I. INTRODUCTION
A quasicrystal is a solid that lacks translational symmetry but exhibits a diffraction pattern with sharp Bragg peaks and a rotational symmetry forbidden in periodic lattices [1,2].The quasicrystalline structure results in exotic electronic states, such as critical states [3][4][5][6][7][8][9], which are distinct from those of conventional periodic crystals.Recently, an experimental work discovered bulk superconductivity in a Bergmann-type Al-Zn-Mg quasicrystalline alloy (T c ∼ 50 mK) [10].More recently, superconductivity has also been reported in a van der Waals layered quasicrystal Ta-Te (T c ∼ 1 K) [11].The superconductivity in quasicrystals poses new questions since they do not possess fundamental prerequisites such as the Fermi surface in the conventional Bardeen-Cooper-Schrieffer (BCS) theory [12] due to the absence of translational symmetry.In previous theoretical works, such superconductivity has been studied by considering the attractive Hubbard model on quasiperiodic lattices.These studies showed that superconducting pairing is inhomogeneous, with real-space distributions of the sitedependent local electron density and superconducting order parameter [13][14][15][16][17][18][19][20][21][22][23][24].More interestingly, it has been pointed out by two of the present authors that non-BCS type superconductivity, comprised of Cooper pairs with finite center-of-mass momentum, exists in the weakcoupling region [17].
The interplay between the quasiperiodicity and superconductivity was studied in earlier works on quasiperiodic pinning arrays in periodic superconductors [25][26][27][28][29], as well as quasiperiodic networks of ordinary supercon-ducting wires [30][31][32][33][34][35].These were studies for superconductivity in artificially fabricated quasiperiodic structures.In contrast, we investigate the electromagnetic response of superconducting quasicrystals at the atomic level.In particular, supercurrent that flows in response to a uniform vector potential such as the Meissner current is a basic property that has however been scarcely explored.In a periodic system with a simple unit cell, it is obvious that the local supercurrents are uniformly distributed in the lattice, due to the homogeneity of the superconducting state in this case.In a simple crystal, each of the paramagnetic and diamagnetic components [36] of the local supercurrent is also uniform.In contrast, we will show that in the quasicrystal, where both the local electron density and superconducting order parameter are spatially varying, the local supercurrent exhibits a nontrivial spatial dependence as well.
In this study, we consider a two-dimensional quasiperiodic structure and the local supercurrent flow induced by an external uniform vector potential.We investigate the attractive Hubbard model on the Ammann-Beenker structure [37][38][39] by means of the Bogoliubov-de Gennes (BdG) mean-field theory.First, we formulate the expression of the local supercurrent under the uniform vector potential and discuss its real-space distribution on the structure.To clarify how the inhomogeneity of the superconducting state affects the supercurrent flow, we further decompose it into the paramagnetic and diamagnetic components.We then discuss the dependence of the local supercurrent on (i) average filling, (ii) temperature, and (iii) angle of the applied vector potential.Interestingly, we find that the non-uniform diamagnetic current can locally violate the current conservation law, i.e. have a non-zero divergence.The paramagnetic current flows so as to re-establish the conservation of the local current.This leads to non-uniform supercurrent distributions which are unique to quasiperiodic systems.We arXiv:2310.03484v1[cond-mat.supr-con]5 Oct 2023 find furthermore that the paramagnetic currents do not vanish at zero temperature, an anomalous property that was observed earlier by Liu et al. [23] in the site averaged value.
The rest of this paper is organized as follows.In Sec.II, we introduce the model Hamiltonian and explain our theoretical approach.We discuss the distribution of the local supercurrent and its dependencies on the average filling, temperature, and the angle of the applied vector potential in Sec.III.A brief summary is given in Sec.IV.The relation between formulations in our previous study [40] and the present one is explained in Appendix.

II. MODEL AND METHOD
This study is carried out on the Ammann-Beenker structure (Fig. 1), which is a two-dimensional quasiperiodic tiling with an eight-fold rotational symmetry [37][38][39].For our numerical calculations, we use a square approximant of the perfect infinite tiling, consisting of N = 1393 sites.This square approximant of the Ammann-Beenker structure was generated by the cut-and-project method [41].Here, we adopt a vertex model, where an atomic orbital is placed on each vertex of the Ammann-Beenker tiling.The coordination number Z i at each site ranges from 3 to 8 and the vertices can be categorized into six classes if one does not distinguish two geometries of Z i = 5 [42,43].
We consider an attractive (U < 0) Hubbard model [44] to study s-wave superconductivity in this Ammann-Beenker structure, as was done in previous studies [17,18,21].To study the local supercurrent in the presence of a vector potential, we include A(r) as the Peierls phase in the transfer term of the model Hamiltonian [45].Thus, the Hamiltonian is given by Here, ĉ † iσ (ĉ iσ ) creates (annihilates) an electron of spin σ at site i.We suppose a finite electron-transfer integral t only between the nearest neighbor sites (denoted by ⟨i, j⟩) connected by an edge of a square or a rhombus and set it as the unit of energy.In the noninteracting limit, the energy width of the density of states is about 8.5t [46].We define a local electron density n i = σ ⟨n iσ ⟩ with niσ = ĉ † iσ ĉiσ .The chemical potential µ is tuned to fix the average electron density n = i n i /N where N is the system size.We fix the attractive interaction strength to U = −3, and select the averaged electron density n = 0.3, 0.5, 0.7, and 0.9 to avoid the delta-function singularity in the density of states at the half-filling due to confined states [47].In the case of a uniform vector potential A, the Peierls phase in Eq. ( 1) can be rewritten as where r ij = r i − r j = ae n is the bond vector between sites i and j.Here a is the bond length and the unit vectors e n = (cos ϕ n , sin ϕ n ) with ϕ n = 2nπ/8 (n = 0, ±1, ±2, ±3, 4) correspond to the eight permitted bond orientations on the tiling, where we take ϕ n = 0 as the x direction.We henceforth assume that the uniform vector potential A is applied parallel to the plane of the Ammann-Beenker structure and evaluate the local supercurrent induced by A on each bond r ij .We control the direction of the vector potential A = −|A|(cosθ, sinθ) by changing the angle parameter θ in the range 0 ≤ θ < π 4 .Using the mean-field approximation, Hamiltonian in Eq. ( 1) is reduced to where K σi,j = −tδ ⟨i,j⟩ e −iA•rij +(U n iσ −µ)δ i,j is a kinetic term with n iσ = ⟨n iσ ⟩ and ∆ i = −U ⟨ĉ i↑ ĉi↓ ⟩ is the sitedependent superconducting order parameter.δ ⟨i,j⟩ is the Kronecker delta that counts only between the nearest neighbor sites [21].We note that the Hartree term U n iσ needs to be explicitly incorporated above because it has a site dependence and hence cannot be absorbed into the chemical potential term [17,18,21,48,49].The Hamiltonian of 2N ×2N matrix in Eq. ( 3) is diagonalized through the Bogoliubov transformation and we obtain the Bogoliubov-de Gennes (BdG) equation [21,40,[50][51][52] j Ĥi,j Here, E ϵ denotes an eigenenergy of the BdG Hamiltonian Ĥi,j and u ϵ (r i ), v ϵ (r i ) denote wave functions on the site i.The index ϵ distinguishes eigenstates of the BdG Hamiltonian, which runs over 1 to 2N , including eigenstates with both positive and negative E ϵ .As the selfconsistent condition, the gap equation and local electron density for each spin are obtained as [53,54] with f (E) = 1/(e E/T + 1) is the Fermi-Dirac distribution function at temperature T .Using only positive E ϵ , we obtain the conventional formula [50], where we have used the particle-hole symmetry, (u ϵ (r i ), (6).
The local supercurrent from a site j to i is given by It can be divided into the paramagnetic current and the diamagnetic current so that J j→i = J para j→i + J dia j→i .Here, we have defined the paramagnetic (diamagnetic) component as an even (odd) function of A in the expression of the local supercurrent.In the weak limit of the vector potential, trigonometric functions in Eqs. ( 11) and ( 12) can be reduced to 1 and A • r ij = |A||r ij | cos α, respectively reproducing the conventional definition of J dia j→i and J para j→i used in previous studies [23].The angle parameter ϕ n of r ij specifies the flow direction of the local current.Here, we define an angle α = θ − ϕ n between the applied vector potential and the bond vector and call cos α a bond factor.We note that Re{⟨ĉ † iσ ĉjσ ⟩} = (⟨ĉ † iσ ĉjσ ⟩ + ⟨ĉ † jσ ĉiσ ⟩)/2 represents the effective bond strength between the site i and j (̸ = i), and Im{⟨ĉ † iσ ĉjσ ⟩} = (⟨ĉ † iσ ĉjσ ⟩ − ⟨ĉ † jσ ĉiσ ⟩)/2i gives the net transfer from the site j to i.In Eqs. ( 10)−( 12), ⟨ĉ † iσ ĉjσ ⟩ is obtained from the eigenstate of the BdG equation (5) as In this study, we have chosen the amplitude of the uniform vector potential to have the value |A| = 0.005.Since the local supercurrents are linear as a function of the vector potential in the weak |A| limit, changing the value of the external vector potential will not result in qualitative changes in our results.

A. Real-space distribution of the local electron density and the superconducting order parameter
We begin by considering the case when the uniform vector potential A is parallel to the x-axis, i.e., A = (−|A|, 0) by setting θ = 0. Before discussing the supercurrent distribution, we study the inhomogeneous distribution of local quantities.Figure 2 presents real-space distribution of the local electron density n i and the superconducting order parameter amplitude |∆ i | after the self-consistent calculation of Eqs. ( 5)-( 8) at T = 0.01 for the fillings n = 0.3 and 0.7.The figure zooms in a small region consisting of about 100 sites, for which one clearly sees inhomogeneous spatial distributions of n i and |∆ i |, which moreover exhibit an approximate eight-fold symmetry as reported in Ref. [18].Classifying the vertices by the coordination number Z i , we plot the distributions of n i and |∆ i | against Z i in the right panels.We note that the values of n i and |∆ i | have variations even among the sites with the same coordination number since such sites have different next-nearest neighbor (or further neighbor) configurations.
For the filling n = 0.3, the right-hand panel of Fig. 2(a) shows that n i tends to increase with Z i .This behavior can be deduced from a property of the non-interacting model: when the Fermi energy lies below the main pseudogap, it is the sites of large Z i which are preferentially occupied [55].This leads, in the BdG equation in Eq. ( 7), to the factors |u ϵ (r i )| 2 being larger for larger Z i values.
In contrast, the local order parameter amplitude, shown in Fig. 2(b) does not increase monotonically with Z i but has a maximum at Z i = 5, as can be seen from the righthand panel.This maximum in |∆ i | can also be simply explained -according to Eq. ( 6), the order parameter amplitude is given by the product of |u ϵ (r i )| and |v ϵ (r i )|, which are increasing and decreasing functions of Z i respectively.This leads to the maximum at Z i = 5, which is the value separating low and high coordination sites in this tiling.As the filling is increased to larger values, for n = 0.7, the n i increases at all the sites, however, the differential increase is largest at the sites with smaller Z i , as shown in Fig. 2(c).Overall, the distribution range of n i for different Z i becomes narrower as filling n is increased, until at half-filling one reaches the uniform state n i = 1 for all sites.For n = 0.7, the local superconducting order parameter is enhanced compared to the case of filling n = 0.3.As shown in the right-hand panel of Fig. 2(d), the values of ∆ i are in this case largest on the sites with Z i = 3 and 4.This reflects the fact that, for a large filling, the non-interacting local density of states at these sites is significantly larger than that of the sites of large Z i [55], for reasons discussed in [56].These sitedependences of n i and |∆ i |, which hold even in the absence of the uniform vector potential A, affect behaviors of local supercurrent flow as we now discuss below.

B. Real-space distribution of the local supercurrent
Figure 3 shows the spatial distribution of the supercurrent in the case of θ = 0. Bonds are either parallel to this direction (ϕ n = 0), perpendicular (ϕ n = π 2 ), or at an angle of π 4 , which leads to large differences in the bond factor.In uniform systems such as a square lattice, the differences in J j→i = |J j→i | among the bonds are attributed only to the bond factor cos α for fixed |A|.It is just because J j→i ∼ J dia j→i at sufficiently low temperature.Therefore, the local supercurrent flows on the bonds where the bond factor cos α is non-zero, i.e., the bonds with |α| ̸ = π 2 .Also, the same current flows for the same α bonds, forming a one-dimensional flow distribution consisting of the respective local currents.These are well-known responses of uniform superconductors [36].
On the other hand, the local supercurrent J j→i in the inhomogeneous superconductor flows non-uniformly as shown in the left panel of Fig. 3(a) for n = 0.3, which is not determined simply by the bond factor.The overall tendency to flow along one-dimensional "channels" is similar to that of the uniform systems.These onedimensional channels having a cross-sectional width of a few lattice spacings are stacked along the y-direction.Notably, J j→i depends on the sites i and j even if the bonds have the same bond factor.In addition, there are small supercurrent flows even in the directions of |ϕ n | = π 2 .These features are characteristic of the quasiperiodic superconductor.
To understand the non-uniform distribution, we decompose J j→i into the diamagnetic current J dia j→i and paramagnetic current J para j→i as shown in the middle and right panels of Fig. 3(a).Since J dia j→i can be considered as a direct response to the vector potential A and has the bond factor cos α, J dia j→i of ϕ n = 0 is larger than that of |ϕ n | = π 4 .We note that J dia j→i = 0 for |α| = |ϕ n | = π 2 , which means that for these bonds J j→i = J para j→i .Such J para j→i flowing on the bonds perpendicular to A is unique to the non-uniform superconductor, and the presence of J para j→i prevents the formation of the one-dimensional channels.
As we show in Appendix B, results for the Penrose structure show that perpendicular currents are likewise present in that case.Indeed, perpendicular local currents can arise in quasiperiodic structures as these systems do not have translation invariance.There are no such currents on the square or honeycomb lattices (see Appendix B).
Furthermore, we have checked that such currents flow even when the local order parameter ∆ i and electron density n i in Eq. ( 5) are assumed to be uniform on all sites (self-consistency is not imposed).This shows that the non-uniformity of ∆ i and n i are not an essential condition for perpendicular currents to flow in this case.
Here, the existence of the paramagnetic current in these directions can be understood in terms of the conservation law of the local supercurrent, which is defined as that the total currents coming in and going out of each site should agree.To see this more clearly, in Fig. 4, we show the divergence (J dia i ) out − (J dia i ) in of J dia j→i at each site with classifying the sites by the coordination number Z i .The diamagnetic currents entering and leaving a site i are expressed as (J dia i ) in and (J dia i ) out , respectively.We note that the local supercurrent J j→i is conserved at any site in both the periodic and quasiperiodic systems.While the local current must be conserved, as required by gauge invariance, this constraint does not apply to the diamagnetic and paramagnetic parts taken separately.In periodic systems with uniform superconducting states, one finds nevertheless that the diamagnetic and paramagnetic currents are separately conserved.That is because they are proportional to scalar products of A and r ij and the summation of the most neighbor r ij is zero at all sites i.In contrast, in the quasiperiodic system, we observe that the diamagnetic currents are not locally conserved.This can be seen from the plot in Fig. 4 which shows that the divergence of the local diamagnetic current is not zero.This leads to the fact that the paramagnetic current is not locally conserved, either, in order to satisfy the local current conservation of J j→i = J para j→i + J dia j→i .In addition to the ϕ n dependence, we see in Fig. 3(a) a trend that J dia j→i becomes larger on bonds connected to the sites with larger coordination numbers such as Z i = 8, 7, and 6 where the local electron density n i is larger.As expected from the physical role of J para j→i , it flows in the opposite direction to that of J dia j→i .Therefore, the paramagnetic current also flows more on the bonds connected to the sites with larger Z i .This point will be further discussed in Sec.IIIC.Importantly, this paramagnetic current remains finite even at zero temperature as pointed out in Ref. [23] (see Sec. III D), contrary to the case of the uniform system.Moreover, J para j→i remains finite regardless of the flow directions.
Summing up, we have described the spatial distribution of the local supercurrent J j→i on the Ammann-Beenker structure for θ = 0 (vector potential along the x-axis).We find that J j→i is inhomogeneous, and takes different values on the tiling, even among the bonds sharing the same bond factor.In contrast to the case of periodic systems, J dia j→i itself is not locally conserved on this structure but is compensated by J para j→i .One of the consequences of this type of compensation is that J para j→i flows on the transverse bonds of |ϕ n | = π 2 where no diamagnetic current flows.We stress that this effect is observable only upon examining current patterns at a given node, that is at a local scale.

C. Filling n dependence
Since the distribution of n i and ∆ i changes significantly with the filling as shown in Fig. 2, the supercurrent distribution is also expected to change accordingly.First, we compare the spatial structure of the supercurrent for n = 0.7 in Fig. 3(b) and n = 0.3 in Fig. 3(a).In the case of n = 0.7, the distribution of J dia j→i becomes relatively uniform for the same |ϕ n |, and each J dia j→i is larger than that for n = 0.3.At the same time, J para j→i around the Z i = 8, 7, and 6 sites is strongly reduced from that for n = 0.3.Moreover, the much larger J para j→i flows in the direction of |ϕ n | = π 2 .We observed that this trend is particularly pronounced on bonds connected to the Z i = 4 sites, where n i increases significantly with n in Fig. 2

(c).
To see the n dependence more systematically, we plot J j→i , J dia j→i , and J para j→i against n in Fig. 5. Since ϕ n and −ϕ n are equivalent for θ = 0, the distributions are grouped by |ϕ n | in Fig. 5.The flow directions ϕ n of J para j→i is rotated by π from those of J dia j→i in Fig. 5(c) since it flows in the opposite direction of J dia j→i .In Fig. 5(b), we see a trend that J dia j→i of |ϕ n | = 0 and π 4 increases with n.Note that J dia j→i of |ϕ n | = π 2 is 0 due to the bond factor.In Fig. 5(c), while we do not find a clear trend in J para j→i of |ϕ n | = 0 and π 4 , we find that J para j→i of |ϕ n | = π 2 and its distribution range increase monotonically with n.This suggests that the conservation law of the diamagnetic current is further violated as the filling n increases.As shown in Fig. 4 for n = 0.7, the deviation from 0 becomes big compared to the case of n = 0.3.Interestingly, the divergence at  c).We note that J para j→i flows in the opposite direction of J dia j→i .The data points for each |ϕn| are slightly shifted in the horizontal direction for the sake of visibility.
distribution J dia j→i after dividing by the bond factor (for bond angles |ϕ n | = 0 and π 4 where cos ϕ n ̸ = 0).The relation J dia j→i ∼ Re σ ⟨ĉ † iσ ĉjσ ⟩ (from Eq. ( 12)) suggests that there may exist two types of simplified dependence.The rescaled supercurrent variable is thus plotted in two different ways in Fig. 6: as a function of √ n j n i (left-hand column), and as a function of |∆ j ∆ * i | (right-hand column), for four different values of the filling.As one can see in the figure, the blue (ϕ n = 0) and red dots (|ϕ n | = π 4 ) overlap, showing that the new variables J dia j→i / cos ϕ n are independent of the bond orientation.The plots show that systematic correlations do exist between the rescaled local currents and the local charge/order parameter in some limits.For small filling, n = 0.3 (top row), J dia j→i is positively correlated with √ n j n i , but is uncorrelated with the local superconducting order parameter amplitudes.For the large filling n = 0.9 (bottom row), the vice versa is true: the current is correlated with |∆ j ∆ * i |, but is uncorrelated with the local charges.Based on these numerical observations we conclude that at low filling the diamagnetic current on a given bond is approximately while at higher filling, the diamagnetic current on a bond is approximately given by Intermediate behaviors can be seen for n = 0.5 and 0.7, showing that both amplitude and phase variations are important in the generic case.The limiting behaviors for small and large fillings help to explain our observations: at small filling n, the first relation attributes the large J dia j→i around the sites with Z i = 8 and 7 in Fig. 3(a) to the large n i at such sites [Fig.2(a)].In the opposite limit of a high n, the second relation accounts for the observation in Figs.3(b) and 2(d) that J dia j→i flows well on bonds connected to the sites with Z i = 3 and 4, where The superconducting transition temperature T c for the Ammann-Beenker structure is found using the condition |∆ i (T c )| = 0 for all the sites.For U = −3 and filling n = 0.5, the value of T c = 0.344 for the Ammann-Beenker structure, which can be compared with a value of 0.333 for the square lattice on the same interaction strength and filling within our framework.
Figures 7(a), (b) and (c) show the temperature dependence of the local current and its dia-and para-components respectively.Black curves show the results for a square lattice (N = 900) with the same parameters.One sees that the local supercurrent in the tiling tends to zero as T approaches T c , in accordance with expectation.Fig. 7(b) shows that J dia j→i is almost constant as a function of temperature for all the bond orientations.Note that the diamagnetic current is identically 0 due to the bond factor for |ϕ n | = π 2 .As J para j→i cancels J dia j→i above T c , the supercurrent vanishes at T ≥ T c as shown in Fig. 7(c).With lowering temperature (T < T c ), J para j→i of |ϕ n | = 0 and π 4 decreases.This tendency is similar to the results on the square lattice.However, while the paramagnetic contribution vanishes at T → 0 in the square lattice, it exhibits a non-zero value even at zero temperature in the quasiperiodic structure, as pointed out in Ref. [23], for the site-averaged values.Our results reveal that it holds for all flow directions ϕ n .Remarkably, in the directions of |ϕ n | = π 2 , J para j→i increases on lowering T (< T c ).The existence of J para j→i even at T → 0 can be quali- Here m * , v, and e * respectively denote the mass, velocity, and charge of the Cooper pairs, and c denotes the light velocity.In the quasiperiodic systems, the Cooper pairs hold finite canonical momentum, as pointed out in Ref. [17].Therefore, the first term in Eq. ( 15) does not vanish and gives a finite contribution to J para j→i even at T = 0. We note that J para j→i flows in the opposite direction of J dia j→i .The data points for each |ϕn| are slightly shifted in the horizontal direction for the sake of visibility.
as θ increases.On the other hand, the θ dependence of J para j→i for n = 0.7 is clearly weaker as shown in the right panel of Fig. 9(b).We note that J para j→i in the direction of ϕ n = − π 2 , which is related to a back-flow as discussed later, decreases with θ for 0 ≤ θ ≤ π 18 , while it increases for π 18 < θ ≤ π 8 to cancel the increase of J dia j→i .The left panels of Fig. 9 show the distributions of the resulting J j→i for various angles θ.After changes of J j→i distribution with increasing θ, J j→i for ϕ n = 0 and π 4 reach the same distribution at θ = π 8 .In the case of n = 0.3, while the distribution of J j→i in the direction of ϕ n = 0 has weak θ dependence, max{J j→i } decreases with θ in the range 0 ≤ θ ≤ π 8 .On the other hand, for n = 0.7, the distribution range of J j→i for ϕ n = 0 expands as θ increases, with the increase of max{J j→i } for 0 ≤ θ ≤ π 8 .We have already noted that for θ = 0 the current can flow in the direction transverse to the vector potential.For example, J j→i for ϕ n = − π 2 is a transverse flow to A, in the results shown in Fig. 3.When θ ̸ = 0, a "backflow", i.e., the current satisfying J j→i • A ∝ cos α < 0 , occurs.Figure 10 shows an example of the back-flow in the case of n = 0.7 and θ = π 36 .In the left panels of Figs.9(a) and (b), except for θ = 0, J j→i for ϕ n = − π 2 is the back-flow.We find that this back-flow decreases with θ up to θ = π 8 .Since J para j→i flowing in the direction of ϕ n = − π 2 is larger than J dia j→i in the opposite direction, the back-flow appears in J j→i .
In this way, the back-flow comes from the paramagnetic component, which flows to satisfy the local current conservation of J j→i , compensating for the broken local current conservation of J dia j→i .Thus, the back-flow of J j→i is also one of the characteristics of the quasiperiodic superconductors.

IV. SUMMARY
We have studied the local supercurrent flow under the uniform vector potential on the Ammann-Beenker structure.To address this problem, we introduced the attractive Hubbard model, where the effect of vector potential is incorporated as the Peierls phase in the transfer term, and numerically analyzed it based on the self-consistent BdG mean-field theory.We decomposed the local supercurrent into the diamagnetic and paramagnetic current in our formulation in order to better understand the non-uniform spatial distribution.Our formulation for the local supercurrent is applicable not only to other quasiperiodic structures but also to general non-uniform structures with the periodic boundary condition.
We confirmed that the local electron density and superconducting order parameter are distributed nonuniformly with approximate 8-fold symmetry, as known in Refs.[17,18,21].The distributions greatly vary depending on the filling n.We clarified a spatial distribution of the supercurrent and its variation depending on (i) the average electron filling n, (ii) temperature T , and (iii) the angle θ of the applied vector potential.
Firstly, the diamagnetic current has a temperature dependence similar to that in the uniform systems, but the paramagnetic current has a finite value even at T → 0 for all flow directions ϕ n .We believe that such a phenomenon can be also realized in the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) states [57,58].However, in the quasiperiodic systems, proper adjustment of the magnetic field is unnecessary, and it would be easier to confirm this phenomenon through experiments.Secondly, as the filling increases, the vertical paramagnetic current increases, which is accompanied by a change in the distribution of the local electron density and superconduct- ing order parameter.Lastly, the local supercurrent flows even in the direction transverse to the applied vector potential.Furthermore, as the angle of the vector potential increases, back-flows are observed where the bond factor cos α is negative.In any case, the paramagnetic current behaves differently from that in periodic systems.This is because the diamagnetic current is affected by the distributions of the local electron density and the superconducting order parameter, and is not conserved locally.As a result, an excess amount of paramagnetic current has to be induced even at zero temperature as a counterpart to recover the local current conservation and contributes to characteristic local supercurrent behaviors on the quasiperiodic structure.
In conclusion, we have presented the first theoretical investigation of real space distributions of the supercurrent in a structure that does not possess translation invariance but is perfectly ordered.The novel spatial distributions of the supercurrent revealed in this study are the first step in understanding the response to magnetic fields and the Meissner effect in quasicrystalline superconductors.More detailed investigations of the distribution of screening currents under external fields are planned for future work.

V. APPENDIX A. Gauge transformation of Hamiltonian
The formulation of the local supercurrent in Sec.II differs from our previous formulation [40].Here we will see that these two formulations are equivalent through gauge transformations.
In Eq. ( 1), the following gauge transformation is applied to the creation and annihilation operators.from Eq. ( 23), the following expressions are obtained.
Thus, we see that the present formulation and that of the previous one coincide through the gauge transformations.Our intention behind the formulation of this study is to decompose the supercurrent into paramagnetic and diamagnetic components.

B. Supercurrents on the Penrose structure and honeycomb lattice
In this section, we show results for the Penrose structure, and for a simple periodic structure -the honeycomb lattice.These examples help to clarify the reasons for the existence of non-zero perpendicular local currents which we have reported in our paper.We would like to thank one of the referees for suggesting these calculations.
To clarify the role played by structure and the difference in the supercurrent distribution between periodic and quasiperiodic systems, we consider the honeycomb lattice (N = 680) and the Penrose structure (N = 644) under periodic boundary conditions.Figure 11 shows the spatial distribution of J j→i in these structures.Figures 11(a) and (b) show the case of n = 0.3 and θ = 0 on the honeycomb lattice and the Penrose structure, respectively.In the panel (a), one sees that currents are uniformly distributed along the zig-zag lines running parallel to the applied potential.One sees that there is no J j→i in the vertical direction with respect to the applied vector potential.This is expected, due to the translational and inversion symmetries of the honeycomb lattice.As shown in the panel (b), currents J j→i flow non-uniformly on the Penrose structure.Figure 11(c) shows current distribution in the Penrose structure for the case of n = 0.3 and θ = π/10.Flows which are perpendicular to the applied vector potential are shown in red.Such perpendicular currents are thus observed in both the Ammann-Beenker and Penrose structures.The above results show that while perpendicular currents are absent in simple periodic systems such as the square or the honeycomb lattices, they can exist in quasiperiodic structures.We note, finally, that for the honeycomb lattice, J para j→i becomes zero at T = 0, as seen already for the square lattice.As noted in the main text, the existence of non-zero paramagnetic currents at T = 0 is another important qualitative difference between periodic and quasiperiodic systems.

FIG. 2 .
FIG. 2. Real-space distribution (left panels) and the coordination number Zi dependence (right panels) of the local electron density ni [(a) and (c)] and the superconducting order parameter amplitude |∆i| [(b) and (d)] on the Ammann-Beenker structure for n = 0.3 [(a) and (b)] and 0.7 [(c) and (d)] for U = −3, T = 0.01, and θ = 0.For the spatial distribution, we show a part of the system consisting of about 100 sites for visibility.

FIG. 3 .
FIG. 3. Real-space distribution of the local current Jj→i(left), diamagnetic current J dia j→i (middle) and paramagnetic current J para j→i (right) on the Ammann-Beenker structure at U = −3, T = 0.01, and θ = 0.The results were obtained for n = 0.3 (a) and 0.7 (b).Length and orientation of arrows represent the strength and direction of the supercurrent on each bond.Black dots show position of the vertex.In each panel, we show a part of the system consisting of about 100 sites for visibility.

FIG. 5 .
FIG.5.Filling dependence of the local current Jj→i (a), diamagnetic current J dia j→i (b) and paramagnetic current J para j→i (c) for U = −3, T = 0.01, and θ = 0.The results in the panel (a) are separated into the diamagnetic and paramagnetic components respectively in the panels (b) and (c).We note that J para j→i flows in the opposite direction of J dia j→i .The data points for each |ϕn| are slightly shifted in the horizontal direction for the sake of visibility.

FIG. 7 .
FIG. 7. Temperature dependence of the local current Jj→i (a), diamagnetic current J dia j→i (b) and paramagnetic current J para j→i (c) for U = −3, n = 0.5, and θ = 0. Results on a square lattice (SL) of 900 sites are shown by black curves.The vertical dotted line represents Tc of the Ammann-Beenker structure.The results in the panel (a) are separated into the diamagnetic and paramagnetic components respectively in the panels (b) and (c).We note that J para j→i flows in the opposite direction of J dia j→i .The data points for each |ϕn| are slightly shifted in the horizontal direction for the sake of visibility.

FIG. 10 .
FIG. 10.Real-space distribution of the local supercurrentJj→i for n = 0.7 and θ = π 36 .The back-flows are shown as the red downward arrows.We show a part of the system consisting of about 100 sites for visibility.

FIG. 11 .
FIG. 11.Real-space distribution of the local supercurrent Jj→i on the honeycomb lattice (a) and the Penrose structure (b) and (c) at U = −3, T = 0.01, n = 0.3.The panels (a) and (b) are for the case of θ = 0 while the panel (c) is in the case of θ = π/10.The red arrows in the panel (c) show the vertical supercurrents against the applied vector potential.In each panel, we show a part of the system consisting of about 100 sites for visibility.