The Fractal-Lattice Hubbard Model

,


Introduction
The one-band Hubbard model [1] is one of the simplest and most fundamental models for describing the effect of electron correlation in solids.In its fermionic version, it physically resembles electrons hopping in a lattice configuration, with hopping parameter t and on-site interaction parameter U .It has been solved analytically in one-dimension [2] using an extension of the Bethe ansatz technique [3,4].However, despite its deceptively simple form, analytical solutions in higher dimensions are yet to be found.In this case, quantum simulations [5,6,7], machine learning [8], and numerical techniques such as Density Matrix Renormalization Group (DMRG) [9,10,11,12], Linked cluster expansion [13], Variational Monte Carlo [14,15,16], Quantum Monte Carlo (QMC) [17,18], among oth-ers [19] must be implemented.
The relevance of the model and its versatility are demonstrated by its numerical applicability to a wide range of two-dimensional lattice structures, leading to the exploration of numerous phases of matter.For instance, the Hubbard model is able to describe spin and density waves [20] and, in frustrated geometries, spin liquids [21,9,15,11,16,22].In addition, the model is related to unconventional superconductivity [23,12,18].
Despite the progress made in more than 60 years of its discovery, the study of the Hubbard Model in a non-integer dimension still remains an open problem.Fractals offer a way to explore lattice configurations with non-integer dimensions.These geometrical structures are characterized by their intricate and self-replicating patterns, obtained through iterative processes, in which a basic shape is repeatedly transformed or replicated.
Recently, fractal geometries have been successfully realized in experimental settings, offering empirical validation of theoretical concepts.A fractal lattice with electrons has been engineered and characterized [24].Quantum transport in fractal photonic lattices [25] revealed anomalous behaviour, deviating from the expected patterns observed in infinite regular lattices.In addition, higher-order topological insulators have been theoretically [26,27] and experimentally [28] observed in acoustic experiments on a Sierpinski carpet, in which the Hausdorff dimension is d H ∼ 1.89.Finally, topological corner modes were shown to emerge not only in the above mentioned metamaterials, but also in real materials such as thin layers of bismuth deposited on InSb substrates [29].
Here, we investigate the Hubbard model in a fractal lattice.We focus on the Sierpinski triangle, which has a Hausdorff dimension of d H = log(3)/ log(2) ∼ 1.58.This lattice geometry is non periodic, which challenges conventional approaches, such as band theory, and requires the exploration of alternative methods to analyze the system.Moreover, the fractal lattice that we adopt contains sites with different connectivity and is bipartite, uniquely influencing the particle behaviour.We first study the model in the limit of zero interaction, since we expect the geometry of the lattice to influence mainly the kinetic term.Afterwards, while white sites have connectivity 2, except for the sites at the corners of the triangle, which have only one neighbour.
We study the Hubbard model with Hamiltonian Here, i, j are site indices and σ is a spin index, which refers to the up or down projection of the electron's spin along a fixed axis.The operator n i,σ = c † i,σ c i,σ is the number operator and one of the sums runs over nearest-neighbouring sites ⟨ij⟩ belonging to the set Λ containing the sites indices.The first term on the RHS describes hopping of electrons on neighbouring sites, with hopping amplitude t.For simplicity, we set t = 1 as the energy unit.The second term accounts for Coulomb interactions between electrons with opposite spin on the same lattice site, described by the Hubbard parameter U .We set open boundary conditions.
Due to the bipartiteness of the lattice, Fig. 2(b), the model is particle-hole symmetric at half-filling [1].In addition, the Hamiltonian has a global SU (2) invariance [1], which reflects the spin-rotational symmetry.
We are interested in the ground-state properties of the system.However, a naive exact diagonalization of the many-body Hamiltonian would require the diagonalization of a matrix of size 4 M × 4 M , with M the number of sites.This is computationally demanding, and we could not go beyond the first generation.In the following sections, we apply various methods to overcome this difficulty.3 Tight-Binding Approach In order to study the many-body ground-state properties, we follow an approach that, in the framework of second quantization, makes use of the TB limit.Therefore, we start by investigating the Hubbard model in the limit when the interaction parameter U = 0.The resulting Hamiltonian represents electrons hopping freely around the lattice, without any energy cost for double occupation of lattice sites.
Given N particles that populate the system, one builds the possible single-particle orbitals as superpositions of single-site wavefunctions.These orbitals, in turn, are used as basis of the N -particle Fock space, where we are able to construct anti-symmetrized many-body wavefunctions making use of Slater determinants.We follow closely the formalism detailed in Ref. [33].
We also study the Hamiltonian at the single-particle level, since the many-body behaviour is strictly related.This means that we perform a spectral analysis of the single-particle orbitals that we use to construct the Slater determinant.Notice that computing these orbitals simply consists in diagonalizing the Hamiltonian H T B in the single-particle basis, on which it reduces to the size M × M .

Tight-binding ground-state properties
In this section, we present results of TB implementations on the fractal lattice.In Fig. 3, we show the first generation of the lattice and the behaviour of the ground-state many-body energy as a function of electronic filling N = 2N σ , where N σ = N ↑ = N ↓ is the number of electrons with spin σ =↑, ↓.
The symmetry in the energy distribution reflects the fact that at zero interaction, the system is particlehole symmetric around half-filling.This means that the system is equivalent upon exchange of particles with holes and vice versa.This symmetry is also evident in the configurations displaying the average density per site on the lattice Notice that, below halffilling, the sites with higher connectivity, i.e. the three sites with three neighbours at the centre of the triangles, have a higher density.We can understand this behaviour by taking into account the fact that hopping is energetically favourable, so in the ground-state it is preferable to store electrons in sites from which they have more probability to hop around.Along this logic, the sites with connectivity 2 are less populated, followed by the sites at the corners, which are connected to the lattice just through one link.If we increase the number of electrons and place N σ = 3 electrons as in Fig. 4(b), the Pauli exclusion principle starts to play a role, meaning that two electrons with same spin cannot occupy the same lattice site.The electrons then must occupy sites which are kinetically less favourable, if the most favourable ones are already occupied.When we exceed half-filling, i.e. the case where the total number of electrons coincides with the number of sites, the pattern is inverted: the three corners are then more populated, followed by the sites with connectivity 2 and finally the sites with connectivity 3, as shown in Fig. 4(d).It seems energetically more favourable to store electrons in the corners and let the ones in the bulk hop.Configurations above half filling follow this pattern, see for example Fig. 4(c).

Compact Localized States and Scaling Properties
In order to understand more in detail the behaviour of the system, we must determine the eigenvalues E n and eigenvectors χ n of the single-particle Hamilto- nian.In the TB limit, spin-up and spin-down sectors are decoupled, and we can thus treat them independently.

