Barriers to macroscopic superfluidity and insulation in a 2D Aubry–André model

We study the ground state phases of interacting bosons in the presence of a 2D Aubry–André (AA) potential. By using a mean-field percolation analysis, we focus on several superlattice and quasicrystalline regimes of the 2D AA model, including generalisations that account for a tilting or skewing of the potential. We show that barriers to the onset of macroscopic phases naturally arise from weakly modulated domains in the 2D AA model. This leads to the formation of extended crossover domains, in which the macroscopic properties are dominated by a minority of the system. The phase diagrams then exhibit substantially different features when compared against crystalline systems, including a lobe-like or wave-like appearance of the Bose glass, sharp extrusions and thin, percolating clusters. By studying the 2D AA model across multiple regimes, we have shown that these extended crossover domains are not distinct to a small set of parameters.


Introduction
The Bose Glass (BG) is a famous example of a quantum phase that stabilises a coexistence of insulating and SuperFluid (SF) ground states on local scales [1].While the BG is viewed as macroscopically insulating, local SF domains will mean that the phase is compressible, unlike a Mott-Insulator (MI).Furthermore, these SF domains do not percolate, leading to the absence of macroscopic phase coherence [2][3][4][5].The BG is perhaps most commonly associated to disordered, crystalline systems.In this case, short-range random disorder is known to introduce atomic localisation for both interacting [6][7][8][9] and non-interacting systems [10,11].The disorder averaged phase is of particular interest in these models, as no single disorder realisation should dominate the overall physics.In other words, the local fluctuations do not play an important role in the underlying phase transitions.For this study, however, we wish to study the properties of models that are not random.In this scenario, local variations of any order parameters are far more important than they would be for disordered systems.This can allow for the formation of mixed phases, where the macroscopic properties of a phase are dominated by a minority of the system.
The system we study is based on a 2D extension of the quasiperiodic Aubry-André (AA) model with on-site interactions.The AA model and its generalisations have been widely studied in the single-particle picture [12], and are known to host self-dualities [13][14][15], novel dynamical properties [16][17][18] and topological phases [19][20][21].While many studies have been conducted in 1D, various 2D extensions of the AA model have also found similar properties [22][23][24].For quasiperiodic systems, we have a distinct scenario in which we have both short-range disorder and long-range order present in the system.Importantly, however, this short range disorder is not random, and the systems structure is normally characterised by long-range self-similarity.Furthermore, this order can naturally lead to the formation of barriers to macroscopic phases.
Barriers to macroscopic phases manifest as weakly modulated domains in the 2D AA potential, which can form in both quasiperiodic and periodic limits.These domains can then stabilise the mixed phases on the lattice.This includes regimes where a small SF domain of a few sites can percolate through the system and support macroscopic superfluidity, despite a majority of the phase possessing insulating characteristics.We can also observe the opposite scenario, where small, percolating domains with insulating behaviour will block the onset of macroscopic superfluidity.
Recently, the BG has been theoretically observed for some 2D quasicrystalline systems through a mean-field analysis of tight-binding models [25] and quantum Monte Carlo studies of continuous systems [26].To-date, the mixed phases have not been studied in detail, which is the central aim of this work.We will explore the local nature of ground state phases for the many-body 2D AA model.We achieve this through a Gutzwiller mean-field analysis, based on percolation methods for inhomogeneous systems [3-5, 25, 27-30].
Here, we will present our results for the 2D AA model as follows.First, in Sec. 2, we will define the 2D AA model that we work with in this study, including a discussion on the mean-field percolation methods and order parameters used to determine different macroscopic and mixed phases.Using these definitions, we then discuss the structure of mixed phases in Sec. 4, and show several examples of the mixed phases on the lattice for different forms of the underlying potential.Finally, in Sec. 5, we present full ground state phase diagrams for the system in different parameter regimes, before ending with our conclusions in Sec. 6.

