The mean-field Bose glass in quasicrystalline systems

We confirm the presence of a mean-field Bose glass (BG) in 2D quasicrystalline Bose–Hubbard models. We focus on two models where the aperiodic component is present in different parts of the problem. First, we consider a 2D generalisation of the Aubry–André (AA) model, where the lattice geometry is that of a square with a quasiperiodic onsite potential. Second, we consider the randomly disordered vertex model, which takes aperiodic tilings with non-crystalline rotational symmetries, and forms lattices from the vertices and lengths of the tiles. For the disordered vertex models, the mean-field BG forms across large ranges of the chemical potential, and we observe no significant differences from the case of a square lattice with uniform random disorder. Small variations in the critical points in the presence of random disorder between quasicrystalline and crystalline lattice geometries can be accounted for by the varying coordination number and the different rotational symmetries present. In the 2D AA model, substantial differences are observed from the usual phase diagrams of crystalline disordered systems. We show that weak modulation lines can be predicted from the underlying potential and may stabilise or suppress the mean-field BG in certain regimes. This results in a lobe-like structure for the mean-field BG in the 2D AA model, which is significantly different from the case of random disorder. Together, the two quasicrystalline models studied in this work show that the mean-field BG phase is present, as expected for 2D quasiperiodic models. However, a quasicrystalline geometry is not sufficient to result in differences from crystalline realisations of the BG, whereas a quasiperiodic form of disorder can result in different physics, as we observe in the 2D AA model.