First generation
The first generation of the fractal lattice has 9 sites and the Hamiltonian is therefore a 9 × 9 matrix.The energy spectrum is symmetric with respect to E = 0, see Fig. 5.This is a direct consequence of the bipartiteness (chiral symmetry) of the system.This bipartite character is maintained for higher generations of the fractal, and consequently their energy spectra are symmetric around E = 0 as well.
Let us now consider the eigenvectors χ 4 , χ 5 , and χ 6 belonging to the three-fold degenerate level E = 0. Figure 6 shows their amplitude on the sites of the lattice.An interesting behaviour becomes evident: the amplitudes on sites that belong to the sublattice B are zero.This is a consequence of the fact that the amplitudes of their neighbouring sites sum up to zero and lead to destructive interference [32].These states are examples of so-called CLS, since they are perfectly localized on a specific set of sites in the lattice, such that their amplitude on the rest of the lattice is exactly zero.
Before continuing, let us remark two interesting points about such states.Firstly, CLS are completely robust against any perturbations-no matter how strong-that only affect the sites on which they vanish [34].This makes these states a candidate for the storage of information in the form of qubits [35,34].Secondly, in periodic systems, CLS usually lead to the emergence of one or more perfectly flat bands, which recently became a topic of intense research interest; see Refs.[32] and [36] for two reviews on flat-band systems.The existence of the three zero-energy CLS, χ 4 , χ 5 and χ 6 , can also be deduced as follows: For a bipartite system with N A (N B ) sites in the A (B) sublattice, it is a classical result that there must be at least ∆N = |N A − N B | zero-energy eigenstates which have vanishing amplitudes on the minority sublattice [37,38].In our case, we have N A = 6, N B = 3, yielding ∆N = 3, which is exactly the number of zeroenergy CLS that we found.
We are now able to better understand the behaviour of zero-energy states: the three CLS χ 4 , χ 5 and χ 6 discussed above form a basis in the degenerate subspace, meaning that a general eigenstate with eigenvalue zero is found by taking linear combinations of χ 4 , χ 5 , and χ 6 .Any pair of these three states share one corner.For instance, eigenvectors χ 4 and χ 5 share the top corner; a linear superposition χ 4 + χ 5 will thus have enhanced amplitude at this shared corner.Equipped with these insights, we can therefore understand why states around half-filling contribute to the density in the corners.As an example to make this point clearer, we show the density configuration χ 4 + χ 5 + χ 6 in Fig. 7.
Importantly, the existence of CLS and the subsequent density-enhancement at the corners of the Sierpinsky triangle is not limited to a simple TB picture.Indeed, as we show in Section 3.3, it persists even when considering intrinsic spin-orbit coupling.Moreover, the scaling of the number of CLS-as investigated further below-shows the same behaviour in both cases, that is, with and without intrinsic spinorbit coupling.
Having dealt with the CLS, let us now turn our attention to the other eigenstates, in particular, to the two-fold degenerate ones.Their degeneracy is a consequence of the C 3v point-symmetry group [39], which has order 6 and contains the identity, two rotations C 3 , C 2  3 and three reflections.This group has an irreducible representation of dimension 2, which explains two-fold degeneracies [40].As a side note, these levels do not posses the property of CLS states.There is not a definite group of lattice sites where the amplitude is zero for every eigenvector of the degenerate level, as it was the case for the zero-energy states.
It is now possible to look at the density configurations in the first generation of the fractal lattice, Fig. 4, from a different perspective.Configurations below the degenerate level around half-filling populate the sites belonging to sub-lattice B with an average density ⟨n i ⟩ = 1, see Fig. 4 (b).This is in agreement with the fact that, due to particle-hole symmetry, the average density per site at half-filling is ⟨n i ⟩ = 1 and the zero-energy CLS only populate sublattice A. Starting from the configuration with N σ = 3 and raising the number of electrons until N σ = 6 means including zero-energy CLS in the many-body wavefunction.The total many-body energy does not increase because the CLS added have zero energy, but sites belonging to sublattice A, in particular the corners, become more populated.
In the following, we analyze the emergence of CLS in higher generations of the fractal.We shall see that there appear more CLS, both at zero and finite energy.

Second generation
The energy spectrum of the Hubbard Hamiltonian on the second-generation of the fractal lattice is shown in Fig. 8. Once again, due the fractal's bipartiteness, the energy spectrum is symmetric with respect to E = 0. Now, there are N A = 15 and N B = 9 sites in the two  respective sublattices, which results in eigenstates with zero energy.Interestingly, these six zero-energy CLS can be generated from just two prototypical states, show in Fig. 9.To obtain six states from these two, each of them must be rotated by 120 • and 240 • .The resulting set of six states are linearly independent, and thus span the entire six-dimensional degenerate subspace of eigenvalue zero.
Apart from these six zero-energy states, the second generation now features a new type of CLS at the non-degenerate energy levels ±1.These states are depicted in Fig. 10.We note that these two states are tightly related to each other: One can be constructed from the other, simply by flipping the sign on all sites of the sublattice B. This is again a direct consequence of the fractal's bipartiteness.
As we shall see, in the third generation, CLS appear at even more energies than just 0 and ±1.To ease the discussion, we will enumerate the CLS according to the absolute value of their energies.We denote the zero-energy CLS as type-1, and the ones at E = ±1 as type-2.

Third generation
The energy spectrum for the third generation TB model is shown in Fig. 11.Apart from the type-1 and type-2 CLS that we encountered in the second generation, there are now five additional types of CLS, namely, type-3 to type-7.Examples of some CLS are  depicted in Fig. 12. Again, the fractal's bipartiteness means that each of the depicted states has a partner whose amplitudes on each site of the B-sublattice is flipped, and which has the same energy, but with an inverted sign.
Let us now discuss some aspects of the CLS in more detail.The (non-depicted) type-1 CLS are simple: they vanish on the entire sublattice B. We note that there are N A = 42 and N B = 27 sites in the two sublattices, resulting in 42 − 27 = 15 type-1 (zeroenergy) CLS.The CLS of type-2, 5, and 6 [Fig.12(a),  (c) and (d) respectively] have a very peculiar property: they vanish on all three other corners of the fractal.On the other hand, the CLS of type-3 [depicted in Fig. 12(b)] and type-7 (not depicted) do not have this property.CLS type-7 exhibit destructive interference on the same set of sites as CLS type-3 and were not included in the graphical representation of Fig. 12 to avoid redundancy.The same reasoning applies to CLS type-4 and we only show CLS type-6 in the figure.

