Electronic instabilities in Penrose quasi-crystals: competition, coexistence and collaboration of order

Quasicrystals lack translational symmetry, but can still exhibit long-ranged order, promoting them to candidates for unconventional physics beyond the paradigm of crystals. Here, we apply a real-space functional renormalization group approach to the prototypical quasicrystalline Penrose tiling Hubbard model treating} competing electronic instabilities in an unbiased, beyond-mean-field fashion. {\color{red} Our work reveals a delicate interplay between charge and spin degrees of freedom in quasicrystals}. Depending on the range of interactions and hopping amplitudes, we unveil a rich phase diagram including antiferromagnetic orderings, charge density waves and subleading, superconducting pairing tendencies. For certain parameter regimes we find a competition of phases, which is also common in crystals, but additionally encounter phases coexisting in a spatially separated fashion and ordering tendencies which mutually collaborate to enhance their strength. We therefore establish that quasicrystalline structures open up a route towards this rich ordering behavior uncommon to crystals and that an unbiased, beyond-mean-field approach is essential to describe this physics of quasicrystals correctly.

Introduction.-The discovery of quasicrystals has triggered exciting, pioneering experimental [1][2][3][4][5][6][7] and theoretical [8][9][10][11][12][13][14][15][16][17] research on the topic.Recently, in experimental studies, superconductivity [18] as well as antiferromagnetic ordering [19] were observed in quasicrystalline systems and their approximants, which are large clusters of quasicrystalline tilings as periodically repeated unit cells.These new experimental findings, especially the reported superconductivity, cannot be explained by our current theory and, therefore, new theoretical approaches are needed.For theory, an intrinsic complication of quasicrystals is the non-layered structure in three dimensions (3D).An exception to this are twisted materials, like twisted bilayer graphene, which form a quasiperiodic lattice in the x − y plane projection for a variety of incommensurable twisting angles [5].Additionally, the lack of translational symmetry in quasicrystals leads to a loss of momentum conservation resulting in severe computational challenges, necessitating the use of simplifying models.A much-studied example, and one of the prototype models of a two-dimensional (2D) quasiperiodic structure is the Penrose model (a Hubbard model on a Penrose tiling) [20,21], which we examine in this letter.It has a single point of global five-fold rotational symmetry and a local ten-fold rotational symmetry [22].The Penrose tiling can be seen as the 2D cut through an icosahedral quasicrystal [1] and constitutes a prototypical model to understand many phenomena in quasicrystalline materials on a qualitative basis, e.g., being cur-rently used to investigate the experimentally observed superconductivity [23][24][25][26][27][28], although a direct connection to the three dimensional materials is less clear.Other platforms which can be used to more faithfully realize the 2D Penrose model are quantum simulators [29] using ultra-cold atomic gases.Possible explanations for the experimentally reported superconductivity in quasicrystals include the existence of unconventional superconducting order generated by spin fluctuations [25].However, the mean-field (MF) theory employed in these studies is biased due to the choice of the decoupling and it remains unclear whether other orderings might prevail.Thus, there is clear demand for an unbiased, beyond-MF study to explore the physics of quasicrystals.
In this letter, we systematically study the electronic instabilities of the quasicrystalline Penrose model employing such a beyond-MF method, the here developed real-space truncated unity FRG (TUFRG) [30,31].The real-space TUFRG offers a versatile tool for the study of translation symmetry-broken models, such as quasicrystals, and scales favorably enough to reach the thermodynamic limit.Utilizing this advance we unveil the surprisingly rich ordering behavior of quasicrystals expanding significantly on what is realized in crystals.Studying different Penrose models and parameters, we find either a mutual suppression of order, similar to a competition of phases known from crystals, or mutual evasion of order as well as mutual collaboration of multiple ordering tendencies.The latter two are both usually not found in crystals.With this we expand the catalog of how phases of matter emerge in quasicrystals and show that in general one requires an unbiased, beyond-MF approach to capture the intricate nature and interplay of orderings in these systems.
Model and method.-We examine a Penrose tiling generated by the substitution method [32] using 10 triangles as initial configuration.Based on this tiling, we construct Hubbard models with sites located on either the vertices or the centers of the rhombi of the lattice, called vertex or center models, respectively (compare Fig. 1).In all simulations presented here, we employ 3126 lattice sites in the vertex model or 3010 in the center model.The phases and critical scales do not change significantly upon further increasing the lattice size.The Hubbard Hamiltonian in second quantization reads with the operators c ( †) i,σ annihilating (creating) an electron on site i with spin σ.We concentrate on two cases for the hopping amplitudes t i,j .First, we apply nearest-neighbor hoppings for which t i,j = t 0 for neighboring sites, shown as red lines in Fig. 1, while for all non-nearest neighbors t i,j = 0. Second, we consider exponentially decaying hoppings using an exponential form t i,j = t 0 e 1− |r i −r j | a .Here, r i is the real-space position of the site i.We choose the minimal distance of any two sites, a = min ij |r i − r j |, as lattice spacing and measure bond lengths relative to it [25].For convenience, we set units by t 0 = 1.We assume spin-independence of the interaction with an onsite repulsion U and a nearest-neighbor repulsion U , resulting in an SU(2)-symmetric Hamiltonian.This is a convenient, but in no means necessary, simplification for our approach (see Appendix A).All of our results are calculated at temperature T = 10 −3 .
To treat possible competing orders in an unbiased fashion, we employ a real-space variant of the TUFRG [30,31,[33][34][35], based on an one-loop formulation of the FRG [36,37].In a FRG scheme, a cutoff function R(Λ) is introduced in the bare propagator, such that the system reduces to a solvable problem at an initial scale.At a final scale, the full solution is recovered.By variations of the cutoff parameter Λ, one obtains an infinite set of differential flow equations.This set of flow equations needs to be truncated in order to be numeri-cally tractable, lending a perturbative motivation to the FRG.Here, we employ one of the most commonly applied truncations, keeping only the effective two-particle interaction without its frequency dependence, which was successfully applied before in 2D for the study of competing instabilities in crystalline structures [36,[38][39][40].The diagrams for the effective interaction can be classified in three groups or channels.Each of the channels is related to a separate effective MF Hamiltonian [36].Thus, a divergence (also called flow to strong coupling) in a certain channel can be directly associated with an emergent order parameter and possibly to a gap opening in the corresponding MF picture, indicative of a phase transition.The emergent phase and its spatial structure are encoded in the type and value of the diverging channel.More specifically, the so-called pairing-channel indicates an opening of the pairing gap, the spin-channel of a magnetic gap and the charge-channel of a charge gap (see Appendix B for details).The real-space TUFRG approach exploits the dependency of each of these channels on so-called "native" indices.These native indices are the dependencies generated in a random phase approximation calculation, which is equivalent to a FRG flow without intra-channel coupling.Dependencies beyond these native indices are only generated at higher orders in the interaction, where the further apart the third and fourth index from the native ones are, the higher the necessary interaction order to generate these contributions gets.Motivated by this, we introduce projections onto each channels' native indices as follows where Γ is a general vertex object and the different projections are marked with a letter to associate them to their respective channel.The pairing-channel is abbreviated as P , the spin-channel as C and the charge-channel as D. The form-factors or bonds f b l (l) form an orthonormal basis on the lattice defined by Due to the pertubatively motivated character of FRG it is reasonable to neglect terms generated at high orders in the interaction, which translates to restricting the bonds used in the expansions in Eq. 2 to a small subset of the full lattice, effectively reducing the size of the projected channels.The flow equations for the projected channels are obtained by inserting the definitions in the original flow equations and introducing a productive unity, see Appendix A for a detailed derivation.Leading to the following flow equations in the SU (2) invariant case where P , C and D are the respective channels projected on their main dependencies and Γ is the effective interaction reconstructed by χpp and χph are the scale differentiated particle-particle and particle-hole propagators which can be calculated from the Greens function G = R(λ)(iω − H) −1 and the single scale propagator where T is the temperature.In this formula we already used the symmetry of the summand w.r.t.frequency to reduce the numerical effort.Throughout this paper we will use the so called Ω-cutoff [41] given by In practice, the spatial ordering is extracted from the leading eigenvectors of the diverging channel.We stop the flow if an eigenvalue of one of the three channels surpasses a threshold corresponding to the stopping scale denoted as Λ = Λ c .Due to the employed open boundary conditions the kinetic energy is reduced at the edges of the lattice, increasing the relative relevance of the interaction there.To avoid phases arising solely due to this boundary effect, in case they do not coexist with a spatially separated bulk phase, we employ a tanh envelope falling off to the boundary: where Introducing this envelop, we can probe the bulk phase-diagram more easily which is related to the thermodynamic limit by the self-similarity of the lattice.
Competing orders in the vertex model.-First, we investigate the vertex model including nearest-neighbor hopping at half filling (µ = 0.0).The model is bipartite and known to show antiferromagnetic ordering in the case of U = 0 and U > 0 [42,43].The density of states (DoS) has a δ-like peak at ω = 0.0 separated by a small gap from the rest of the spectrum (see Appendix D).This peak does not arise due to a Van-Hove singularity (which are very relevant in the context of ordering in crystals), but instead occurs due to a macroscopic number of degenerate states with eigenvalue zero [44][45][46].To access the bulk phase-diagram of the vertex model we employ the envelope Eq. (C2).
We encounter two different phases.For U U , an antiferrromagnetic spin-density wave (SDW) instability prevails with a similar ordering pattern (see lower left plot in Fig. 2) to the ordering pattern at half-filling without nearest-neighbor interactions (U = 0) [42].For increased U /U , the leading divergence changes to a chargedensity wave (CDW) (see upper left plot in Fig. 2).We We employing Eq. ( 7) for the interaction.The upper left plot visualizes the charge order parameter at U = 1.5 and U = 1.The lower left plot visualizes the spin order parameter (or magnetization pattern), at U = 0.5 and U = 2.0.On the right side, the critical scales depending on the nearestneighbor interaction are shown for U = 3.0 in the vicinity of the phase transition, to highlight the suppression of the critical scale upon approaching the phase transition.We find a competition of the two phases at the transition, similar to the behavior prototypically found in conventional crystals.
do not discuss these orderings in detail here, as on the one hand such discussions can be found elsewhere [43,44] and on the other hand we want to focus on the interplay in between the two ordering types.The transition between SDW and CDW is accompanied by a reduction of the critical scales, as can be seen in the right plot in Fig. 2. The critical scale without coupling the different channels with each other is Λ c ≈ 0.4 at the transition, but the ordering vectors are predicted correctly by this simplified calculation.This implies that the divergences are generated within each of the channels individually, in contrast to fluctuation driven divergences which are generated by the feedback of one channel to another.In the simulation incorporating inter-channel feedback, the critical scale is reduced to Λ c ≈ 0.18.Thus, we conclude that the ordering tendencies compete, leading to a mutual suppression of the phases.
To sum up, the physics of the vertex model is similar to the one of crystals with translational invariance, where multiple mutually competing bulk orders take center stage.
Spatial coexistence of order in the nearest neighbor center model.-Next, we examine a center model including nearest-neighbor hopping and interactions without employing Eq. (C2), as here a bulk-boundary coexistence of order is observed.This model is not bipartite and, in contrast to the vertex model, each site has the coordination number four.Its DoS (see Appendix D) displays two main peaks consisting of a macroscopic number of degenerate states (as discussed above): One at ω ≈ 2.35 and one at ω = 2.00.
We concentrate on a Fermi level at the main features at ω = 2.00, which arises in parts due to so-called string states, which are self-similar states with a fractional dimension of 1.44 [9,10].Here, we find a CDW as well as two different SDW divergences depending on the values of U and U .At U U we find an on-site SDW ordering, whereas at U U we find an on-site CDW ordering (see Appendix F).In the transition region between the two, a SDW phase with a more complex ordering pattern emerges.This leads to two distinct phase transitions, one is a smooth interpolation from a bulk ordered on-site SDW to a boundary pinned phase, the second one is a transition between the latter phase and a bulk ordered on-site CDW.The second transition has a mutually evasive nature, as the two occurring orderings have vanishing spatial overlap as we will analyze next.The region in which the transition between SDW and CDW occurs is of special interest to us.In contrast to the vertex model, the strongest two channels, namely charge and spin channel, diverge on equal footing and have main weight in separate regions.This separation is leading to a gap opening in different spatial regions of the lattice, i.e., a coexistence of two different orderings.To show this more clearly, we proceed with a MF decoupling of the effective FRG interaction at the final scale Λ c (see Appendix B or [47]).The spatially separated support of the non-zero gaps for charge and spin order are shown in Fig. 3.The charge order forms in the center as a self-similar structure and shows alternating signs.The spin order emerges close to the boundary of the system and has strongly peaked maximal values separated by an order of magnitude from the lower gap values.In the transition region separating the phases both have low weight.This behavior arises generically for U and U combinations at the transition line between the charge and spin ordering.The coexistence is quite sensitive to changes in U and U , as moving away from the transition line rapidly promotes a single phase.Due to the spatially separated nature, the two orders coexist without affecting each other.
Summarizing these findings, we have identified a coexistence of mutually evasive orders with clear spatial separation in the center Penrose model.Such instabilities can likely only be found in structures with (infinitely) large unit cells such as true quasicrystals, their large-scale approximants, or twisted moiré materials.
Collaborative order in the center model with exponentially decaying hopping amplitudes.-Finally, we investigate the phase-diagram of a center model with hopping amplitudes decaying exponentially in real-space.For such models, recently an unconventional superconducting order has been predicted [25].The DoS in this setup has three main peaks at ω ≈ 0.83, 0.99, 1.23 (see Ap- pendix D).Between the second and the third peak, there is a small energy gap.The most interesting physics is again expected for the Fermi level at the DoS peaks or in their vicinity.For the simulations we chose U = 0 and again employ Eq. (C2) in order to suppress the boundary states.
A SDW is the dominant phase in the whole interaction region at all three peaks of the DoS.Its ordering pattern is dependent on the filling as well as the interaction strength U (see Appendix G).Upon decreasing U to 0.3 we observe the emergence of a sub-leading pairing divergence indicated by a slower decay of the fraction of maximal eigenvalues of the pairing-and spinchannel λ pair max/λ spin max than at higher interactions.This is accompanied by the loss of the standard behavior for an emergent spin ordering, namely that λ charge max /λ spin max = 1 /2 [36].Instead we find that λ charge max /λ spin max decreases slightly upon approaching the critical scale, see Fig. 4.These findings imply only slightly different critical exponents for the three different ordering tendencies.Analogous multi-channel instabilities are known to signal non-trivial non-MF ground states in one-dimensional models for correlated fermions [48][49][50] and were argued to indicate Fermi surface truncations in the two-dimensional Hubbard models [51][52][53].Yet, in this case as well as in ours order collaborate and a proper classification of the true ground state needs further work.In any case, it is exciting to see that quasicrystalline systems offer another playground to exhibit such physics at the frontier of our understanding.
To sum up, in the center model with exponentially decaying hopping amplitudes we find a delicate mutual collaboration of ordering tendencies, pointing to an unconventional ground state.Here MF decouplings are highly biased and a more sophisticated approach needs to be employed.Recent theoretical reports of superconductivity should therefore be taken with caution if they focus on one particular channel in a MF treatment [25].
Conclusions.-We used a newly developed real-space TUFRG formalism to study the electronic instabilities in quasicrystals.With their infinite size unit cells, the spatial degree of freedom in these systems opens up even greater variety of ordering tendencies beyond competition, allowing for spatial coexistence by mutual evasion and joint, collaborative instabilities.We expect similar findings for large but finite unit cell systems, with twisted van der Waals heterostructures as most prominent [54] but certainly not last example.
An intriguing avenue of future research concerns the study of disordered vertex arrangements or quasicrystalline approximants which could help in understanding the interesting interplay between interactions and lattice geometry in quasicrystals.Additionally, an in depth study of the observed orderingss robustness towards rotational symmetry breaking of the lattice could reveal new emerging orderings.As a next step, it is desirable to keep a higher accuracy, e.g. by taking into account the self-energy effects [55] and the frequency dependence of the vertex [31,[56][57][58].This allows to include interactions mediated by phonons or photons and thus to study conventional superconductivity in quasicrystals [24].On the level of real materials the next step is the development of a feasible beyond inter-orbital-bilinears [59] approach, combining the here employed real-space TUFRG with the momentum-space TUFRG [30].With such a scheme 3Dquasicrystalline approximants can be addressed, opening the ground for predictions of phase-diagrams of real materials if combined with density functional theory in-puts [60,61].
Appendix A: Derivation of the flow equations We now sketch the derivation of the flow equations.The effective vertex can be separated in three different channels, each of them related to a specific fermionic bilinear [36].A divergence of their eigenvalues indicates a flow to strong coupling, which in turn can lead to the opening of a specific gap in the self-energy.The P -channel indicates pairing, the C-channel gives information of magnetic ordering tendencies, and charge ordering information is contained in the D-channel.We called the P -channel pairing-channel, the C-channel spin-channel and the D-channel charge-channel in the main body to highlight their physical meaning.Technically the charge ordering information has to be extracted from the physical charge channel [41], but if no other channels diverge, the ordering can be extracted directly from the D-channel or as called in the main body the charge-channel.This can be seen by a mean-field decoupling in the three native channels, see below.In general, we can write the effective two-particle interaction Γ (4) as (suppressing the dependency on the flow parameter Λ) Γ (4) (1, 2; 3, 4) = U (1, 2; 3, 4)+ Particle-Particle (PP) Each number represents a multi-index consisting of orbital, spin and frequency index.We will use quoted numbers, or quoted indices, for all degrees of freedom which are summed.The flow equation for the effective Interaction can be separated in those three channels [30] as can be seen in Eq. ( A2a)-(A2c).As all occurring quantities are Λ-dependent we will from now on leave out the subscript.
dΦ C (1, 2; 3, 4) ))Γ(3 , 2; 3, 2 ), (A2b) ))Γ(4 , 2; 1 , 4).(A2c) We will neglect the frequency dependence of the vertex, therefore we switch to letters as indices, each of which describes a lattice point with associated spin.The truncated unity approach we develop is analog to the momentum-space one.There, the main idea is that without inter-channel coupling the equations amount to a standard RPA series which will only produce a dependency on two spatial indices, or a single momentum (if the interaction is local).If we now include the coupling, we will technically generate dependencies on all indices, but the "native" ones will contain the main features, as configurations beyond those are only generated at higher order in the interaction strength.This indicates that only two orbital indices are of central importance and should be treated with high accuracy and for the others a lower accuracy is sufficient.We want to exploit this by an expansion of each non-native indices in a basis centered at a native one.In a general real-space setting the basis is dependent on the position it is centered around.
Mathematically speaking we use a form-factor expansion of pairwise orthonormal functions which form a complete basis (see Eq. (A3)): Later the included bonds in the unity are restricted to a small subset.We start with general bonds for the derivation and specify them later on.With the usage of a basis we define the projections onto the main dependencies of each channels in Eq. (A4) as This is exact and just a rewriting of the vertex as long as we do not truncate the basis.The full vertex can be recovered by a unprojection, the inverse projection, using the completeness relation in Eq. ( 3).Further we define the channels projected on their respective native indices in Eq. (A5).
Thus we recover the full vertex by applying The flow equations for the projected channels can be derived from Eq. ( A2a)-(A2c) with the help of an insertion of an unity, here shown exemplary for the P -channel The other equations can be derived analogously, resulting in the following coupled set of differential equations where we defined the particle-particle and particle-hole propagator as The equations can be reinterpreted in terms of multi-index blockmatrix-products which are well optimized numerically.
As already mentioned we now want to truncate the basis expansion used in the derivation.This will in fact not change the projected flow equations, but it will change the vertex reconstruction as the inverse projection is not exact anymore.Therefore we obtain projections in between the three channels.The way how the basis is chosen as well as how it is truncated is not unique.The choice of a specific truncation can either be physically or systematically motivated.For example, the first one can be applied in Moiré-lattice models where we have two competing length scales.It is believed that such models can exhibit superconductivity at the scale of unit cell vectors which then could be included explicitly in the expansion.The second approach excludes contributions above a certain bond length or distance in a controlled manner.We can express this truncation as where we defined the set of all allowed bonds L. We define δ L b k to be one if b k ∈ L and zero otherwise.To obtain the explicit projections needed for the flow equations we now insert this into the full vertex projections and keep track of all occurring indices.Without specifying a specific basis or writing out spin dependencies, we obtain the general form of the projections as Before specifying our basis we take care of the spin degrees of freedom.