Introduction
The Bose Glass (BG) is a well-known phase that arises in disordered Bose-Hubbard models [1][2][3][4][5][6], which separates a direct Mott Insulator (MI) to Superfluid (SF) transition [1,7].The BG contains small regions of a SF nature, but it is a macroscopic insulator; with no transport exhibited across the full system.The phase is characterised by a finite compressibility, the lack of a gap, and an infinite SF susceptibility [1].This means its properties are different from that of the SF or MI states present in the disorder free system.In this work, by using common mean-field methods and a percolation analysis, we will confirm the existence of a mean-field BG phase in quasicrystalline systems.
A ubiquitous example of a quasiperiodic system is the 1D Aubry-André (AA) model, where it is known that the underlying self-duality separates regimes where states are either localised or extended [23], which prevents the formation of a single-particle mobility edge in the spectrum [24].The expansion of the AA model itself to a 2D square lattice has also recently been gaining interest [25][26][27] and shares many of the same properties as that of the 1D AA model.Another method to model quasicrystals is by considering a vertex model of an aperiodic tiling [28][29][30][31][32].In a vertex model, edges of tiles are taken to be bonds and vertices of the tiling pattern are taken to be lattice sites.These kinds of quasicrystalline lattices have been used extensively in prior works to study single-particle physics on aperiodic systems, including their transport [33,34], and topological properties [35][36][37].
The concept of disorder in quantum physics represents the absence of some intrinsic symmetry or structure to a particular system.In this context, impurities and defects to a crystalline material can dramatically change observable properties.In the singleparticle picture, it is known that random variations can result in Anderson localisation [38,39].When coupled with interactions, competing energy scales can introduce further novel properties in the form of exotic phase transitions [40][41][42][43], many-body localisation [44][45][46][47], and glass states [1,[48][49][50].Hubbard models describing random parameters were first introduced in the context of disordered fermionic and spin systems [51][52][53], but have also been studied in the bosonic scenario.In particular, the disordered Bose-Hubbard model has been used to describe superfluid helium in porous media [54], granular superconductors [55], Josephson junctions [56], and ultracold atoms in speckle potentials [57,58].
In this work, we will confirm the presence of a mean-field BG in two distinct quasicrystalline models.First, we will consider a quasiperiodic potential supported on a square lattice, in the form of a 2D AA model.This model will show that the meanfield BG can indeed exhibit self-similar properties of the underlying potential.Second, we will consider quasicrystalline lattices generated from vertex models of aperiodic tilings.In these models -which we will label as disordered vertex models -we have a quasicrystalline lattice with a fully disordered potential.From these two models, we can probe the influence of the quasiperiodic nature arising from either the potential or lattice structure independently.We will show that the mean-field BG is only substantially changed when the quasiperiodic component of the model is present in the potential, with a quasicrystalline form of the lattice only shifting transition points in a minor way.

Quasicrystalline Models
For all models, we will consider an inhomogeneous Bose-Hubbard Hamiltonian, with nearest-neighbour tunnelling, given by where N is the total number of lattice sites, U is the onsite interaction strength, i is an onsite energy offset, J is the tunnelling coefficient, i, j denotes the sum over nearestneighbours, and µ is the chemical potential.The operators bi ( b † i ) are the individual bosonic annihilation (creation) operators with ni being the number operator.

Interacting 2D Aubry-André Model
The first model we will look at is an interacting version of the AA model generalised to 2D.We consider a 2D square lattice with unit lattice spacing and assume that the i of Eq. ( 1) are distributed according to the AA quasiperiodic potential , with λ denoting the modulation strength, β is the wavenumber, and x i (y i ) are the position indices.A phase can be incorporated into the distribution of i and averaged over to allow for a direct analogy of the quasiperiodic system with random disorder [57,[59][60][61][62].However, we are interested in Eq. ( 2) for its quasiperiodic properties, and therefore the introduction of a phase is irrelevant.Note, alternative extensions of the AA model to 2D have been considered, including coupled 1D chains [25] which can exhibit mobility edges [26].
For illustrative purposes, we plot a visualisation of the AA potential in Fig. 1 for a small lattice.As we can clearly see in Fig. 1(c), the discretised potential can form a pattern, with weakly modulated lines appearing throughout the lattice.In a recent study within the single-particle picture, it was shown that these weakly modulated lines prevent a mobility edge from separating localised/extended states and can support ballistic transport [27].As we will see later, the weakly modulated lines will also play an important role in the formation and structure of phases in the mean-field scenario.

Disordered Vertex Model
We will also consider the Bose-Hubbard model with random disorder on lattices defined via the vertex model of aperiodic tillings.The lattices we work with are generated from rotationally symmetric quasicrystalline plane tilings [14,63].Note, the strict group theory definition of a lattice is that of a repeating set of points, and in general the models we study in this work, and any model with open boundary conditions, are labelled as graphs.However, the use of the term lattice to describe models with open  boundary conditions and quasiperiodic sets of points is now ubiquitous in physics and, hence, we will label these sets of points as lattices throughout this work.
The link between aperiodic tilings and quasicrystals was further emphasised from Shechtman's original discovery [64][65][66], allowing for quasicrystals to be interpreted as canonical projections of higher dimensional crystalline lattices [67,68].Cut-andproject sets, in particular, provide a useful framework that we will use to generate quasicrystalline lattices.For further details on this approach, we refer the reader to Appendix A. By using the cut-and-project method, we generate the Ammann-Beenker (AB) and Moore-Penrose (MP) tilings from 4D and 5D hypercubic lattices respectively.In Fig. 2, we show the real-space structure of the tilings, which clearly illustrates the self-similarity and aperiodicity of quasicrystals over large length scales.In a vertex model, we take edges of tiles to be bonds and vertices to be lattice sites.To better visualise the lattice geometry, smaller tilings and their respective vertex models are considered in Fig. 3.Both plots also show the 8-fold and 5-fold rotational symmetry of the AB and MP tilings respectively.
In the disordered vertex model, the potential will be random and will contain no self-similarity, unlike the 2D AA model already defined.In particular, we assume that the i in Eq. ( 1) are drawn from a uniform random distribution, with where ∆ is the disorder strength.However, the underlying lattice geometry is aperiodic, with local disorder in the coordination number and a long-range self-similar structure.

Gutzwiller Mean-Field
To study the phases that may arise in these models, we employ a self-consistent meanfield approach.This is done by using the well-known Gutzwiller ansatz [52,[69][70][71] for the many-body wavefunction, which takes the form where z is the maximum number of particles per site, f i n is the amplitude for n atoms in site i and |n i is the corresponding number state.By separating the many-body wavefunction into onsite products, we neglect quantum correlations, which is valid in the limit for small J/U and a large number of nearest-neighbours between sites.In this form, we can explicitly calculate observables from Eq. ( 4), which will serve as order parameters that govern the systems phase.
First, we can define an order parameter that acts as an effective measure for the correlation properties across the lattice as Next, the density properties of the system can be captured as By using these relations, Eq. ( 1) can be written as a sum of on-site contributions that depend on the previously defined order parameters with each on-site Hamiltonian being given by where ϕ i is taken to be real without loss of generality for the considered phases.In order to then determine the ground state of Eq. ( 8), we use the method of self-consistency for the local problems.Strictly speaking, z should be infinite for a bosonic system.For numerical purposes, however, we limit the maximum z to a sufficiently large cutoff to obtain convergence of the mean-field phases.
In the initial step of the self-consistent loop, we define a set of uniformly distributed random order parameters for each site.The local Hamiltonians are then diagonalised such that the global ground state can be identified.From this ground state, new order parameters are evaluated using Eq. ( 5) and Eq. ( 6) for each site.The loop continues with these new order parameters and the process is repeated until the ground state energy converges to a specified accuracy.At convergence, the mean-field ground state and phase of the system is known.
We finally note here that, within the context of mean-field theory, the way in which we will define phases necessarily differs from exact approaches.Rigorously, one should inspect correlation functions and winding numbers to check for the onset of superfluid order.Since we can not reliably define correlators and other such macroscopic order parameters in a site decoupled model, we instead work with local order parameters.As we will discuss later, however, the application of a percolation analysis can be used as an effective measure for the onset of true superfluidity and allows for the definition of a mean-field BG.

Defining the Mean-Field Bose Glass
In the homogeneous Bose-Hubbard model, we can directly characterise two phases -the MI and SF -based on whether ϕ i ∀ i is zero or non-zero respectively.When we have a non-zero onsite potential present, it is possible to have a third intermediate phase to separate the insulating and superfluid domains, which is the BG.The methods we use here to determine different mean-field phases are based on percolation approaches used in the study of disordered phase transitions [72][73][74][75][76][77][78][79].In particular, we use the approach outlined by Ref. [77], and analyse percolations of superfluid clusters on a discrete function S to define the three phases of the inhomogeneous Bose-Hubbard model.We label sites with integer fillings of density ρ i as MI sites with S i = 0, and non-integer as SF sites with S i = 1.We also define an additional function F in order to determine the onset of the mean-field BG phase, i.e the presence of correlations.For this function, we have F i = 1 for sites with finite ϕ i , and F i = 0 otherwise.
First, the standard MI can be characterised by a uniform, zero S and F across all lattice sites, with no transport.As the tunnelling becomes more significant, a BG phase is expected to appear, which is instead defined by some finite, localised transport properties.In other words, clusters of SF sites will begin to form, but do not percolate across the entire system.These clusters will of course have some finite correlations associated to them, captured by the ϕ order parameters.Furthermore, the phase will have some finite mean-field correlation fraction F, which we define as where N ϕ is the total number of sites with finite ϕ i and N sites is the total number of sites.As the tunnelling is further increased, macroscopic phase coherence is expected to appear, giving rise to the SF phase.Unlike the BG phase, the SF has at least one percolating cluster of SF sites appearing in the discrete function S, allowing for a percolation probability P to be defined as where N span is the number of SF sites in a percolating cluster.The phases and above relations for each case are summarised in Table 1.
Table 1.Phases of the inhomogeneous Bose-Hubbard model Strictly speaking, the BG is defined as an insulating phase, which would imply vanishing ϕ i across the lattice.However, due to the presence of onsite fluctuations from the energy potential, there will always be sites with finite ϕ i in the mean-field BG, with each site essentially having its own MI to SF critical point.Indeed, the corresponding average order parameter φ may then fall below some threshold and be sufficiently close to zero, allowing for the phase to physically be viewed as macroscopically insulating.
We also remark that these definitions differ from other mean-field methods [80][81][82], which usually rely solely on the average order parameter φ and average local compressibility κ to distinguish different phases.With these, the MI and SF can be defined by ( φ, κ) ≈ (0, 0) and ( φ, κ) = (0, 0) respectively.The mean-field BG now has φ ≈ 0 and κ = 0.However, as discussed in Ref. [77], the use of averaged order parameters leads to an exaggeration in critical points for finite disorder strengths, and will fail to predict an intermediate BG phase between a MI and SF phase in certain scenarios [7].
By instead employing a percolation analysis, we can find signatures of global phase coherence in the system, with results from crystalline disordered systems in good agreement to those of quantum Monte-Carlo [3,7,83,84] and tensor network [85,86] approaches.This is because the presence of energy modulations induces pronounced fluctuations to onsite particle numbers for finite J/U , giving rise to non-integer densities within SF clusters.The percolation of these clusters is intrinsically linked to the SF density calculated in quantum Monte-Carlo simulations, and is hence a very reliable measure for the onset of true superfluidity, even in the mean-field limit [77].

Characterising the Phase Transitions
Due to the continuous behaviour of the local order parameters we consider here, we must define threshold conditions for the phases so that we can label different regions accordingly.We will use similar thresholds to those in Ref. [77], but we will motivate this choice in this section by investigating the transition points themselves.To check for integer densities, we evaluate the difference of ρ i with the nearest integer for each site.If this difference is less than 5 × 10 −3 , then the density is considered integer, hence S i = 0, otherwise it is non-integer and S i = 1.In a similar manner, to define F , we look for sites with ϕ i > 10 −2 , i.e sites with finite mean-field correlations.We evaluate P by looking at the maximum extent of each SF cluster.If this extent is greater than βK, where K is the maximum size of the lattice and β = 0.99 is a numerical scaling prefactor, then P takes a finite value as defined by Eq. (10).Otherwise, P will be zero if there is no percolation of the clusters.In the case of a square lattice, K is simply the side length of the square.For the quasicrystalline tilings, K is defined by 2R m , where R m is the maximum radial coordinate of the tiling with respect to the origin (central site).
We show an example slice of the order parameters in Fig. 4 for the 2D AA model.First, it can be clearly observed that the mean-field BG to SF transition is relatively sharp.This reflects the enhanced definition of the mean-field BG via the percolation analysis.Therefore, any variations in the thresholds considered for the BG to SF transition point will result in little change to the critical points considered in this work.The MI to BG transition is not as well defined for the mean-field approach, as the percolation analysis does not give this transition point.This can be seen in the slower transition exhibited by F in Fig. 4. Strictly speaking, the MI phase is destroyed as soon as a single site exhibits non-zero transport.However, this is not a practical definition due to the influence of finite size effects and numerical noise.This requires the setting of a numerically non-zero but effectively zero threshold for the ϕ order parameters in the definition of F, as discussed earlier.However, due to the slowly increasing nature of F observed in Fig. 4, any changes in the threshold conditions can result in moderate changes in the MI to mean-field BG transition points.These changes will never be substantial enough to remove the mean-field BG from the phase diagram of the quasicrystalline models considered here, as F ∼ 1 before P > 0. From considering various cuts in constant µ/U , we have determined that a lattice size of ∼ 5000 sites is sufficient to obtain convergence in the mean-field critical points to an average order of 10 −4 for the BG to SF transition and 10 −3 for the MI to BG transition.The convergence of the transition points with system size is shown in Fig. 4 by the inclusion of four different system sizes, and we will consider full mean-field phase diagrams for different lattice sizes in Sec.5.2.

Interacting 2D Aubry-André Model
We will now turn our attention to the mean-field phase diagrams of the 2D AA model.To reiterate, this model has a crystalline lattice geometry (that of a square) and an onsite potential which is quasicrystalline and given by Eq. ( 2).In Fig. 5 we show the mean-field phase diagrams for a 99 × 99 lattice, which is sufficient for good convergence of the transition points as analysed in Sec.5.1.As can be seen from these results, we that the mean-field BG phase clearly separates the MI and SF domains, as expected.By increasing λ/U , the mean-field BG slowly reduces the extent of the MI phase in Fig. 5(b) and 5(c), with the BG forming lobe-like structures, similar to the MI in a homogenous system.While the behaviour we observe here is comparable to what is seen in other inhomogeneous and disordered models [73,74,77], we also find several key differences between both results as well.First, we find that extruding features can emerge from the mean-field BG lobes at different µ/U , as observed in all three phase diagrams of Fig. 5.We plot the discrete function S i , which is non-zero if a site is in a SF state, across one of these features in Fig. 6 with fixed µ/U = 0.39 and λ/U = 0.15.The three example states in Fig. 6 thus reflect moving further along the extrusion of the mean-field BG.From these figures, the reason for this curious feature on the phase diagram becomes immediately clear.The weak modulation lines discussed in Sec.2.1 prevent the percolation of a macroscopic SF.This is shown clearly near the limits of the extruding feature in Fig. 6(b), where the state is a set of large SF domains, but no single cluster percolates across the full system, resulting in a weak BG phase.Furthermore, the BG has a more intricate and self-similar structure near the start of this extrusion as shown in Fig. 6(a), with the state being constructed of many self-similar SF structures.The extrusion of the mean-field BG is truly an artefact of the quasicrystalline potential, as such weak modulation lines would not exist if the potential was fully disordered.
In normal disordered realisations of the BG, the phase will replace MI lobes and be present across the full range of µ/U , thus it lacks the lobe structure that is associated to the MI domain [1].Furthermore, intermediate disorder strengths will decrease the width of the MI lobes across µ/U [73,74,77].Interestingly, with a quasicrystalline potential, the mean-field BG retains a lobe-like structure even for large modulation strength λ/U , as shown in Fig. 5.However, we do observe the decrease in width of the MI phase normally associated with the appearance of the BG state.
To investigate the retention of the lobe-like structure, we again consider discrete functions S i at fixed µ/U .In this case, the region of interest is at integer chemical potentials, where the SF phase is stabilised to allow for the lobe structure to form.We show a variety of S i for small J/U at µ/U = 1 with λ/U = 0.2 in Fig. 7.The reason for the SF stabilisation and BG suppression can again be linked to the weak modulation lines.Unlike the previous feature in Fig. 6 where the BG was extended due to lines of insulators, the BG is now suppressed at integer µ/U due to the SF phase percolating across lines of weak modulation.Similar behaviour is also observed at other λ/U and integer µ/U .An interesting question is then, how does the appearance of weakly modulated lines change global properties of the phase diagrams for different system sizes?To answer this, we look at two specific λ, one relatively moderate and one large, and plot the corresponding mean-field phase diagrams in Fig. 8 for system sizes between 33 × 33 and 99 × 99.We can immediately see that the mean-field critical points have very little variation within this range of system sizes, and the majority of features are invariant against changes to the system size.This is to be expected from the previous analysis in Sec.5.1.However, for the smallest of lattices, i.e. 33 × 33, the features corresponding to the weak modulation lines are less pronounced.This is again not surprising due to the smaller number of sites in general and, hence, fewer regions in the quasicrystalline potential that will exhibit the behaviour of weak modulation lines.We also note that this is reflective of other quasicrystalline properties that can often be dependent on the system being of a sufficient size [36,87].In other words, the partial dependence on the 100 0 10 0 10 -3 system size is a result of the long-range order present.

Weak Modulation Lines
Our results have shown that the mean-field BG does appear in the presence of a quasicrystalline potential.The quasicrystalline nature results in the stabilisation and suppression of the mean-field BG in different scenarios.To better understand the origin of the weak modulation lines, we now turn our attention to how they arise in the potential.First, we rewrite Eq. ( 2) as We know from the form of the weak modulation lines that they can appear in single rows or columns.By then fixing one position index y i = Y , the appearance of a weakly modulated line where i = 0 can be defined as where k is an odd integer.In other words, if then the corresponding set of sites (x i , Y ) will be a weakly modulated line with i = 0. Note, we could have fixed x i in the beginning and obtained the same condition for a set of points with fixed x, meaning the derivation confirms the presence of weakly modulated lines in both dimensions parallel to the boundary of the lattice.We could also fix a relationship between x i and y i which would allow for the same derivation to show that general weak modulation lines are possible, e.g. at some angle to the boundary of the lattice.Similar weak modulation lines could appear in other potentials of a quasicrystalline nature.Taking β = 1/ √ 2, we plot the function d(Y ) in Fig. 9 by considering the difference of each Y of the lattice with the closest point in the set k/4β, with k an odd integer.In practice the condition given by Eq. ( 13) can be relaxed to allow for weak modulation lines to effectively occur whenever which allows some small fluctuation along the lines due to the cosine potential, but these modulations would be weak.Here, we see that there are indeed a few specific Y that meet the condition d(Y ) ≈ 0, which directly correspond to the weakly modulated lines observed in Figs. 6 and 7. Therefore, the weakly modulated lines are predictable from the form of the potential, and any other quasicrystalline potentials could be checked to see if such lines exist prior to any calculations.From Fig. 9, it is also clear that if the system is not of sufficient size, then there will be fewer weak modulation lines to stabilise/destabilise the mean-field BG, as d(Y ) ≈ 0 will be satisfied less often.This predicts the behaviour exhibited for the smallest system shown in Fig. 8, where the weak modulation line features are less pronounced.
If weak modulation lines are present in the potential, then we would expect the same physics to be observed for the mean-field BG as has been shown in Sec.2.1.This includes the stabilisation of the mean-field BG around non-integer µ/U due to insulating lines in Fig. 6, and the suppression of the BG due to SF percolation across the same lines in Fig. 7.

Disordered Vertex Model
We have so far focused on the quasicrystalline nature of the problem being present in the onsite potential of the model.In considering the disordered vertex model, we now turn to a different scenario, where the onsite potential is disordered and the underlying lattice geometry is that of a quasicrystal.The lattice is constructed as discussed in Sec.2.2 and illustrated in Fig. 3, with the Hamiltonian still given by Eq. ( 1), but now with a uniform random disorder in the onsite potential.Due to the presence of random disorder, multiple mean-field critical points are calculated and averaged over for different random initialisations of the onsite potential.We will consider 250 disorder realisations for each critical point, with the error being sufficiently low from this choice.All critical points are plotted with error bars from the analysis of the disorder realisations, however, the errors are of the order of the point size plotted.
For comparison, we will consider an 8-fold Ammann-Beenker (AB) quasicrystal, a 5-fold Moore-Penrose (MP) quasicrystal, and a square (SQ) lattice.All three lattices have ∼ 3000 sites, with the AB having 3041, the MP having 3056, and the square having 3025 sites.Note, the rhombic tilings of the AB and MP lattice result in vertex models with average coordination number in the bulk ≈ 4, meaning that the square lattice is Based on these points, the phase regions are marked accordingly.The well-known Mott lobe structures are observed in the disorder free scenario when ∆/U = 0.0.As we enlarge ∆/U , we find that the BG phase clearly separates SF and MI domains, before enveloping the MI completely in the strong disorder limit.
a good crystalline comparison.The mean-field phase diagrams for all three of these models across a range of random disorder strengths are shown in Fig. 10.We observe a surprisingly good convergence of the mean-field critical points for all three models.In addition, we do not observe the suppression or stabilisation of the mean-field BG due to the quasicrystalline structure.
It is to be expected that the mean-field phase diagrams would converge between periodic and quasicrystalline lattices with no random onsite disorder as shown in Fig. 10(a) [88,89].When random disorder is included, we observe the same features as expected from previous results on disorder in periodic lattices [1].The results also show that there is no direct MI to SF transition in the presence of random disorder in a quasicrystal, as is the case in crystalline lattices [2,3,7,83,85,86].Small fluctuations in the critical points of the AB, MP, and square lattices are a result of the varying coordination number in the quasicrystals and the difference in the rotational symmetry that must be exhibited by the ground state, i.e. 8-fold, 5-fold, and 4-fold.

Conclusions
We have confirmed the presence of a mean-field BG phase in quasicrystalline models using the Gutzwiller mean-field method for a Bose-Hubbard Hamiltonian.The characterisation of the MI to BG transition within the mean-field approach is introduced by studying the correlation fraction of the state, and suitable thresholds investigated to consistently define the critical point.On the other hand, the BG to SF transition was identified through the use of a percolation analysis, which allows for a consistent definition of a BG like phase in the mean-field.
Two different quasicrystalline models were considered.First, the 2D AA model, which has an underlying periodic square lattice geometry but a quasicrystalline onsite potential.Second, the disordered vertex model, which has a quasicrystalline lattice geometry with a random uniform disorder in the onsite potential.Both these models contain a quasiperiodic component, however, we observe that the different areas in which they manifest is of vital importance.In both models we confirm the presence of a mean-field BG phase, which is an intermediary for the macroscopic MI and SF.
The disordered vertex models show effectively no difference between periodic and aperiodic lattices.The critical points for the MI to SF, MI to BG, and BG to SF show only marginal variation between quasicrystalline and crystalline models.The small fluctuations exhibited in the critical points are accounted for by the non-constant coordination number for the quasicrystals and the differing rotational symmetries that the ground state must exhibit.As all rhombic aperiodic tilings have an average coordination number of ≈ 4, we would expect the mean-field results of this model to hold for other rhombic tilings as well.Furthermore, we would also expect that the critical points will scale with the average coordination number for any disordered vertex model.This would be true as long as the distribution of coordination numbers does not become large and uniform, which would result in a different critical point.However, we note here that we do not know of any aperiodic tiling that would result in such a distribution of the coordination number.
For the 2D AA model, we confirm that the mean-field BG phase is exhibited by periodic lattices with quasiperiodic potentials.In this case, the mean field BG phase is stabilised or suppressed in different scenarios due to the presence of weak modulation lines.We showed that the weak modulation lines can be predicted from a simple analysis of the quasiperiodic potential.This could be useful for their prediction in other models outside of the mean-field analysis considered in this work.For non-integer µ/U , the weak modulation lines can result in a stabilisation of the BG phase over a larger range of J/U than expected, giving the protruding features of the mean-field BG shown in Fig. 5  and 8.We showed that this feature is due to the weak modulation lines blocking the formation of a macroscopic SF state.This can result in large regions of SF (but not of the same size as the state) separated by thinner regions of MI.For integer µ/U , we observe the opposite scenario.The mean-field BG is now suppressed due to the percolation of thin SF layers that form on weakly modulated lines resulting in a macroscopic SF.The combination of this suppression and enhancement at different µ/U conspires so that the mean-field BG gains a lobe-like structure, reminiscent of the ubiquitous MI lobes.This is significant, as the BG is usually formed across the majority of the small tunnelling domains; with little structure as a function of the chemical potential, as is exhibited in the disordered vertex model.The weak modulation lines of the 2D AA model are a result of the quasicrystalline structure of the potential, and could be exhibited by other quasiperiodic potentials.This means that the lobe-like structure exhibited by the mean-field BG could be a signature that the BG phase is due to a quasicrystalline potential instead of a uniform random potential.
Comment -During the finalisation of this manuscript, a study considering optical lattice quasicrystals using the Gross-Pitaevskii equation confirmed the presence of the BG phase in the continuous regime [90].In that work, it is shown that due to the inhomogeneous and quasiperiodic nature of the potential itself, a BG phase can indeed be stabilised in potentially experimentally relevant regimes.While the considered models and approaches are different, the results in this manuscript complement those in Ref. [90].
By inspection, simply projecting all W to T will result in a dense set points in R 2 , with no particular aperiodic structure.The idea is then to limit the projected set of points in R using some cutoff in I.In most cases, this is defined through the unit cell of Z, which we will label as Λ.We now take the points spanned by Λ, apply the previous transformation in Eq. (A.1) and project Λ to I. In forming the tiling pattern for T , we now only accept points in T whose dual in I falls within the convex hull spanned by the projected Λ to I. By using the unit cell as the bounding volume in I, we ensure that multiple tiles/edges do not overlap.
In our work, we have considered the AB and MP tilings, which are generated from 4D and 5D hypercubic lattices respectively.

Appendix B. Rotation Matrices
Here, we give the precise forms of R that were used to construct the quasicrystalline tilings used in this paper.

Appendix B.1. Ammann-Beenker Tiling
The AB tiling can be generated as a canonical projection of a 4D hypercubic lattice Z.A general position in Z can be written as where each a i ∈ Z.The rotation of V to W is defined by a 4 × 4 rotation matrix R and can be written as Appendix B.2. Moore-Penrose Tiling The MP tiling can be generated as a canonical projection of a 5D hypercubic lattice.Similar to before, a general position can be expressed as

,Figure 1 .
Figure 1.Visualisation of the 2D Aubry-André model on a finite lattice.Here, we show (a) the underlying square lattice with unit spacing and (b) the continuum limit of Eq. (2) with period β −1 in the same spatial interval.(c) The 2D AA model can be formed from a combination of (a) and (b) on a larger lattice.The dashed (red) box highlights a particular line of weak modulation for the 2D AA potential.

Figure 2 .
Figure 2. Tilings generated from the cut-and-project method, showing (a) the Ammann-Beenker tiling and (b) Moore-Penrose tiling.In both cases, the tilings are composed of two rhombic unit cells (shaded and unshaded rhombi on the tilings) that fill all of space continuously.For practical purposes, a circular cut-off is set in real space to limit the total number of tiles / lattice sites and preserve rotational symmetry with respect to the origin.The corresponding vertex model of both tilings has (a) 3041 and (b) 3056 lattice sites respectively.

Figure 3 .
Figure 3. Smaller tilings generated from the cut-and-project method, showing the (a) Ammann-Beenker tiling and (b) Moore-Penrose tiling, with the vertex model for each on the foreground.Lattice sites correspond to the intersection of tile edges, where the edges are bonds between sites.The reduced lattices contain (a) 192 and (b) 196 sites.

Figure 4 .
Figure 4. Plots of the mean-field correlation fraction F (blue curves) and percolation probability P (red curves) for system sizes of N = 33 × 33, N = 55 × 55, N = 77 × 77 and N = 99 × 99.The darker curves correspond to the larger system sizes.For each row, we fix the modulation strength to (a,b) λ/U = 0.15 and (c,d) λ/U = 0.525, whereas each column corresponds to a fixed chemical potential of (a,c) µ/U = 0.6 and (b,d) µ/U = 0.2.

Figure 5 .
Figure 5. Phase diagrams of the 2D AA model for fixed modulation strengths λ/U .Here, we consider (a) λ/U = 0.15, (b) λ/U = 0.20 and (c) λ/U = 0.35.As we increase λ/U , the MI lobes are slowly reduced in extent, leaving behind larger regions of the BG phase.

Figure 6 .
Figure 6.Plots of the (a) phase diagram in a reduced interval and (b-d) discrete functions S on a line across the extruding feature with µ/U = 0.39 for λ/U = 0.15.The labelled green circles in (a) have their corresponding discrete functions S plotted for (b) J/U = 0.035, (c) J/U = 0.043 and (d) J/U = 0.056.As we increase J/U , lines of weak modulation throughout the system prevent true percolation until a critical tunnelling rate.

Figure 7 .
Figure 7. Plots of the (a) phase diagram in a reduced interval and (b-d) discrete functions S to show the SF persistence at µ/U = 1 for λ/U = 0.2.The labelled green circles in (a) have their corresponding discrete functions S plotted for (b) J/U = 0.01, (c) J/U = 0.005 and (d) J/U = 0.002.At small J/U , lines of weak modulation allow for percolation to occur, which stabilises the SF phase.

Figure 9 .
Figure 9. Plots of the difference function d(Y ) for the first 99 Y columns.Throughout this range of Y , several points will minimise the difference function to less than 10 −2 , which directly correspond to the columns of weak modulation appearing in the energy distribution.

Figure 10 .
Figure 10.Phase diagrams of the disordered Bose-Hubbard model for (a) ∆/U = 0.0, (b) ∆/U = 0.6 and (c) ∆/U = 2.0.The coloured points indicate the numerically obtained critical points, with blue circles, green squares and red diamonds denoting the square (SQ), Annmann (AB) and Penrose (MP) lattices respectively.Based on these points, the phase regions are marked accordingly.The well-known Mott lobe structures are observed in the disorder free scenario when ∆/U = 0.0.As we enlarge ∆/U , we find that the BG phase clearly separates SF and MI domains, before enveloping the MI completely in the strong disorder limit.

R
will now be a 5 × 5 rotation matrix, which can be written as