Higher generations
The number of CLS grows considerably at larger generations.For the type-1 CLS-which is caused by the fractal's bipartiteness-, for instance, we can see that for the n-th generation, the number of sites in the two sublattices is given by N zero-energy CLS, as predicted by Lieb's theorem.For instance, in the 8th generation, there are N = 16404 sites in total, with an imbalance in the two sublattices given by N A = 9843 and N B = 6561.The number of zero-energy CLS is thus ∆N = 3282.The plot of density of states (DOS), in Fig. 13, provides a clear visualization of the zero modes present in our fractal lattice.An analysis of the DOS in the tight-binding regime on the Sierpinski lattice was done in Ref. [24] to show the self-similarity at different scales.For large n, the first term in Eq. ( 2) becomes dominant, from which we conclude that the ratio α = ∆N (n + 1)/∆N (n) thus equals 3. Since the sides of the fractal are duplicated whenever going from one generation to the next, the dimension is defined as and we thus observe that the dimensionality we assign to number of zero-energy CLS tends to the dimension of the fractal, d H = log 3/ log 2.
By explicitly diagonalizing the system for the first eight generations, we also studied how the number of the type-2 to type-7 CLS changes when considering higher generations of the lattice.The results are presented in Fig. 14.The number of CLS of types 1, 2, 4, 5 and 6 tends to scale as the Hausdorff dimension when the generation number increases.
We notice that CLS of type-4, 5 and 6 share the same scaling, even if they belong to different energy levels, see Fig. 11.Regarding CLS of type-3 and type-7, we found that the number of states in the corresponding levels ± √ 2 and ± √ 5 does not increase but remains the same for odd generations, and is zero for even generations.Motivated by these scaling behaviours, we also studied how the many-body ground-state energy E T around half-filling (i.e. the value of energy corresponding to the degenerate level in Fig. 3(b) changes when increasing the generation.As shown in Fig. 15, these values also scale as the Hausdorff dimension.Therefore, they can be derived by recursive computations upon increasing the generation, without the need of diagonalizing the Hamiltonian.This also holds for the number of CLS in some energy levels, and constitutes a great advantage because the com-g2/g1 g3/g2 g4/g3 g5/g4 g6/g5 g7/g6 1.50 putational cost increases rapidly as we increase the generation.
As a last addition to the analysis of the zero-energy level, we investigated the average density on the corners in the last configuration of the zero-energy level of the many-body spectrum.This means, for example, filling N σ = 6 for the first generation, N σ = 15 for the second, and so on.We noticed that its value is ⟨n⟩ = 1.5909 for every generation that we implemented, very close to the Hausdorff dimension.

Compact localized states when considering intrinsic spin-orbit coupling
So far, we have seen that CLS naturally appear in the Sierpinsky triangle when treated within a simple tight-binding formalism.To emphasize the importance and universality of these states, we will demonstrate in the following that they are still present in the model even when taking intrinsic spin-orbit coupling (ISOC) into account.
To introduce ISOC into the model, we start from the TB Hamiltonian of the Sierpinsky triangle and add the term to the Hamiltonian.Here, the sum goes over pairs of sites that fulfil both of the following conditions: (i) they both are in the majority sublattice, and (ii) in the original Hamiltonian without ISOC, they are two sites apart from each other, i.e. the next-nearest neighbor (NNN).The coefficient v i,j = +1 when the two-hop path from site i to site j is in the clockwise direction, and v i,j = −1 when it is in anti-clockwise direction.The Hamiltonian includes intrinsic spin-orbit coupling (SOC), Rashba SOC, and a staggered mass in a honeycomb tight-binding lattice [41].Since our focus is on generating topological features, we considered only the intrinsic SOC term in the Hamiltonian because this is the term that leads to the quantum spin Hall effect.The staggered mass opens a trivial gap and the Rashba SOC tends to close the topological gap.It is important to note that, because spin is a conserved quantum number, it can be omitted from the Hamiltonian during calculations.At the conclusion of the calculations, the spin can be reintroduced, with the understanding that all results are valid for the opposite spin as well, but with an inverted sign for edge propagation.Effectively, Eq. ( 3) adds complex-valued coupling (with amplitude b) between NNN sites-similar to Ref. [42]-, but only if the coupling does not cross the "empty regions" of the Sierpinski triangle.In Fig. 16(a), the Hamiltonian is shown in pictorial form, with complex couplings denoted by dashed, arrowed lines.Let us now analyse the situation in detail.In contrast to the case without ISOC, the first few generations (we have checked until the fifth) only show one type of CLS, which lies at zero energy.Moreover, in contrast to the case without ISOC, the CLS is nondegenerate in the first generation; the state is depicted in Fig. 16(a).Remarkably, it is exactly the state from Fig. 7, so that the amplitudes are enhanced at the three outer corners of the fractal.The ISOC entangles the three degenerate eigenstates into a sum of the them.
In the second generation, the CLS become two-fold degenerate.The first of the corresponding eigenstates is shown in Fig. 16(b).We note that it is similar to the one from Fig. 16(a), and indeed can be obtained from the latter by simply copying it three times.In the following, we will call this eigenstate the ubiquitious state.It sill vanishes on a large number of sites, which are marked in Fig. 16(b) in black.Again, as with any CLS, the cause is destructive interference; see inset of Fig. 16(b).The second eigenstate with zero energy is what we call a closed-loop state.It loops around the central gap (white space) of the fractal and vanishes outside of it.We further note that the amplitudes of this closed-loop state depend on the coupling t and the strength b of the intrinsic spin-orbit coupling.
Similar to Fig. 7(d), we can superpose the two degenerate zero-energy states to obtain a state which has a larger enhancement on some sites.This is demonstrated in Fig. 16 (d) and (e).In higher generations, the overall picture does not change much.For instance, in the third generation, there are five zero energy eigenstates: One ubiquitous (Fig. 17(a)), one large closed-loop state for the central gap (Fig. 17(b)) and three smaller closed-loop states for the smaller gaps (Fig. 17 (c), (d), and (e)).Once again, one may superpose these states to obtain states with enhanced amplitudes on certain sites (corners, for instance).
We have checked the number N of zero-energy CLS until the fifth generation.For each generation, this number is equal to one third of the number of zeroenergy states without ISOC.Thus, the scaling of the number of such states is the same, with or without ISOC.

The Hubbard Model
We now introduce interaction to the model in the fractal lattice and present both the numerical methods used and the results of our implementations.
The full Hamiltonian in Eq. ( 1) contains a twobody term, which prevents us from diagonalizing the system at a single-particle level, as done for the TB model.In the next subsections, we present three distinct methods that we use to investigate this Hamiltonian numerically.

Exact Diagonalization
The first method consists in naively diagonalizing the Hamiltonian in a many-body basis of the Fock space and identifying the ground-state energy as the lowest energy eigenvalue and the ground-state wavefunction as its corresponding eigenvector.To diagonalize the Hamiltonian, we first need to choose a basis and our convention for the basis choice is, in the Fock space, strings with 2 × M entries that can take values of either 0 or 1 and that identify the positions of the electrons on the lattice, 1 means that there is an electron, 0 that there is none [43].We also need to define an order: the first half of the entries represents spin up particles populating the numbered lattice sites, while the second represents spin down particles with the same lattice site ordering.The number of such basis states is 4 M .Since this number grows exponentially with the number of sites, exact diagonalization is a computationally demanding method to solve the Hamiltonian and it will only be used to solve small-size systems with a small number of electrons.

Mean-Field Approximation
One way of circumventing the two-body issue is to perform a mean-field Hartree-Fock approximation.This consists in rewriting the interaction term as [44] ) which allows us to express the Hubbard Hamiltonian in terms of one-body operators and proceed by constructing the many-body ground state in the same way as we did in the TB limit.The difference is that the mean-field term is not known, since the averages present in Eq. ( 4) require the knowledge of the solution, which is, in turn, what we want to compute.The way we proceeded is by applying a self-consistent iterative scheme.Starting from a randomly generated Hamiltonian, we use the diagonalization approach introduced in the previous section to compute the many-body ground-state energy, average density per site, and wavefunction.The latter is then used to build a mean-field Hamiltonian H M F , from which we again compute the above mentioned quantities.This scheme is iterated until the solutions converge.We show later that this approximation is valid for very small values of the interaction parameter U .