SU (2) symmetric formulation
To get rid off the spin degrees of freedom we will now assume that we have a SU (2)-symmetric interaction and Hamiltonian, which allows to simplify the channels with the help of the crossing relations [37] Φ P (ik; jl) The spin degrees of freedom can now be eliminated by choosing a specific spin configuration.The full spin dependence can be reconstructed by applying the relations in Eq. (A20).For the sake of simplicity we choose (σ i σ k σ j σ l ) = (↑↓↓↑) and redefine P bibj ij = V P (i, i + b i ; j, j + b j ) and analogously for C and D. The diagrammatic flow equations now read To write out the equations explicitly, we use Kronecker-deltas spanning the lattice around a specific site as a basis, defined as f bi (k) = δ i+bi,k .We additionally assume a density-density type interaction.The sums in Eq. (A11, A12, A13) can be reinterpreted as matrix-products in the multi-indices consisting of orbital and bond index.Thus, the flow equations now simplify to where χpp and χph are redefined in terms of the spin summation.The Matsubara sum can be calculated either analytically which scales like O(N 4 N 4 b )(with N the number of orbitals) or numerically as a summation over the fermionic frequencies which scales proportional to O(N 2 N 2 b N f ) (with N f the number of frequencies and N b the average number of bonds included) with an additional factor of O(N 3 N f ) for the calculation of the Greensfunction, which can be cached and therefore are only calculated once.The number of frequencies can be reduced by employing a non-uniform Matsubara grid to approximately 10 2 −10 3 for reaching an accuracy below 10 −5 at T = 10 −3 , rendering the numerical approach faster than the semi-analytical.The propagators are given as χbi,bj The equation for the D-channel can be simplified by a completion of the square in its flow equations resulting in (defining reducing the imbalance between the three channels greatly.With the before defined Kronecker basis, the projections simplify due to the cancellation of sums.This results in Eq. (A25)-(A27) (compare to [31,33]) introducing the difference set δ d i,j being one if the bond connecting i and j is included in the truncation.

P [Γ]
bi,bj i,j = P bi,bj i,j Appendix B: Decoupling of Vertex Function At the final scale, we obtain an effective interaction for the low lying energy degrees of freedom, with which we can reformulate the Hamiltonian as: where we already went back to the Grassmann notation of the fermionic operators and introduced the multi-indices α = (i, σ i ).A post-FRG mean-field theory can enable a differentiation between competing orders, which are not separated by the FRG.But as the two approaches are decoupled, the resulting gap magnitudes depend on the stopping scale and the results should be therefore only be seen as a qualitative ordering pattern.We now want to derive mean-field equations for the charge and spin gap.We restrict the derivation to a general, SU (2)-symmetric vertex and rewrite the effective interaction in the channel decomposed form: where we neglected the pairing channel as for the cases in which we applied the mean-field decoupling it was suppressed with respect to the other channels.Additionally, we neglect the bare interaction, as in a flow to strong coupling its influence is small and therefore negligible.In order to sum up the spin indices we again apply the relations in Eq. (A20).C and D are the channels we obtain as a result of our FRG scheme.We start with the Grassmann path integral for the partition function: We now introduce a matrix notation for the vertex components without spin degrees of freedom: note that here we still have four spatial indices and no truncation has been applied yet.We will now write out the spin dependence explicitly using the relations from Eq. (A20) In order to have one type of index ordering per channel we commute the Grassmann variables and relabel the indices.The commutation of the fields results in an additional minus sign.We define the fermionic bilinears for each channel as: Note that we have to use the physical channels in order to correctly assign all contributions to the respective mean-field.
Those physical channels are defined as [41] Φ M defines the magnetic-channel, a divergence of it thus describes transitions to a magnetic phase.Φ K describes the charge-channel indicating charge ordering.The charge-channel does not contain a spin divergence part anymore, unlike the D-channel.Thus if both channels, C and D, diverge the charge-channel does only contain the charge divergence which leads to an unambiguous decoupling.Thereby, the full SU (2)-symmetric effective interaction is given as Upon inserting and reordering of the terms we obtain: For the decoupling we now introduce two bosonic fields using a Hubbard-Stratonovitch transformation [36].In general this transformation reads (neglecting all occurring constants as we are only interested in the saddle point): With the restriction that b 2 q = 1 a .We label the bosonic fields φ with the subscript of their respective channel, K for the charge-channel, M for the magnetic-channel.Additionally to the bosonic fields, we introduce the energy gap as condensation term, defined by ∆ I = I • φ I (• denotes a channel specific contraction of the tensor with the two fields, it will be specified later on, the charge gap contains an additional minus sign in its definition).
By introducing the Hubbard-Stratonovitch field, the fermionic part of the action is reduced to a Gaussian and thus integrable.In the absence of magnetic fields we assume that φ M = φ M,↑↑ = −φ M,↓↓ , φ K = φ K,↑↑ = φ K,↓↓ and φ M,↑↓ = φM,↓↑ .Additionally, the diagonal gaps must be self-adjoint.To simplify the notation we introduce the Nambu spinor with which we can rewrite the fermionic action in a compact form (note the reordering of Grassmann variables performed in order to have no minus signs occurring due to the integration): with the Matrix M defined as Here we assumed that the magnetic gap will break the spin rotational invariance.The gaps will be initialized as a superposition of the leading eigenvectors of their respective relevant channel.We now carry out the functional integrals, a Gaussian Grassmann integral, which results in (ignoring the prefactor as it has no effect on the saddle point equation) where we identified the Greens-function G.The mean-field equations are obtained by a saddle point approximation, defined as: The variation of each field needs to vanish individually, leading to a coupled set of self-consistent equations.We sketch the derivation for the charge gap in the following (we suppress spin indices for brevity).
The spin indices are marking the spin block indices of the Greens-function.This equation can now be solved for the charge gap: The equation for ∆ M follows analogously The channels we obtain from the TU calculation have the following index structure: which leads to the following form of the physical channels We now return to the real-space TUFRG notation we introduced earlier, thus the index ordering is changed.Using the eigenvector-decomposition which we obtain as result of our FRG flow we reconstruct the approximate vertex as: This sum will be truncated to a few of the largest eigenvalue/eigenvector combinations indicated by a tilde from now on.Here we still need to sum out the inner spin index as well as the Matsubara sum.As the vertex is frequency independent we can solve the Matsubara sum by complex integration The sum over all leading eigenvectors gives the effective channels resulting in the following set of self-consistent equations If only a single channel is diverging, we can reduce the number of mean-fields to obtain a more efficient description.
In practice we keep the particle number constant during the self-consistency iteration by adding an appropriate, adaptive chemical potential.Additionally, we add and subtract the "zero gap charge gap" and absorb one part into this constant, the other part is used to redefine the charge gap in order to make its value vanish if no divergence is encountered in the physical charge-channel.The fixing of the particle number introduces convergence issues which are addressed by a mixing parameter (see Eq. (B24)).
Additionally, the fixing of the particle number introduces indirect coupling of the two regions, as particles cannot simply be pushed out of the one region.They then flow into the other region creating a charge displacement there.We always performed both simulations, once for fixed particle number and once for non-fixed particle number to check convergence.The non coupled MF/FRG leads to very large gap magnitudes which introduce additional convergence issues.This is resolved by a rescaling of the vertices by a factor, we used 1 /10.As the resulting gap magnitudes are only qualitative this is not introducing a bias.In Fig. 5 we show the gap magnitudes for the non particle number fixed case of the coexistence of magnetic and charge ordering in the center Penrose model.

Appendix C: Interaction envelop
For the interaction envelop we apply as scaling factor for the interaction (d i is the distance of site i to the center of the tiling and d max = max i d i ).For the nearest-neighbor interactions, this factor is applied too.To keep the interaction symmetric and C 5 -invariant we need to symmetrize the scaling factor (see Eq. (C2)), thereby we introduce a slight bias as explained below.
This creates a region of higher U /U-ratio than in the bulk.This slight deviation seems not to have any influences as each interaction term is lower individually.The effective lattice size is of course reduced due to the application of the screening.The distance dependence of the screened on-site and nearest-neighbor interaction is shown in Fig. 6.For pure interaction, the main results are summarized in Fig. 9.We find a charge-density wave phase at dominant nearest-neighbor interactions (U ) and a spin-density wave phase at dominant on-site interactions.The latter is expected as in the limit of weak U we should recover the known antiferromagnet for the vertex model.The SDW phase has its main weight at the boundaries, which is a consequence of the open boundary conditions and can be understood as follows.At the boundaries the local kinetic energy bandwidth is reduced, which results in a larger U /W D ratio, which in term leads to faster divergences of diagrammatic re-summations at the boundaries.The sub-leading eigenvectors of the SDW divergence also show this bulk behavior.Upon decreasing U , this effect is also vanishing as it should be.The CDW phase loses weight towards the boundaries and thus is seen as a bulk phase.Due to the reduced boundary weight, it could be possible that a coexisting boundary phase emerges at lower critical scales or lower temperatures.However, we did not find a phase coexistence.At the transition between CDW and SDW we find a mutual suppression of the phases resulting from their non-zero overlap in space.This behavior can be understood as follows.In a spin-channel (SDW) antiferromagnetic divergence the chargechannel (CDW) has an eigenvalue, describing the strength of the order, half as large as the one of the spin-channel, which already follows from the RPA.This eigenvalue is in our sign convention positive.The diverging eigenvalue of the CDW-RPA itself is negative, thus the two will suppress each other in the vertex reconstruction.This partial cancellation is the source of the reduced critical scales and thus reduced transition temperatures.As soon as the CDW divergence is dominant the critical scales grow rapidly again, as shown in the right inset of the left plot in Fig. 9.

FIG. 1 .
FIG. 1. Illustration of quasicrystal lattices.Four times iterated Penrose tiling with the vertex model left and the center model right.Nearest neighbors are marked as red bonds.The rhombi referred to in the main text are shown in grey for the center model.The lattice sites are marked as blue dots.The hoppings in the Hamiltonian are set to t0 only if a red line connects the two sites in question.

FIG. 2 .
FIG.2.Vertex model at µ = 0.0: Competition of order.We employing Eq. (7) for the interaction.The upper left plot visualizes the charge order parameter at U = 1.5 and U = 1.The lower left plot visualizes the spin order parameter (or magnetization pattern), at U = 0.5 and U = 2.0.On the right side, the critical scales depending on the nearestneighbor interaction are shown for U = 3.0 in the vicinity of the phase transition, to highlight the suppression of the critical scale upon approaching the phase transition.We find a competition of the two phases at the transition, similar to the behavior prototypically found in conventional crystals.

FIG. 3 .
FIG.3.Center model at µ = 2.0: Coexistence of order.Gap values of the decoupled effective interaction at Λc for U = 1 and U = 1.015.Leading ordering strengths of the charge and spin ordering channels have approximately the same magnitude.The calculated gap magnitudes are shown in a red-blue color scheme for the spin gap, which is equivalent to the magnetisation.In the violet-lime color scheme the charge gap is shown.The two faces have no spatial overlap of magnitudes above 10 −3 max(∆).The coexistence of the two phases is found to be a general feature of this phase transition.

FIG. 4 .
FIG. 4.Center model with exponentially decaying hopping amplitudes at µ = 0.99: Collaboration of order.In the upper left, the maximal eigenvalues of the pairingand charge-channel relative to the one of the spin-channel are shown for U = 0.3.The expected ordering pattern for the spin-channel is visualized in the upper right plot, the one for the charge-channel in the lower left plot and the one for the pairing channel in the lower right plot.We observe a slow decay of the relative maximal eigenvalues pointing at a possible unconventional groundstate.The ordering patterns of the three channels have spatial overlap, therefore there is no spatial evasion of order.

FIG. 5 .FIG. 6 .
FIG.5.Gap magnitudes of SDW and CDW in the center Penrose model at µ = 2.0, U = 1.0,U = 1.015 and T = 10 −3 without fixed particle number.It compares well to the gap predicted by the fixed particle number calculation.he calculated gap magnitudes are shown in a red-blue color scheme for the spin gap, which is equivalent to the magnetisation.In the violet-lime color scheme the charge gap is shown.The two faces have no spatial overlap of magnitudes above 10 −3 max(∆).

FIG. 7 .
FIG.7.Non interacting DoS for the three types of Penrose models used.Left, we have the vertex-type Penrose model with 21106 sites with the δ-like peak at µ = 0.In the middle, we show the center-type model with hoppings defined using the graph, the main peak is at µ = 2.0.On the right, we show the center-type model with exponential decaying hoppings, with three main peaks.In both center type models the lattice contains 20800 sites.

FIG. 8 .
FIG. 8. Phasediagram of the vertex Penrose model employing Eq.(C2).we find a SDW (star marker) and a CDW (triangular marker) phase.The critical scales are given by the color-code.

FIG. 9 .
FIG. 9.The left plot shows in the upper left inset the charge order parameter at U = 1.5 and U = 1.The lower left inset visualizes the spin order parameter (or magnetization pattern), at U = 0.5 and U = 2.0.On the right side, the critical scales depending on the nearest-neighbor interaction are shown for U = 3.0 in the vicinity of the phase transition, to highlight the suppression of the critical scale upon approaching the phase transition.The right plot shows the phase diagram of the vertex model with a SDW (star marker) and a CDW (triangular marker) phase, the phase diagrams for the pure case and the one employing an envelop differ only very slightly.

FIG. 10 .
FIG. 10.Phase-diagram of the center model at different chemical potentials, we have µ = 1.9 in the upper left, µ = 2.0 in the upper middle, µ = 2.1 in the upper right, µ = 2.3 in the lower left, µ = 2.35 in the lower middle and µ = 2.4 in the lower right plot.U and U are varied at T = 10 −3 , nearest neighbors are included in the calculation.

FIG. 11 .
FIG. 11.Ordering predictions for each some chosen data-points at each chemical potential, where we show µ = 1.9 in the upper left, µ = 2.0 in the upper middle, µ = 2.1 in the upper right, µ = 2.3 in the lower left, µ = 2.35 in the lower middle and µ = 2.4 in the lower right plot, from FRG plus post-production mean field theory at chosen data points.U and U are varied at T = 10 −3 , nearest neighbors are included in the calculation.