2D Aubry-André Bose-Hubbard Model
In the following, we will consider an inhomogeneous Bose-Hubbard model on a 2D square lattice with unit spacing and Hamiltonian where N is the total number of lattice sites, U is the on-site interaction strength, i is an on-site energy, J is the tunnelling coefficient, i, j denotes the sum over nearestneighbours, µ is the chemical potential, bi ( b † i ) are the bosonic annihilation(creation) operators and ni is the number operator.We will consider the i to be distributed according to a 2D AA quasiperiodic potential of with λ denoting the modulation strength, β 1,2 are the wavenumbers and x, y are the 2D spatial coordinates of the lattice.Throughout this work, we will consider a N = 99 × 99 lattice.
If the wavenumbers in Eq. ( 2) are irrational, they will be incommensurate with the lattice spacing.The distribution of i will then be quasicrystalline, and contain longrange order.On the other hand, if the wavenumbers are rational, we will instead have a commensurate, crystalline distribution of i on the lattice.Throughout this study, we will refer to the potentials with irrational wavenumbers as quasicrystalline distributions and rational wavenumbers as superlattice distributions.In Fig. 1, we plot visualisations of the 2D AA potential for the three different regimes used in this study.Starting with equal wavenumbers in Figs.1(a,d), we can see the appearance of weakly modulated lines with i = 0, as observed in previous studies [23,25].These lines are predictable from the form of Eq. ( 2) as where k is an odd integer, β is a fixed wavenumber when β 1 = β 2 and Y can be either the x or y coordinate.If d(Y ) is sufficiently close to zero, then the corresponding row/column of sites will be weakly modulated.In the single particle picture, it is known that these weakly modulated lines can destabilise the mobility edge that is typically seen in AA models, and can even support ballistic transport [23].
In this work, we will also study other quasicrystalline distributions of the potential that, in general, no longer stabilise precise lines of weak modulation that match the geometry of the underlying lattice.Of particular relevance to experimental protocols is the consideration of tilted quasicrystalline potentials.For example, this would amount to the rotation of some bichromatic quasiperiodic potential on top of an optical lattice, which could describe alignment errors.We will consider a rotation of the 2D AA potential, which transforms the spatial coordinates of Eq. (2) according to where θ denotes an anti-clockwise rotation angle from the x-axis.We show an example of a tilted 2D AA potential in Figs.1(b,e).Finally, we will also consider the case of unequal wavenumbers.There is of course many different choices one can make for the wavenumbers, and we shall consider a particular set in order to demonstrate the rich physics that is present in the manybody 2D AA model.We will therefore consider the following parametrisation for the wavenumbers where φ is an effective skew "angle" between β 1 and β 2 .This ensures at least one irrational wavenumber for φ > 0 • .Given the reflectional symmetry about φ = 45 • , we therefore consider a range of φ between 45 • and 0 • .A skewed 2D AA potential is depicted in Figs.1(c,f), which shows a similar structure to that of the tilted potentials.

Gutzwiller Mean-Field
In order to study the many-body properties of the 2D AA model, we consider a meanfield percolation analysis.First, we decouple correlators according to a Gutzwiller ansatz where ϕ i = bi is the mean-field order parameter at site i, which is taken to be real without loss of generality for the considered model.By substituting the above relation into Eq.( 1), the Hamiltonian for the inhomogeneous Bose-Hubbard model can be written in the local number basis for each site as The ground state for the entire system can then be found by diagonalising Eq. ( 7) for each site and converging.During this process, the order parameters are updated according to the local number basis where z is the maximum number of particles per site and f n are elements of the lowest energy eigenvector for site i.Given an initial set of order parameters ϕ i , this process can be repeated in a self-consistent manner, until convergence to the true ground state.We can also calculate the local density for each site As a practical note, z should in principle be infinite for bosonic systems, but it is set to a finite value for numerical purposes.The ground state will converge with increasing z.We take z = 10 for the results presented here, which is well into the region of convergence.