CP-AFQMC
The method that we largely implemented to study the Hubbard Hamiltonian is CP-AFQMC [33,45,46,47].The starting point is a mean-field ansatz, found by following the procedure outlined in Section 4.2.The imaginary-time evolution of this ansatz tends to the ground-state of the Hubbard Hamiltonian as the imaginary time becomes large.Implementing imaginary-time evolution requires a Hubbard-Stratonovich transformation, which introduces external auxiliary fields at each lattice site.This enables to go from an interacting system to a non-interacting system living in a space of fluctuating external fields.The wavefunction is written as a linear combination of many wavefunctions, the walkers.Evolving the full system in imaginary time means evolving each of the walkers.
In addition, more techniques are required to render the method efficient.Importance sampling is used, and the simulation is guided by the initial ansatz throughout the evolution in imaginary time.Back propagation is required to compute observables that do not commute with the Hamiltonian.A constrained path approximation makes it possible to deal with the sign problem that stems from the fermionic nature of the particles populating the system.For a more detailed description of the method, we refer to Appendix A.
Our starting point was the validation of the QMC method by comparing its results with exact ones, i.e. with results found by performing exact diagonalization (ED).We also compared these two methods with the mean-field (MF) approach to understand to which extent this approximation is appropriate.Since the simulation time required by ED scales exponentially with the system size and the electron number, we considered a small system (the first generation of the fractal lattice), with a couple of electrons (N σ = 1).The quantities computed were the total, kinetic, and potential energy of the many-body ground-state, for values of U ∈ [0.1,9].The potential energy is the contribution to the energy given by the Coulomb interaction in the Hubbard Hamiltonian in Eq. ( 1).In Fig. 18, we present the behaviour of the energies computed using these three different methods, as the interaction strength increases.Strictly speaking, the mean-field (MF) approximation is valid only for U = 0, but it is a reasonable approximation for very small values of the interaction U when U → 0. On the other hand, the QMC approach follows very accurately the exact behaviour.For this reason, we focus on the QMC implementation and use MF only to make it more efficient.We also notice that the value of the kinetic energy computed by QMC for an interaction U = 6 deviates from the exact solution and coincides with the MF solution.This could be a consequence of the fact that the MF ansatz, used in an importance-sampling scheme and in the CP approximation, is biasing this result.
To understand the extent of this influence, we per- formed simulations in the first generation with N σ = 6.We studied the behaviour of the QMC energies as the interaction strength is increased, with the simulation being guided by different MF ansätze.In particular, we performed simulations where the MF ansatz is generated using U eff = U , where U eff is the interaction strength used in the MF simulation, referred to as effective potential.This behaviour is indicated by the inverted triangles in Fig. 19.Then, we performed simulations with U eff = 0.1.The simulation that provides the smallest total energy is the best approximation of the ground state, due to the variational principle.It is possible to observe that the latter ansatz gives better results for every value of interaction.From this analysis, we can conclude that a preliminary study on the influence of the MF ansatz on the QMC simulation has to be performed before proceeding with the implementation.In fact, since the path of the walkers is influenced by the MF trial wavefunction, a bad ansatz affects negatively the paths, making them deviate from the exact solution.Moreover, to handle the sign problem, we introduce a constrained-path approximation that renders systematic errors.However, those errors are small since the results coincides very closely with exact diagonalization, as shown in Fig. 18.

Interaction and density distribution
We now consider the configuration where N σ = 6, which corresponds to the last energy state in the degenerate energy-level found when U = 0, see Fig. 3

(b).
In the TB studies, we understood that this state is a consequence of the CLS, with destructive interference happening at the sites in the centre of the triangles.We also observed that the density at the corners has a value close to the Hausdorff dimension.Studying this configuration allows us to tackle the consequences of interaction in this type of states.We can divide the lattice into three groups of sites: sites 1, 6, 9 called corners, sites 2, 5, 7 called center and sites 3, 4, 8 called connections, see Fig. 3(a) for the site indexing.We expect the density to be the same in each of these groups, since there is no reason for an imbalance.Figure 20 shows how the average density per site ⟨n i ⟩ changes in these groups of sites when increasing the interaction strength.We observe that for small interactions, both corners and connection sites are more populated than the center sites.As the interaction increases, their density decreases, while it grows in the sites at the center.The electrons start spreading towards the center of the lattice, and they keep spreading until the density is approximately homogeneous on every site, at a value of approximately ⟨n i ⟩ ∼ 1.35 and U ∼ 8. Thereafter, the density in the center sites continues to increase, while decreasing slowly on the rest of the lattice sites, meaning that electrons start to accumulate in the center sites.For very strong interaction values, it becomes more favourable to have electrons in sites with more connectivity, where hopping is more probable.
Relating this study to the TB considerations and the CLS type-1 discussed in Section 3.2, we expect that those types of states with destructive interference on center sites get destroyed as soon as the interaction is turned on, since the density on those sites increases.

Quantum phases in the second generation of the fractal lattice
The quantum phases of the Hubbard Model have been intensively studied on various lattices, in particular for the case of half-filling.We consider here the second generation of the fractal lattice, Fig. 2, which has 24 sites, at half-filling, N σ = 12.
Let us start by studying the behaviour of the mag- netisation.We define the local magnetization m i and the total magnetization m α per sublattice as where the index α refers to one particular sublattice, α = A, B, and Λ α refers to the set of site indices belonging to the α sublattice.Fig. 21(a) shows the magnitude of the local magnetization for weak interaction on the lattice structure.We see that the magnetization on different sublattices has opposite sign, and the magnitude in one is larger than the magnitude in the other, thus characterizing a ferrimagnetic phase.A mean-field study of the Hubbard model in fractal-honeycomb lattices also finds spontaneous spin polarization [48], which further corroborates our observation of a ferrimagnetic state.The local magnetization is related to the projection of the spin along a certain axis.Since the sum of local magnetizations on the lattice is clearly not zero, we find that the system equilibrates to a spin imbalanced configuration.In particular, we find that N A = 15 and N B = 9.Without loss of generality, we investigate the case with N ↑ = N A and N ↓ = N B .The opposite case is obtained by inverting the positive direction along the axis where the spin is projected.Figs.21(b) and (c) show the distribution of the 15 spin-up and the 9 spin-down electrons on the lattice, in the weak interacting regime.We can interpret Fig. 21(b) by considering that the 9 spin-up and spin-down electrons fill the single-particle energy levels with negative energy, see Fig. 8.The remaining 6 spin-up electrons are placed in the degenerate zero-energy level.In this level, the CLS have destructive interference on the sites in the centre (B sublattice); sites on sublattice A get more densely populated.Therefore, the origin of imbalanced magnetic properties at half-filling and weak interaction can be reconnected to the zero-energy CLS that we found at a TB level, which also vanishes at the B sublattice.Now, we want to investigate the behaviour of the system when increasing the interaction strength.First, we notice that the imbalance between spin-up and spin-down electrons remains when increasing the interaction strength.Regarding the local magnetization, in Fig. 22 we show the total magnetization in the two sublattices as a function of interaction strength.We observe that the sum of the magnetization in the two sublattices almost does not change.This quantity represents the total magnetization along a projection axis, and since the ratio of number of spin-up and spin-down electrons remains unchanged, we expect a constant behaviour.However, we need to take into account the fact that the number of spin-up and spindown electrons is an output of the QMC simulation and, thus, subject to subtle fluctuations.For a regular periodic lattice, the theory predicts antiferromagnetic behaviour in the strong coupling regime, where the average local magnetization per site in both sublattices is m α .The reason why we are not able to see this behaviour is due to the imbalance in the number of electrons that populated the system which, in turn, is a consequence of the geometry of the lattice.
Another quantity that can be studied to determine the phases of the system is the doublon density, defined as D = 1 where the sum runs over the lattice sites.While being easy to compute, the double occupancy has the advantage of being related to the metallic behaviour of the system.Moreover, studies on periodic lattices suggest that the behaviour of D as a function of interaction should differentiate metallic and insulating phases [18,49].In particular, Brinkman and Rice obtained, through a variational calculation, that the doublon density behaves linearly in the metallic phase [50].Even though this result was derived in the context of conventional band theory, we also observe in critical interaction is the same as the one describing a Mott transition from a paramagnetic conducting state at small values of U , to an antiferromagnetic (AFM) insulator at U > 4.5 in the 2D honeycomb lattice [51].
After this value of interaction strength, the behaviour deviates from the linear one.
To further investigate the magnetic nature of the system, we consider the magnetic correlation function which quantifies the correlation between the local magnetization of a pair of electrons, one placed at site i and one at site j.
From the magnetic correlation, one can compute the magnetic structure factor, defined as where r ij = r i − r j .The symbol k indicates the momenta in the reciprocal space.It is possible to define a similar quantity, where instead of considering magnetic correlations, one considers density correlations.The computation of magnetic and density structure factors is a common method used to investigate both the metallic and spin phases of the system.The magnetic structure factor is used to detect magnetic order since it shows peaks at the K points of the Brillouin zone [18].This method can be used here because it does not require knowledge of allowed momenta in the reciprocal space, which is not well defined for fractals.In Fig. 24 (a), we show the magnetic structure factor in the strong interaction regime, where one observes the formation of peaks, which are a signature of magnetic order.We then computed the magnetic structure factor for various values of interaction and found a similar structure for each value.This seems to suggest that the system has spin order for any value of interaction strength.Before continuing, let us mention that the six-fold rotational symmetry visible in Fig. 24 (a) can be explained using the symmetries of the setup; see Appendix B.
To validate the presence of magnetic order, we look at the radial behaviour of the magnetic correlation defined in Eq. ( 5).In particular, we define the angleaveraged radial magnetic correlation, where n r is the number of pairs {i, j} that have the same distance r.This quantity averages the magntic correlation functions of pairs of electrons located at sites r far from each other.The set of possible distances r is discrete, and it is important to notice that, for a given distance, all the pairs connect sites belonging either to the same sublattice or to different sublattices.In Fig. 24(a), we show the dependence of this correlation on the radial distance r.The behaviour shows positive and negative correlations, which indeed can indicate antiferromagnetic or ferrimagnetic order.One of the properties of anti and ferrimagnetism is spin orientations that alternate in the lattice.This means that electrons on sites belonging to the same sublattice should have spins oriented to the same direction (the correlation is positive) while electrons on sites belonging to different sublattices should have spins pointing in opposite direction (the correlation is negative).To verify this behaviour, we multiply the spin radial correlation by a factor (−1) α , where α is even (odd) if sites i, j belong to the same (different) sublattice.The quantity obtained will be referred to as radial staggered magnetization, The result is shown in Fig. 24(b).We can conclude that the property of opposite spins in different sublattices is satisfied, as expected from the analysis of the local magnetization for the different lattice sites.
In order to study the magnetic order at different values of interaction, we investigate the total staggered magnetization, where the sum is over the number i r of possible distances N r .This quantity represents the average staggered magnetization over different values of radial distances.Its value for U = 6 is plotted in Fig. 24(b) as a reference to the eye.Its behaviour as a function of interaction is shown in Fig. 25, where we see that the value increases with the interaction before the critical value around U c ∼ 4.5.After it seems to oscillate without an overall growth or reduction for larger U .The fact that it always assumes a vanishing value for every value of interaction proves the existence and resilience of magnetic order.In Fig. 25, we also show the staggered magnetization obtained with QMC calculations in a recent work on a honeycomb lattice [52], the regular version of our fractal lattice.
In the regular case, the phase before the transition is not magnetic, as the total staggered magnetization is zero.Finally, the yellow triangles in Fig. 25 show the total staggered magnetization computed using the MF approach under the same conditions.The comparison with MF results at half-filling underlines the necessity of powerful methods such as QMC.In fact, QMC simulations equilibrate to a steady state with an imbalanced number of spin-up and spin-down electrons, even when initializing the calculation with balanced trial wavefunctions.This results in a finite magnetic order, also at weak interactions.In contrast, MF calculations at half-filling preserve the number of spin-up and spin-down electrons balanced, giving a non-magnetic state at U = 0.As expected, the MF result shows a spurious transition to a magnetic phase.
Overall, we observe the formation of a ferrimagnetic order, which can be related, at weak interaction, to the zero-energy CLS states unveiled at a TB level and, in general, to the imbalance in the number of spin-up and spin-down electrons.The latter is a consequence of the geometry of the lattice.Therefore, we observe a phase transition from metallic to insulating phase, suggested by the change in behaviour of the doublon density and the maximum of the magnetic correlation function in k space.Since this phase transition is driven by interaction between electrons, it is a Mott transition.