Percolation and Mixed Phases
The simplest case of the homogeneous Bose-Hubbard model with on-site interactions normally permits the existence of two unique ground state phases; the MI and SF.In the context of mean-field theory, the presence of this macroscopic order can be based on whether the average order parameter φ is finite.This can be captured with the correlation fraction F of a state, which we define as where N ϕ is the total number of sites with finite ϕ i and N is the total number of sites.
In other words, F = 0 for the MI and F ∼ 1 for the SF.
If the model has disorder, then a BG is expected to prevent a direct MI to SF transition [55].Unlike the SF, a BG no longer supports macroscopic phase coherence, and is hence also an insulating phase.However, the BG will have small, isolated SF domains, which means that the phase is compressible, unlike the MI [55].We can identify a BG-SF transition based on the percolation of sites with a local SF character, as has been done in prior works on percolation based methods [5,[27][28][29].To do this, we define a percolation probability P as where N span is the number of sites in a percolating cluster.For the MI and BG, we now have P = 0, whereas for the SF P > 0. The application of a percolation analysis with conventional mean-field approaches has been shown to produce comparative results to those obtained by quantum Monte-Carlo for disordered systems on a square lattice [3,5,56,57].
For the 2D AA model, there is the possibility that short-range, non-random disorder can result in phases that are dominated by local properties of the system.An example of this could be a thin line of SF sites percolating across an otherwise insulating system.This state, while being for the most part insulating, could be considered to host a macroscopic SF.The macroscopic properties of the ground state are then dominated by the local properties of comparatively few sites, e.g. for the N = 99 × 99 system we consider, a single percolating line of SF sites could be supported by as little as 1% of the sites.We will refer to ground states dominated by local properties as mixed phases.
We will characterise mixed phases by considering two local discrete functions that measure the phase at each site and its surroundings.These will be referred to as locality functions.First, we define the insulating locality function S MI i , which is 1 if a MI site (integer density and zero ϕ i ) is surrounded by MI sites across each of its bonds, and 0 in all other scenarios.The complimentary function 1 − S MI i is therefore 1 if a MI site has at least one neighbouring SF site.Similarly, we can also define the SF locality function S SF i , which is 1 if a SF site (non-zero ϕ i ) is surrounded by SF sites across each of its bonds, and 0 otherwise.The complimentary function 1 − S SF i is then 1 if a site has at least one neighbouring MI site.By looking at the average values of these distributions, S MI and S SF , we end up with a measure between 0 and 1 which is more sensitive to the local structure and fluctuations of clusters.
From this, we can first define a weak SuperFluid (wSF) phase, which has a majority of sites being insulating, but has a small percolating SF domain.This thin SF cluster acts as a barrier to macroscopic insulation within the system.In other words, more than 50% of sites will still have a MI character.Note, at the 50% threshold of MI to SF sites, S MI and S SF can change significantly, depending on the local distribution of clusters.For example, if the system is composed of two large clusters of MI and SF character on the square lattice, the maximum values of the locality functions will be which is 0.4899 for a N = 99 × 99 lattice.On the other hand, if there is checkerboard pattern of MI and SF sites, each site will be surrounded by neighbours with different phases, resulting in S MI = S SF = 0.A similar result also follows if the rows/columns on the lattice possess an oscillating MI/SF character.The wSF ground state phase supports macroscopic superfluidity, and will therefore be characterised by S SF ≈ 0, S MI ≥ 0, and P > 0.
The opposite scenario is that of a weak Bose Glass (wBG), which has a majority of sites being in the SF phase, but has small MI domains which prevents a SF percolating through the full system, meaning the state lacks macroscopic phase coherence.This As expected, N V is minimised around rational β, but quickly inflates as we approach quasicrystalline distributions.N L , on the other hand, can either be peaked or minimised to zero around rational β.However, N L will remain largely finite across a large range of irrational β.
ground state will then be characterised by S SF ≥ 0, S MI ≈ 0 and P = 0. Finally, we note that for the macroscopic SF and MI phases, the corresponding locality functions will converge towards S SF = 1 and S MI = 1 respectively.

Weak Modulation Lines
From Eq. (3), we can expect weakly modulated domains to frequently appear in the 2D AA model for a large range of wavenumbers β 1,2 .These domains can act as barriers to macroscopic percolation, and will therefore influence the formation of mixed phases.In order to show that this is indeed the case, we will consider the relation between the wavenumbers and the total number of weakly modulated lines.First, we take a weakly modulated line to be defined when the difference function d(Y ) < 10 −2 for a given row/column, and will consider the limit of β 1 = β 2 ≡ β.For each row/column of the potential, we then count the total number of weakly modulated lines N L that fall below this threshold as a function of β for a 99 × 99 lattice, shown by the line plot in Fig. 2. Due to the quasiperiodic nature of the potential, N L will rapidly oscillate with small changes in β.Importantly, however, there are many regions at irrational β in which N L is finite.We also colour regions of the plot according to the number of unique values N V (to 3 significant figures) in the 2D AA potential.As expected, the number of unique values will be minimised around rational wavenumbers, such as β = [2/3, 3/4, 4/5 ...], and the resulting potential will take a superlattice form.

Example Mixed Phases
In order to better understand the influence of mixed phases, we will now consider the structure of wSF and wBG phases on the lattice.First, we will discuss the case of equal wavenumbers β 1 = β 2 ≡ β, which results in the stabilisation of weak modulation lines.
For this case, we can observe both the wSF and wBG phase in Figs. 3 and 4, where we have plotted the local density ρ i and order parameter ϕ i .In Fig. 3, we can clearly see the presence of a wSF phase, in which we have narrow percolating lines of SF states with non-zero ϕ.The wBG phase is a little more subtle, but it can be seen in Fig. 4 that it is possible to have large numbers of sites with non-zero ϕ that do not percolate  through the full system.Fig. 4(d) is a particularly clear example of this.
Next, we now turn to the more general cases of the tilted and skewed potentials.As we have shown in Fig. 1, it is also possible for both the tilted and skewed potentials to exhibit similar lines of weak modulation.These domains will of course no longer form as precise vertical or horizontal lines on the square lattice, but rather as correlated patches of weak modulation.The strength of modulations across these patches will be more significant than what is typically observed on the weak modulation lines, which will therefore affect the stability and structure of mixed phases.In Fig. 5, we show examples of the mixed phases for the tilted potentials.Here, we can clearly observe the formation of a wSF phase across diagonal lines of weak modulation in Fig. 5(e).Furthermore, the wBG can also form for certain regimes in Figs.5(d) and (f).These mixed phases possess smaller clusters of MI sites across specific diagonal/zig-zag patterns on the lattice, which will again prevent the onset of SF percolation.We also consider several example phases for the skewed potentials in Fig. 6.For these potentials, the wSF can sometimes appear in a more trivial manner across diagonal regions, with some examples in Figs.6(d) and (e).The wSF phase in Fig. 6(e) is a special case in which SF clusters only form across diagonal lines, giving rise to a 1D structure of the order parameters.Finally, the wBG can also form with interesting local properties and patterns of MI clusters, with an example in Fig. 6(f).Mixed Phases In Fig. 7, we plot phase diagrams for the system as function of J/U and µ/U for β = 1/ √ 2, with different λ/U .We segment the ground state phase into SF, MI, and BG domains, as has been previously considered.However, we can now confirm the presence of wSF and wBG domains for the 2D AA model.The average of the discrete locality functions, S MI and S SF , are shown in Fig. 8 for the phase diagrams of Fig. 7.We also plot additional boundaries on these phase diagrams, with the red line indicating when SF sites account for more than 50% of the overall phase, and the black line for the onset of SF percolation.The average locality functions show that there are domains where S MI and S SF are non-zero but not unity, meaning that all sites are not yet surrounded by the same local phase.From these average locality functions, we can define the wSF and wBG phases as previously discussed.The presence of the wSF and wBG phases is far clearer in the case of larger λ/U , shown by the extended transition in the locality functions in Fig. 8(b,d).This shows a direct relation between the strength of the 2D AA potential and the presence of the mixed phases, which is highlighted further by the growth of the wSF and wBG domains in Fig. 7 for increasing λ/U .
A particularly noteworthy feature in Fig. 7 is the presence of lobe-like BG domains, with peculiar extruding features of a wBG character.In Ref. [25], these features were speculated to form due to the presence of weakly modulated lines.We can now confirm -through the measurement of the wBG -that this is indeed the case.The sharp protruding features of the wBG are due to the support of an insulating phase from the weakly modulated lines, and such a structure would most likely be destroyed in a fully disordered system.We also find that the wSF is stabilised by the weak modulation lines near integer µ/U .This is due to small SF clusters forming and percolating throughout the system, but with support from relatively few lattice sites.

Mixed Phases and Critical Points
Due to the inhomogeneous nature of the 2D AA model, calculating the full ground state phase diagrams over large ranges of β would be inefficient.For the sections that follow, we will instead focus on the properties which will be most β dependent, i.e. the mixed phases.Varying β in general can only impact two of the properties of the system: (i) β can be tuned to a quantity commensurate with the underlying lattice, therefore, away from being a quasicrystal, or (ii) β can alter the appearance/location of the weak modulation lines that support the mixed phases.In order to study the effects of both scenarios, we will consider the behaviour of the BG-SF critical points at different µ/U for varying β.We choose two critical points, which correspond to the lobe-like behaviour of the BG and stabilisation of a wSF at an integer µ/U = 1, and the extrusion of the wBG at µ/U = 0.39.It is at these points of the ground state phase diagram that the presence of mixed phases is particularly pronounced.In Appendix A, we also consider full phase diagrams towards a flat, commensurate limit of β = 1.In Fig. 9, we plot the BG-SF critical points as a function of β.Starting with the red line, where µ/U = 1, we can immediately see that the critical points show little variation in this interval, with most J C /U ∼ 10 −3 .The largest fluctuations occur for the superlattice limits of the 2D AA potential, which destabilise the lobe-like structure of the BG.This behaviour can be linked to the underlying structure of mixed phases, as observed in Fig. 3.For quasicrystalline distributions with irrational β, the local SF clusters account for a very small fraction of the overall phase.Despite this small fraction, the local SF clusters still percolate due to their formation across weakly modulated lines.This is in contrast to what is seen for superlattice potentials with rational β, where the SF clusters generally account for a much larger fraction.The structure of local SF clusters can take a more crystalline form in these limits, and may even introduce larger fluctuations to the weakly modulated lines.This short-range order will influence the  showing the J C /U BG-SF critical point.For certain β, we also show the rational fraction above J C /U peaks and on the β grids as vertical lines.The structure of local SF clusters can be divided into two distinct regimes based on whether or not J C /U is maximised or minimised for each µ/U .These then directly correspond to superlattice and quasicrystalline realisations of the AA potential.
percolation of local SF clusters, and therefore destabilise the lobe-like structure observed on the phase diagrams.
We also plot the BG-SF critical points for µ/U = 0.39, which is the blue line in Fig. 9, corresponding to the extruding feature for the wBG.Here, we see that this particular feature is quite sensitive to smaller changes in β.The positions of the local minima and maxima of J C /U shows several analogous properties to what was observed previously.For quasicrystalline β, we know that the phase will possess many SF sites, as per Fig. 4. Small MI domains across weakly modulated regions will, however, prevent percolation until a large J C /U .When considering superlattice potentials, the onset of short-range order can remove such regions in the system, leading to a reduction in the J C /U critical point.
Depending on the rational fraction of β in Fig. 9, we can also link the structure of mixed phases to the range of short-range order.Generally speaking, for the smaller numerators/denominators, the critical point will be significantly shifted from what is seen at irrational β.The reason for this is due to the absence of long-range variations in the on-site potential.As the numerators/denominators are enlarged towards irrational limits, we introduce long-range variations to the on-site energies i , and hence we observe less dramatic shifts in critical behaviour.These properties effectively mimic what is seen in Fig. 2 for the number of weakly modulated lines and unique values.When β is irrational, long-range variations in energy can allow for the appearance of very weakly modulated lines in Eq. ( 3) at specific rows and columns.If β is rational, then it is possible that the condition in Eq. ( 3) may not be satisfied for any row or column, leading to the absence of barriers to macroscopic percolation.For superlattice realisations of the potential with smaller fluctuations in energy, wSF and wBG phases will appear less frequently on the phase diagrams.Furthermore, it is also possible to stabilise macroscopic DW order at finite J/U .

Density Waves from Superlattice Potentials
We now focus on the case of superlattice potentials when β is rational, which is an interesting limit of the 2D AA model.We plot the ground state phase diagrams for a number of different rational wavenumbers in Fig. 10.With the superlattice potentials, it is still possible to observe wSF and wBG phases in certain regimes, but their domain is severely reduced compared to the case of quasicrystalline potentials.
Interestingly, for some of the superlattice potentials, we can also observe the formation of macroscopic DW phases at finite tunnelling strengths, as shown in Fig. 10.Similar to the MI, the DWs are macroscopically insulating and possess no transport.They will, however, contain non-uniform densities.Each distinct DW lobe therefore corresponds to a different non-integer filling of the lattice.The size and structure of DW lobes and intermediate wBG/wSF phases depends on the rational wavenumber and extent of long-range fluctuations.If the numerator/denominator of β is small, there will only be a few distinct on-site energies in the 2D AA potential.The number of distinct DW lobes that appears on the phase diagram will then also be small, as seen in Figs.10(b,c).Furthermore, the width of these DW states is comparable to the size of the MI lobes in J/U .For the larger numerators/denominators of β in Figs.10(a,d), fluctuations in on-site energies will become more pronounced, and the number of unique DW states will increase.The width of each DW lobe in J/U will decrease, however, indicating a stronger sensitivity of these states to particle number fluctuations.
To better illustrate the stability of macroscopic DW order in the 2D AA model, we plot the average density in Fig. 11, as a function of µ/U around a rational wavenumber of β = 3/4, when J/U = 0 and λ/U = 0.15.Here, we see that when β is varied between rational and irrational limits, the properties of the DW states begin to drastically change.Notably, we observe that long plateaus of average density no longer form in the DW Average density as a function of µ/U around β = 3/4, when J/U = 0 and λ/U = 0.15.All states are insulating for these parameters.We observe that the average density plateaus in the DW domains become unstable as we move further away from the rational limit in β.This is a consequence of the number of unique values in the energy potential increasing.
regions when β starts to become incommensurate with the lattice.In other words, the total number of distinct DW states has significantly increased, but their width in µ/U is vanishingly small, giving rise to an almost continuous behaviour of ρ at J/U = 0.It is important to note that away from the rational limit of β = 3/4, the number of unique values in the potential will quickly inflate, as seen in Fig. 2. The small width of each DW in µ/U implies that the system will be very sensitive to particle number fluctuations, and hence why even small tunnelling strengths can immediately destroy in these regimes.Similar properties are also around other rational β.

Ground State Phase Diagrams
In Fig. 12, we plot phase diagrams for two modulation strengths λ/U over a range of rotation angles, with inset figures showing the tilted potential.Starting with small rotation angles in Figs. 12, we observe that the tilted potential no longer stabilises weak modulation lines, but rather weakly modulated zig-zag patterns on the lattice.The overall structure of the phase diagrams is now far more reminiscent of those observed in randomly disordered systems, but now with the underlying phases possessing long-range order.We do observe regions where the ground state is a wSF or wBG for larger λ/U , but this is in relatively narrow regions for small θ.This is evidence that weak modulation regions and mixed phases are playing a smaller role, but can still dominate the ground state properties in specific regions.
As θ increases, the AA potential will generally have less correlated regions of weak modulation, exaggerating the loss in observable effects and mixed phases.A special case is seen for θ = 15 • , where the extruding feature of the BG reappears on the phase diagram.By inspecting the AA potential in Fig. 12(c), it can be seen that there are now diagonal lines of weak modulation, which will act as a barrier to superfluid percolation.SF clusters can also form on these diagonal lines, as shown in Figs.5(b,e).The SF or wSF phase is not stabilised around integer µ/U , however, since the diagonal lines are not connected by bonds, and hence no percolation can occur through them.This again leads to the destruction of the lobe-like BG pattern on the phase diagram.By further increasing θ to 30 • in Figs.12(d,i), we observe similar properties to before, with the mixed phases no longer dominating large regions of the phase diagram.
For larger angles, e.g.θ = 40 • , the mixed phases are actually enhanced as shown in Fig. 12(e,j).This is due to the large rotation angles tending the potential towards a more uniform limit, with large regions of positive/negative on-site potential connected by extended flat regions; see the insert of Fig. 12(e).We also note that for all considered phase diagrams here, no DW phases are found to be stable at finite tunnelling strengths.
From these results, we can see that properties of the mixed phases are strongly dependent on large, correlated domains of weak modulation.For θ = 0 • , weakly modulated lines have a significant impact on the percolation of local SF clusters.When considering tilted potentials, we generally can no longer form regions of weak modulation across the entirety of the lattice.Even if smaller domains of weak modulation exist, these are connected by paths of stronger modulation, as is the case for the zig-zag patterns in Fig. 12(b).For smaller rotations, we can still observe differences in the formation of mixed phases, which is a consequence of these weakly modulated patterns influencing 0.06 0 0 45 local SF percolation.

Mixed Phases and Critical Points
As we have seen, the study of J/U critical points at µ/U = 1 and µ/U = 0.39 have been very important measures in the characterisation of mixed phases.In Fig. 13, we again plot the BG-SF critical points in J/U for a range of tilt angles θ.Here, we immediately observe several profound differences to what has been seen previously.First, the critical point at µ/U = 1 quickly rises for small rotations and oscillates around J/U = 0.008 for the majority of θ, which implies that the lobe-like structure of the BG will be partially destroyed.This is to be expected from the previous phase diagrams.Compared to the critical points with no rotations, there are no significant fluctuations in J C /U for the majority of tilt angles.
As we approach θ = 45 • , the critical point will slowly decrease towards 0, indicating the return of weakly modulated domains.At θ = 45 • , it can be seen that there is a sudden jump in the critical points.The reason for this can be inferred from the quasicrystalline distribution.Taking θ = 45 • , the 2D AA potential reduces to flat distribution.The critical points at θ = • are therefore those of the disorder free system at an effective chemical potential of µ − 2λ.At µ/U = 0.39, the BG-SF critical points also show significant differences.As before, the extruding feature effectively vanishes, and the critical point will remain stable at J/U = 0.04 for most of the tilt angles.Curiously, however, at θ = 15 • , the extruding feature does return for a small range of tilt angles.The reason for this change is due to the existence of diagonal weakly modulated lines, as shown in Fig. 12(c).As we approach the limit of θ = 45 • , the critical point slowly increases, before another sharp change when the 2D AA potential becomes more uniform.The tilted potentials substantially impact the formation and structure of mixed phases in the 2D AA model, but importantly do not completely destroy them.(b,g), the AA potential forms similar looking zig-zag patterns of weak modulation.This will generally destabilise the lobe-like structure of the BG, and remove the extruding features.As before, the wSF and wBG phases will become less pronounced, highlighting changes to the percolation of local SF/MI clusters.However, we do observe that the wSF and wBG phases are more robust for stronger AA potential strengths.This is also reflected in the structure of mixed phases, as seen in Figs.6(a,d).
For φ = 30 • in Figs.14(c,h), we have a special case that has one rational and one irrational wavenumber, leading to a potential with lines of constant i along the diagonals.Different clusters of phases will, therefore, only form across these diagonals, as per Figs.6(b,e).This results in a phase diagram that possesses quite a different structure from what has been seen previously.For instance, the BG now has multiple smaller lobes at certain µ/U .The wSF shows even more profound differences due to the unique characteristics of this potential.Effectively, in this scenario, the wSF is a measure of a 1D MI-SF transition in the 2D structure.Furthermore, we also note that the number of unique values in the potential is comparable to the lattice side length 0.06 0 0 45 due to the diagonal lines of constant i , leading to the appearance of DW states at finite tunnelling strengths around J/U ≈ 10 −2 .At several of the intermediate skew angles, we also see phase diagrams comparable to Figs. 12(d,i) and randomly disordered systems.Interestingly, the wSF can be stabilised across a larger set of chemical potentials, away from the BG domains.Finally, as φ → 0 • , the potential will again become more crystalline and uniform, similar to what was seen in Figs.12(e,j) but now with a diagonal, effectively 1D, form.

Mixed Phases and Critical Points
Here, we will again plot how certain J/U critical points are influenced across a full range of φ in Fig. 15.From these results, we observe some very immediate similarities to what was seen with the tilted potentials.Moving away from φ = 45 • , the lobe-like structure and extruding features of the BG are again destroyed, and the critical points remain stable across the majority of φ.This is be expected, as the skewed potentials will also no longer stabilise precise lines of weak modulation throughout the lattice.As we approach φ = 0 • , both critical points will start decreasing towards a fixed value, before a sudden jump.This behaviour is again due to the emergence of crystalline properties within this interval, with a flat distribution of on-site energies at φ = 0 • .We finally note that both the skewed and tilted potentials share many similar properties in their critical behaviour and features on phase diagrams.This is due to the fact that both kinds of potentials can stabilise similar kinds of weakly modulated zig-zag patterns on the lattice, which allows for the formation of mixed phases.

Conclusions
In summary, we have shown the presence of intriguing mixed phases in the manybody 2D AA model.The AA potential is markedly distinct from randomly disordered systems, as it may permit the formation of large, correlated domains that possess weak modulation.These domains can then act as barriers to macroscopic superfluidity and insulation, with the long-range order playing a key role in the percolation of MI and SF clusters throughout the lattice.This can dramatically shift critical behaviour, with the most striking feature being the appearance of BG lobes with sharp, wBG extrusions, which we have now confirmed to exist across a wide range of quasicrystalline distributions.At integer µ/U , the insulating behaviour of the BG is destroyed through the percolation of small SF clusters, resulting in a wSF phase.We have also linked the appearance of these features to changes in the underlying structure of mixed phases.In particular, we find that local SF clusters are either the majority or minority of the phase at irrational wavenumbers, due to localisation across precise lines of weak modulation.
By considering more general AA potentials, we have also shown that unique properties of the mixed phases can still be observed, provided that there are smaller domains of weak modulation.If these domains do not exist in the 2D AA potential, we actually find results comparable to those in randomly disordered systems, but with local structures now possessing long-range order.
Furthermore, we have also studied the importance of long-range variations within the on-site potential in stabilising precise regions of mixed phases.At the superlattice limits of the 2D AA model, there is a small number of unique values in the energy distribution, which implies the presence of both short-and long-range order.As a result, the percolation and transition of mixed phases will become more correlated and uniform, leading to greater shifts in critical behaviour.On a global scale, phase boundaries can dramatically change, and DW phases at finite tunnelling strengths can also form.On the other hand, when we consider rational wavenumbers that approach quasicrystalline distributions, the influence of long-range variations in on-site energy becomes clear.The number of unique values in the energy distribution will be comparable to the total number of lattice sites, leading to the formation of much fewer regions of weak modulation that stabilise SF percolation.The number of unique DW states will also be vast, and possess a strong sensitivity to particle number fluctuations.

Figure 1 .
Figure 1.Plots of the 2D AA potential for the three different regimes used in this study, including (a,d) equal wavenumbers without a rotation θ of the potential, (b,e) equal wavenumbers with rotation and (c,f) unequal wavenumbers without rotation.The first row of figures (a-c) denote the continuum limits of Eq. (2), whereas the bottom row (d-f) shows the 2D AA potential on a discrete lattice.

Figure 2 .
Figure 2. The red line shows the number of weakly modulated lines N L as a function of β.The curve has been smoothed to better illustrate the most prominent features.We also include a background plot which shows the number of unique values (up to 3 significant figures) N V in the energy potential as a function of β.Several examples of the AA potential are plotted at the bottom of the figure for visualisation purposes.As expected, N V is minimised around rational β, but quickly inflates as we approach quasicrystalline distributions.N L , on the other hand, can either be peaked or minimised to zero around rational β.However, N L will remain largely finite across a large range of irrational β.

Figure 3 .βFigure 4 .
Figure 3. Plots of the order parameters at µ/U = 1.0 and λ/U = 0.15, for different wavenumbers β and tunnelling J/U .Each column corresponds to a fixed wavenumber and tunnelling J/U of (a,d) J/U = 0.004, (b,e) J/U = 0.008 and (c,f) J/U = 0.004, where we plot the local (a-c) density ρ i and (d-f) correlation ϕ i order parameters.The wavenumber in (b,e) corresponds to a rational wavenumber of β = 5/7, whereas the other wavenumbers are irrational.In the quasicrystalline limits, we see percolating lines of densities around ρ i = 3/2, which indicate an enhanced transition of the local SF clusters.

Figure 7 .
Figure 7. Phase diagrams of the 2D AA model for fixed modulation strengths λ/U .Here, we consider (a) λ/U = 0.1, (b) λ/U = 0.15, (c) λ/U = 0.2, (d) λ/U = 0.35 and (e) λ/U = 0.525.As we increase λ/U , the MI lobes are slowly reduced in extent, leaving behind larger regions of the BG phase.At strong λ/U , the extruding features appearing from the BG lobes are still persistent.Furthermore, at integer µ/U , the onset of the SF phase occurs at a small J/U < 0.01, with the wSF denoting regions of weak percolation.

Figure 8 .
Figure 8. Plots of the (a,b) S M I and (c,d) S SF local measures for the characterisation of mixed phases.Here, we consider the two phase diagrams in Figs.7(b,d), which correspond to modulation strengths of (a,c) λ/U = 0.15 and (b,d) λ/U = 0.35.The black line in each plot is the percolation transition, while the red line indicates when SF clusters account for more than 50% of the overall state.

Figure 9 .
Figure 9. Plots of the critical behaviour at fixed chemical potentials for λ/U = 0.15, showing the J C /U BG-SF critical point.For certain β, we also show the rational fraction above J C /U peaks and on the β grids as vertical lines.The structure of local SF clusters can be divided into two distinct regimes based on whether or not J C /U is maximised or minimised for each µ/U .These then directly correspond to superlattice and quasicrystalline realisations of the AA potential.

Figure 10 .
Figure 10.Phase diagrams of the 2D AA model for rational wavenumbers β and a fixed modulation strength of λ/U = 0.15.The inset figures represent small portions of the AA potential.We also express the wavenumber to 4 decimal places, for comparison with Fig. 2. The exact wavenumbers are given as (a) β = 29 41 , (b) β = 5 7 , (c) β = 3 4 and (d) β = 11 13 .For superlattice realisations of the potential with smaller fluctuations in energy, wSF and wBG phases will appear less frequently on the phase diagrams.Furthermore, it is also possible to stabilise macroscopic DW order at finite J/U .

Figure 11 .
Figure 11.Average density as a function of µ/U around β = 3/4, when J/U = 0 and λ/U = 0.15.All states are insulating for these parameters.We observe that the average density plateaus in the DW domains become unstable as we move further away from the rational limit in β.This is a consequence of the number of unique values in the energy potential increasing.

Figure 13 .
Figure 13.Plots of the critical behaviour at fixed chemical potentials for λ/U = 0.15, showing the J C /U BG-SF critical point.The critical points remain stable for a large range of θ.As θ → 45 • , the potential becomes more crystalline/uniform, which shifts the critical behaviour.

Figure 15 .
Figure 15.Plots of the critical behaviour at fixed chemical potentials for λ/U = 0.15, showing the J C /U BG-SF critical point.The critical points are again stable for a large range of φ.As φ → 0 • , the potential becomes more crystalline/uniform, which shifts the critical behaviour.