Conclusion
In this work, we solved the Hubbard Model for fermions on a fractal lattice, using various methods.
Initially, we examined the non-interacting limit of the model, employing a TB approach.It success- For the second generation of the Sierpinski fractal, the calculations were performed using MF (yellow triangles) and QMC (orange crosses).For the regular honeycomb lattice, we extracted the QMC data from Ref.
fully confirmed the presence of particle-hole symmetry through the examination of the average density distribution on the lattice and the symmetry of the total many-body energy.Additionally, by analyzing the energy spectrum of the Hamiltonian at the singleelectron level, we identified the emergence of CLS at various energy levels.In the first generation, we observed the formation of a specific type of CLS, characterized by zero energy and destructive interference on sites with a connectivity of 3. Moving to the second generation, we witnessed the appearance of an additional type of CLS, manifesting at energy levels of t and −t, and displaying destructive interference along the reflection axes.Progressing to the third generation, we observed the emergence of more diverse types of CLS, forming at new energy values and exhibiting destructive interference on other sublattices.Remarkably, the number of CLS at each energy level increased in higher generations of the fractal, scaling with the Hausdorff dimension of the fractal d H = log 3/ log 2 ≃ 1.58.Moreover, we identified CLS of type-3 and 7, emerging on energy levels ± √ 2 and ± √ 5, respectively.The number of these two types of CLS remains the same for odd generations and is zero for even generations.Finally, we discovered that the density at the corners, for configurations with a number of electrons that fills the degenerate level at zero energy, closely approximates the dimension of the fractal across all generations of the lattice.This result resembles the findings of higher-order topological insulators realized using acoustic quantum simulators.Indeed, the outer corner modes were found to exhibit the same dimension of the Sierpinksi carpet [28].Unlike in their framework, this interesting connection between outer modes and fractal dimension arises in the absence of any external magnetic flux.We note that CLS in Sierpinski fractals have been previously con-sidered, both theoretically [53,54,55] and experimentally [56,57].In addition, superconductivity [58] and plasmonic properties [59] were investigated in fractal structures.However, these works use a different arrangement and number of sites, since they do not consider the central site in each small filled triangle [cf.Fig. 2(a)].Their results are thus not directly applicable to the setup considered here.While constructing fractal structures, one may consider sites at the vertices of the fractal, sites at the center of the "bulk" parts of the fractal, or both.The first is the usual "fractal lattice" construction and the second leads to the dual structure, which has the same Hausdorff dimension, but a different distribution of voids.A combination of both cases is considered here.Although a direct comparison with Refs [53,54,55,56,57] is not possible, previous studies on fractal lattices showcase the relevance of fractal geometries, demonstrating, for example, their ability to host anyonic excitations [60,61].We investigated the model in the presence of intrinsic spin-orbit coupling.We found that the zero-energy eigenstates become entangled and lead to high-intensity corner modes in the inner and outer triangles.Subsequently, we numerically implemented the full model, including the interaction term, using three distinct numerical simulations.The first approach involved exact diagonalization, which is an exact method but limited in computational scalability.Therefore, we employed it primarily for validating the results obtained from the other two numerical methods, specifically for smaller system sizes and lower electronic fillings.The second approach entailed employing a mean-field Hartree-Fock approximation, which offers analytical insights but is observed to be valid only for relatively weak interaction strengths.The primary method, extensively employed in this study, is CP-AFQMC.
When examining the effects of interaction, the application of CP-AFQMC enabled a comprehensive exploration the density configurations.By concentrating on the electronic filling value that displayed density patterns closely aligned with the dimension of the fractal at the corners, we observed notable changes in the density distribution upon introducing interaction.We demonstrated that under strong interaction, sites with higher connectivity, where destructive interference occurred under weak interaction, became more populated as a countermeasure against the effects of interaction.We observe a decrease in the density at the corners as soon as interaction is introduced, leading to the conclusion that the CLS of type-1 are not robust in the presence of interaction.
We further employed CP-AFQMC to investigate the quantum phases of the system within the second generation of the lattice at half filling, considering non-zero interaction strengths.Initially, we observed that the simulation reached an equilibrium state char-acterized by an unequal number of spin-up and spindown electrons, corresponding to the sizes of the two sublattices.Our subsequent focus shifted towards studying the magnetic order within the system.To accomplish this, we computed the local magnetization per site and discovered that neighbouring sites exhibited opposite signs and varying magnitudes.This was consistently observed across all interaction strengths and indicated the presence of a ferrimagnetic phase, see Fig. 21(a).We were able to establish a connection between the imbalance in local magnetization and the type-1 CLS under weak interaction.Subsequently, we verified the magnetic order by observing prominent peaks in the magnetic structure factor and long-range order in the magnetic correlations.Lastly, we detected a Mott transition by analyzing the behaviour of the tail of the magnetic correlations and the doublon density.This transition occurs at an approximate value of U C ≃ 4.5, which intriguingly corresponds to the critical value detected through QMC methods for a Mott transition on a honeycomb lattice [51], the periodic lattice configuration most similar to our chosen fractal structure.The difference is in the magnetic order, which we find to be ferrimagnetic for every value of interaction strength.The latter result, in turn, agrees with studies of Hubbard model on non-periodic lattices, such as the two-dimensional hexagonal golden-mean tiling [62].
Concerning the effect of spatial dimensionality in the phase diagram of the ground state of the Hubbard model at half-filling, important differences are observed.In one dimension (1D), no Mott transition at finite interaction is predicted, with the system being metallic only for U = 0 and an AFM insulator for U > 0 [63].Moreover, the low-energy excitations of the 1D Hubbard model bear resemblance to Tomonaga-Luttinger liquid theory [64,65,66].In 2D, the aforementioned phase diagram changes depending on the geometry considered.For the honeycomb lattice, which is the most similar to the fractal lattice considered in this work, a Mott transition is predicted at finite interaction [52].Below the critical interaction U C , the system describes a non-magnetic semimetal, whereas above U C , the honeycomb-Hubbard model describes an AFM insulator.Our calculations consider a spatial dimension that lies between 1 and 2, and new physics was observed.Similar to the 2D honeycomb model, a Mott transition at U C ≃ 4.5 was observed, but no AFM state was found in the range of interactions considered.Instead, we observed that the system is always ferrimagnetic: a ferrimagnetic metal below U C and a ferrimagnetic insulator above it.Therefore, our work shows that new physics emerges at fractal dimensions, which is different from the lower and upper integer boundaries.It remains to verify whether this behavior is specific to the Sierpinski triangle, with dimension 1.58, or whether it is generic to other fractal lattices with dimension be-tween one and two.We hope that these findings will stimulate both theoretical and experimental research in interacting systems at non-integer dimension.

A CP-AFQMC method
Consider an initial state |Ψ(0)⟩ that is not orthogonal to the ground state |Ψ 0 ⟩ of the Hamiltonian H. Through imaginary-time evolution, |Ψ(τ )⟩ = exp(−τ H)|Ψ(0)⟩ asymptotically converges to |Ψ 0 ⟩ as τ (a real number) increases.In the Auxiliary Field Quantum Monte Carlo (AFQMC) method, the antisymmetric wave function is expressed as a linear combination of Slater determinants, In our simulations, the coefficients ξ(Φ k ) are not explicitly considered.As imaginary-time evolution proceeds, Slater determinants are either replicated or eliminated.The number of a given |Φ k ⟩ in the sum reflects ξ(Φ k ) [68,69].Typically, the process starts at τ = 0 with all Slater determinants equal to a given trial state |Φ T ⟩, which approximates the ground state and is usually derived from mean-field theories.These determinants are then updated via the application of the imaginary-time evolution operator in a stochastic process.For large systems, direct diagonalization of the Fermi-Hubbard Hamiltonian becomes impractical.Consequently, the Trotter formula [70] is employed to factorize the evolution operator into a product of three terms where K and V represent the hopping and interaction terms, respectively, in the Fermi-Hubbard Hamiltonian.By choosing a sufficiently small δτ , the error introduced by neglecting the O(δτ 2 ) terms in Eq. ( 7) can be minimized to be less than the statistical uncertainty inherent to Monte Carlo calculations, thereby maintaining numerical exactness.The desired limit τ = nδτ ≫ t −1 , where t denotes the hopping strength, is achieved after n successive applications of this small-δτ approximation to |Ψ(0)⟩.A specific iteration on the Slater determinants is described by where the superscript n indicates the imaginary time τ = nδτ .The application of one-body operators on |Φ n k ⟩ results in another Slater determinant.Consequently, exp(−δτ K/2) propagates |Φ n k ⟩.However, since V is a sum of two-body operators, the remaining exponential presents a challenge.To address this, the Hubbard-Stratonovich transformation is employed to convert the two-body interaction into one-body interactions between electrons and auxiliary fields x.We adopt the spin discrete decomposition [71] e −δτ n i↑ n i↓ = e − δτ 2 U (n i↑ +n i↓ ) x=±1 p(x)e γx(n i↑ −n i↓ ) , (9) where p(x) = 1/2 and γ is defined by the relation cosh(γ) = exp(δτ U/2).
For the Fermi-Hubbard Hamiltonian H, this is expressed as where x = (x 1 , x 2 , . . ., x M ) represents a configuration of auxiliary fields, with M being the number of lattice sites, to be sampled within Monte Carlo calculations.Here, p(x) = (1/2) M is a probability distribution function (pdf), and B V (x) is a product of one-body exponentials Direct simulation of Eq. ( 8) is inefficient due to the constant nature of p(x), necessitating an importance sampling technique [68,69].Importance sampling also aids in defining an estimator for system properties and determining constraints that eliminate the sign problem.The importance function implemented is O T (Φ n k ) = ⟨Φ T |Φ n k ⟩, leading to the modified imaginary-time evolution with the modified pdf p(x) = O T (Φ n j )p(x)/O T (Φ n−1 j ).
Since the pdf p(x) is usually not normalized, the normalization factor for each Slater determinant N (Φ n k ) is defined, transforming the iterative projection equation into To manage normalization, weights for each Slater determinant | Ψn ⟩ = k ω n k |Φ n k ⟩ are introduced, updating the weights in each iteration as , with ω 0 k = 1.In practice, the pdf p(x) is sampled by considering each auxiliary field in the configuration x individually.Detailed implementation of the sampling and weight update process can be found in Ref. [33].
The equivalence |Φ n k ⟩ = − |Φ n k ⟩ results in a sign problem that impedes numerical convergence.To mitigate this, auxiliary-field paths are constrained to a region of the configuration space where O T (Φ n k ) > 0, similar to the fixed-node approximation [72].Since the method remains numerically exact if the nodal structure of the trial wave function matches the ground state, this constraint is effective.However, the exact nodal structure of the ground state is generally unknown, necessitating the use of |Φ T ⟩ as an approximation, introducing a small systematic error [68,69,33,73].
Ground-state estimates of the total energy are obtained using the mixed estimator with E n k = ⟨Φ T |H|Φ n k ⟩/O T (Φ n k ) and sufficiently large n.This mixed estimator is exact only if the operator in the numerator of Eq. ( 14) commutes with the Hamiltonian H.
Estimates of other physical observables require the back-propagation technique [68,74,75].The backpropagation estimator is derived from which asymptotically reaches the average value of the observable O in the ground state for τ −τ bp and τ bp ≫ t −1 .The numerical evaluation of Eq. ( 15) is efficiently performed by storing the auxiliary fields sampled during the forward propagation exp(−τ bp H)|Ψ(0)⟩ and using them for back propagation ⟨Φ T |.For a detailed description of this estimator, see Ref. [75].

B Explanation of the six-fold symmetry of the structure factor
To understand the six-fold rotational symmetry of the magnetic structure factor e ik•rij c m (i, j) , as visible in Fig. 24 (a), let us start by noting that our Hamiltonian enjoys a three-fold rotational symmetry.Without loss of generality, let us restrict ourselves to an arrangement of four sites; the generalization to an arbitrary number of sites is then straightforward.The setup is depicted in Fig. 26.To use the symmetry, we introduce the rotation operator R(x, ϕ) = cos(ϕ) − sin(ϕ) sin(ϕ) cos(ϕ) which rotates the vector x counterclockwise by an angle ϕ around the origin (which is here taken to be the location of site 4 in Fig. 26).Explicitly written out, the structure factor of this arrangement reads S(k) = 2c m (1, 2) (A 1,2 + A 2,3 + A 3,1 ) + 2c 1,4 (A 1,4 + A 2,4 + A 3,4 ) + with A i,j = cos(q • r i,j ), and where we have used the fact that c m (i, j) = c m (j, i), and additionally that the setup is three-fold rotationally symmetric.Additionally, due to this symmetry, the following relations (with α = 2π 3 ) hold r 2 = R(r 1 , α), r 3 = R(r 2 , α), r 1 = R(r 3 , α) and r 4 = R(r 4 , α). (18)

FIG. 2 :
FIG. 2: (a) Second generation of the fractal lattice, (b) division of the lattice into two sublattices.

FIG. 3 :
FIG. 3: (a) First generation of the fractal lattice.(b) Its corresponding many-body ground-state energy for different values of Nσ.The x axis in (b) indicates the number of electrons of one kind σ.The total number of electrons in the system is 2Nσ.

6 FIG. 4 :
FIG. 4: Average density configuration of the TB ground-state on the first generation of the fractal lattice with filling (a) Nσ = 1, (b) Nσ = 3, (c) Nσ = 8 and (d) Nσ = 6.Notice the different scale in the colour bar for each individual picture.

1 FIG. 5 :
FIG.5: Spectrum of the single-particle TB Hamiltonian for the first generation of the fractal lattice.The spectrum shows different energy levels: one three-fold degenerate level at zero energy, two two-fold degenerate levels at En = ± √ 2 and two non-degenerate levels.

FIG. 6 :FIG. 7 :
FIG.6: Amplitude scaled to unity of the three basis eigenvectors in the zero-energy level of the TB Hamiltonian, in (a) χ4, in (b) χ5 and in (c) χ6.Dark-green and pink dots on the lattice represent, respectively, negative and positive amplitude of the wavefunction on the sites.On sites with white dots, the amplitude is zero.

FIG. 10 :
FIG. 10: Eigenstates corresponding to the eigenvalues of the two-fold unitary degenerate level CLS type-2.(a) One of the eigenvectors with energy E = −1.(b) One of the eigenvectors with energy E = 1.They are the second type of CLS found by spectral analysis.

FIG. 12 :
FIG. 12: Four out of the seven types of CLS on the third generation of the fractal lattice.(a) One of the type-2 CLS with energy E = −1, (b) one of the type-3 CLS with energy E = − √ 2 (c) one of the type-5 CLS with energy E = − √ 3, (d) one of the type-6 CLS with energy E = −2.149.The dark-green (pink) dots indicate sites with negative (positive) parity, where the size of the dots represents the value of the amplitude: the larger the dot, the higher the amplitude.

FIG. 13 :
FIG. 13: Hisotgram of the density of states in the 8th generation of the fractal lattice.The energy interval is discretized with an energy step of ∆E = 5 • 10 −2 .

FIG. 16 :
FIG.16: Graphical representation of both the Hamiltonian and specific eigenstates for the first and second generation Sierpinsky triangle with intrinsic spin-orbit coupling of strength b.Each circle represents a physical site, with lines denoting coupling between these sites; solid grey lines denote the normal coupling (with strength t), while the dashed lines denote the intrinsic spin-orbit coupling (which is complex and thus directional) of strength ib.The size and color of the circles encode the amplitude of the specific eigenstate shown.(a) The only zero-energy eigenstate for the first generation of the Sierpinsky triangle with intrinsic spin-orbit coupling of strength b.(b) and (c) the two degenerate zero-energy eigenstates for the second generation.Superposing these, as shown in (d) and (e), gives a state with high enhancement at certain sites.

FIG. 17 :
FIG. 17: Graphical representation of both the Hamiltonian and specific eigenstates for the third generation Sierpinsky triangle with intrinsic spin-orbit coupling of strength b.Each circle represents a physical site, with lines denoting coupling between these sites; solid grey lines denote the normal coupling (with strength t), while the dashed lines denote the intrinsic spin-orbit coupling (which is complex and thus directional) of strength ib.The size and color of the circles encode the amplitude of the specific eigenstate shown.(a) to (e) show the five zero-energy eigenstates of the third-generation Sierpinsky triangle (see text for details).

FIG. 18 :
FIG.18: Comparison of the ground-state many-body energies as a function of interaction.We consider the first generation of the fractal lattice with Nσ = 1.These energies are computed with three different implementation methods: ED (green line), MF (yellow triangles), and QMC (red dots).(a) The total energy, (b) the kinetic energy, and (c) the potential energy.The error bars associated to the QMC results are smaller than the pointer's size.

FIG. 19 :
FIG. 19: Comparison of the ground-state many-body energies as a function of interaction in the first generation of the fractal lattice with Nσ = 6.These energies are computed with three different implementation methods: MF (yellow triangles), QMC (U eff = U , red dots) and QMC (U eff = 0.1, green triangles).The error bars associated to the QMC results are smaller than the pointer's size.

FIG. 20 :
FIG. 20: Average density per site ⟨ni⟩ as a function of interaction strength U in the three different groups of states, on the first generation of the fractal lattice with Nσ = 6.

FIG. 21 :
FIG. 21: Local magnetization and average density at halffilling on the second generation with interaction parameter U = 0.1.(a) Local magnetization per site, it shows a ferrimagnetic state; red (blue) arrows pointing up (down) indicate positive (negative) local magnetization.The length of the arrows represents the intensity of the local magnetization.(b) Average density per site of spin-up electrons n ↑ .(c) Average density per site of spin-down electrons n ↓ .The magnitude of the dots is representative of the intensity of the density.

FIG. 22 :
FIG. 22: Average of the local magnetization on the two sublattices for different values of interaction parameter U .

Fig. 23 a
FIG. 23: Doublon density as a function of the interaction parameter U .The behaviour is linear until a critical value of UC ∼ 4.5.

FIG. 24 :
FIG. 24: (a) Magnetic structure factor in the k = (kx, ky) space for U = 6.(b) Radial staggered magnetization (orange stars), radial magnetic correlation (blue triangles) and value of total staggered magnetization (orange horizontal line) for U = 6.The horizontal light-blue line highlights the position of the origin along the y axis.The radial distance r is normalized by a, the distance between two neighbouring sites in the lattice.

FIG. 25 :
FIG.25: Total staggered magnetization as a function of the interaction strength U .For the second generation of the Sierpinski fractal, the calculations were performed using MF (yellow triangles) and QMC (orange crosses).For the regular honeycomb lattice, we extracted the QMC data from Ref.[52] (gray circles).