Robustness of chiral edge modes in fractal-like-lattices below two dimensions: A case study

One of the most prominent characteristics of two-dimensional Quantum Hall systems are chiral edge modes. Their existence is a consequence of the bulk-boundary correspondence and their stability guarantees the quantization of the transverse conductance. In this work, we study two microscopic models, the Hofstadter lattice model and an extended version of Haldane's Chern insulator. Both models host Quantum Hall phases in two dimensions. We transfer them to lattice implementations of fractals with a dimension between one and two and study the existence and robustness of their edge states. Our main observation is that, contrary to their two-dimensional counterpart, there is no universal behavior of the edge modes in fractals. Instead, their presence and stability critically depends on details of the models and the lattice realization of the fractal.


I. INTRODUCTION
The last decade has seen a lot of activity in the classification and realization of insulating electronic states that can be characterized by topological indices. One of the milestones in the field has been the classification scheme known as the 'ten-fold way' or the 'periodic table' of non-interacting insulating states [1][2][3]. The two key ingredients in the scheme are the spatial dimension and symmetries, such as chirality, time-reversal, and parity. For a certain combination of these, it allows to determine whether a topological index can be defined or not, and of which type it is. Importantly, the periodic table does not require the existence of specific lattice symmetries, implying that it also works in the presence of symmetry preserving disorder.
In more recent years, numerous extensions of this classification scheme have been introduced: crystalline topological insulators [4], non-Hermitian Hamiltonians [5], driven non-equilibrium systems [6], or the so called higher order topological insulators [7,8], to name some of the prominent ones.
An alternative path, pursued here, is to extend the classification scheme to allow for non-integer dimensions. The 'tenfold-way' holds only for integer dimension, and an important question is what happens in between. Or more directly, when precisely does a topological state cease to exist? Fractal structures can be defined as possessing non-integer Hausdorff dimension d f [9]. Recently, lattices with fractal features have been manufactured in the laboratory in a variety of ways. This includes the use of molecular assembly [10][11][12][13][14][15][16], templating, and coassembly methods [17,18]. Furthermore, fractal lattices were created by scanning tunneling microscope techniques [19] and with arrays of waveguides [20].
A special place in the 'periodic table' is taken by the class A. In two dimensions, this corresponds to the class of the Integer Quantum Hall (IQH) effect, characterized by the absence of all of the above mentioned symmetries. In the IQH in two dimensions, it is possible to define a bulk topological invariant, the Chern number, which is an integer and directly related to the measurable transverse conductance σ xy [34]. IQH systems are exceptionally robust due to their lack of symmetries. They are not protected by symmetries but instead by the bulk gap. This robustness of bulk properties directly translates to the edge properties. The edges of the otherwise gapped IQH system host one dimensional chiral modes. Each of these edge modes carries one quantum of conductance, e 2 /h, and admits ballistic transport, meaning that they are protected against backscattering from impurities. As a consequence thereof, the transverse Hall conductance is quantized in units of e 2 /h, i.e. σ xy = ne 2 /h. The integer n is related to the Chern number of the bands below the Fermi energy, but also counts the number of protected chiral edge modes in the finite system.
In this work, we study the existence and the robustness of edge modes when two-dimensional models of the IQH effect are transferred to fractal geometries with a Hausdorff dimension between one and two. In practice, this means that we consider lattice implementations of fractals that are embedded in two dimensions and study microscopic tight-binding models belonging to the IQH class. Specifically, we study two microscopic models, the Hofstadter model [35] and a generalized version of the Haldane Chern insulator [36] inspired by Ref. [37]. Both models exhibit stable edge modes in two dimensions, but they differ in one main aspect: the flux pattern in the Hofstadter model corresponds to a physical magnetic field, but the one in the Haldane Chern insulator model does not.
Main result: We study both models on two different fractal lattice implementations, the Sierpinski carpet and the Sierpinski gasket (or triangle). Our main finding is that there is no generic and universal stability of the edge modes, in contrast to their two dimensional counterpart. Instead, the existence of current carrying edge modes strongly depends on microscopic details of both the fractals and its edge construction, as well as the models themselves.
Organisation of the paper: The paper is organized as follows: We introduce the construction principle and the lattice implementations of the Sierpinski carpet and gasket/triangle in Sec. II. Subsequently, we discuss the two microscopic models and their key properties in Sec. III. We then present the computational method to calculate the transport properties in Sec. IV. In Sec. V and Sec. VI, we present results related to Hall transport on the respective carpet and gasket fractal lattice implementations and finish with a conclusion in Sec. VII. The more technical parts are relegated to appendices.

II. FRACTALS AND THEIR LATTICE IMPLEMENTATION
The fractals considered in this paper are defined in the limit of an infinite iteration of a dilution scheme in the continuum that is specific to the fractal. On a lattice, this iteration comes to a natural halt, once the lattice scale is reached. This implies that all the fractals considered in this paper are only approximate fractals, on a scale that is larger than the underlying lattice scale. We study two different implementations of fractals in this paper, the Sierpinski carpet and the Sierpinski gasket. They differ in two important aspects: their fractal dimension and their edge connectivity. Especially the latter plays a crucial role in our studies and is an important factor in determining the stability of the edge modes.

A. The Sierpinski Carpet
The Sierpinski carpet is a fractal that is constructed starting from a simple square through the iterative application of the cutting rule illustrated in the left-hand side of Fig. 1. One first divides a regular square into 9 squares of equal size. Then the central square is removed and one is left with 8 smaller squares surrounding the empty central one. The procedure is then repeated on each of the remaining squares. The full fractal is obtained in the limit of an infinite repetition of this procedure.
Fractal Hausdorff dimension: An illustrative and simple way to determine the Hausdorff dimension of the Sierpinski carpet is to investigate its scaling properties. We consider the scaling exponent of the area under a single cutting procedure. The square count is 8 instead of 9 because the centre was removed. We then proceed to relate this rescaled area to the rescaled length of the outer  square via the effective Hausdorff dimension d H A graphical illustration of the scaling procedure applied to the Sierpinski carpet is presented in Fig. 2.
Lattice implementation: In order to connect the fractal with a tight-binding model, we have to define a lattice version of the above procedure. For obvious reasons, the square lattice provides the most straightforward starting point for this. On a finite lattice, the fractal always has a maximal depth, i.e. a maximum possible number of iterations for the cutting procedure. This is illustrated in the right-hand side of Fig. 1, where we regularize the lattice by placing sites on the centers of the squares in the Sierpinski carpet. Technically, this construction should be named the dual Sierpinski carpet [38], but for brevity we will simply refer to it as the Sierpinski carpet.
Since the maximal depth is related to the size of the lattice, we can parameterize the side length of the square lattice as where a is the lattice spacing. The variable G denotes the size generation, and is the maximal depth to which the fractal can be cut. Given a lattice corresponding to a specific size generation G, we are always free to not cut to the maximal depth and instead stop at any earlier iteration of the cutting procedure. This proves to be a useful strategy for studying the edges. Therefore, we introduce the quantity  F , named fractal generation, which denotes the number of cuts actually applied to our system. It is obvious that F ≤ G holds at all times.

B. The Sierpinski Gasket
The Sierpinski gasket (SG) is constructed starting from an equilateral triangle in a very similar manner. We first connect the centers of all three sides. This forms another equilateral triangle sitting in the inside of the original triangle, but upside down. In a next step we remove it, which leaves us with three triangles. Afterwards, the procedure is iterated on each of the remaining triangles, see Fig. 3.
Fractal Hausdorff dimension: The fractal Hausdorff dimension of the Sierpinski gasket, analogously to the discussion above, can be obtained by considering the scaling exponent of the area with a single cutting procedure according to Lattice implementation: In the SG (left), the underlying triangular lattice is centered on the corners of the resulting triangles, thus, under the cutting procedure, some sites needs to be removed, and some bonds need to be cut, indicated by dashed red lines. Consequently, the remaining triangles share corners. One could also build the dual SG which has the same Hausdorff dimension, but different connectivity and lacunarity. The construction proceeds by assigning a lattice site to the centre of each triangle in the SG.
While both procedures produce the same thermodynamic limit, the effect on the finite systems is very relevant, as the lacunarity is distinct. In Sec. VI, we compare the two different lattice versions. We find that the microscopic choice changes details of the physics, but not the overall conclusion that is drawn in this work.
Note that these are by no means the only ways to implement the Sierpinski gasket on a triangular lattice. One may for instance place the lattice points on both the corners and the edges [10,19], leading to a set of lattices displayed in Fig. 5. This lattice will, however, not be studied in this work.

III. THE MODELS
We consider two prototypical models that both exhibit the IQH effect in two dimensions. Specifically, we consider the Hofstadter lattice model and a version of Haldane's Chern insulator, which we represent on a square lattice for computational convenience.

A. The Hofstadter model
We only shortly review the salient features of the Hofstadter model and defer the interested reader to the existing literature for more details [35,39,40].
In its original version, the Hofstadter model is formulated on the square lattice [35]. Spinless electrons, created (annihilated) by a † i (a i ), hop on a square lattice under the influence of a homogeneous magnetic field that pierces all of its plaquettes. The corresponding tightbinding Hamiltonian reads where the magnetic field is implemented by means of the Peierls gauge connection A ij = e/h rj ri A · d l. Here, L is the set of nearest neighbors with support on the lattice. For concreteness, the magnetic field and the associated flux can be implemented in the Landau gauge A = B(y, 0). We parametrize its strength via B = 2πΦ/a 2 , where Φ is the magnetic flux piercing every plaquette. For future reference, Φ 0 = h/e is the quantum of flux, such that a flux of Φ 0 can be trivially gauged away. The Hofstadter model is most famous for its spectrum versus flux diagram, which reveals the butterfly structure, Fig. 6. Additionally, the gaps have been labeled by their respective Chern numbers, which have been calculated following the TKNN formula [34]. They are in one-toone correspondence with the Hall conductivity. Importantly, by virtue of the bulk-boundary correspondence, the Chern number is equal to the number of protected chiral edge modes.

B. A Haldane type Chern insulator
The Haldane Chern insulator is a paradigmatic model that describes electrons hopping on the honeycomb lattice. Its main feature is that it exhibits Quantum Hall physics without Landau levels, as well as without a net magnetic flux [36]. Here, we study a variant of this model introduced in Ref. [37]. Contrary to the aforementioned paper, we consider a spinless version which explicitly breaks time-reversal symmetry. This is facilitated by complex next-nearest-neighbor hoppings in addition to nearest-neighbor hopping. The model is again formulated on the square lattice and features two degrees of freedom per site, called a and b. The hopping pattern, shown in Fig. 7, does not correspond to a net flux through the plane (the two degrees of freedom are graphically represented as black and red).
In real space the Hamiltonian reads H = H • + H ↑ + H ↓ + H → + H ← + H + H + H + H , and the arrow on each element of the total Hamiltonian indicates the direction in which an electron hops (and agrees with the sketch in Fig. 7). The term H • describes the on-site energy, whereψ is a two-component wave function, whose entries correspond to the two degrees of freedom, i.e. ψ i = (a i , b i ) T and σ x,y,z are the standard Pauli matrices acting in this subspace. The nearest-neighbor hopping terms are given by and the next-nearest-neighbor hopping terms by The index i (j) describes the corresponding site's x-coordinate (y-coordinate). The parameter M is related to the difference in the local potentials between the a and b degrees of freedom, while B (B) is proportional to the (next-)nearest-neighbor hopping between the degrees of freedom of the same type. The absolute amplitude of both the nearest-neighbor and next-nearest-neighbor hoppings between different degrees of freedom are in the following set to B =B = 1, leaving M the only tunable parameter. While the net flux through each unit cell is zero, the model involves complex hopping parameters in a way to cause a quantum anomalous Hall effect (see Fig. 7).
The Chern number C can be computed in momentum space [40]. After imposing periodic boundary conditions and performing a Fourier transformation, the Hamiltonian can be written in the more compact form with and The wave function reads Ψ k = (a k , b k ) T , and σ is the vector composed of the Pauli matrices, i.e. σ = (σ x , σ y , σ z ).
The Chern number for this model can be easily calculated according to C = 1 4π π −π dk x dk yd · ∂ kxd × ∂ kyd , withd = d/| d|. We can identify three topologically distinct regions of the phase diagram, depending on the value of M , i.e. The number of protected chiral edge states is identical to the absolute value of C (the sign of C indicates the direction of propagation) for a given set of parameters. We illustrate the connection between the Chern number and the number of edge modes on a cylinder geometry in Fig. 8. In the trivial phase (C = 0), we see the band structure of an insulating system with no in-gap states. In the C = −1 phase, there are two dispersive in-gap states, one living on the upper boundary of the cylinder, while the other one lives on the lower boundary. The states have opposite group velocities, i.e. the slope of the eigenvalues have opposite signs, which means that the states travel in opposite directions. For the C = −2 phase, we find two additional in-gap states. An analysis of the group velocities and localization reveals that edge states that have the same directionality are also localized on the same edge, meaning they are chiral and there are two of them per boundary.

C. The Hofstadter model on the triangular lattice
This model is formally equivalent to the one discussed in Sec. III A, in that it is described by the same Hamiltonian. The main difference, however, is that the smallest closed loop is not given by a plaquette but by a triangle, and the choice of gauge is slightly more involved. The smaller loop means that the system is not 2π periodic in the phase picked up per plaquette, but is 4π periodic, see Fig. 9 (a). The spectrum of the 'butterfly' as a function of the flux Φ through a lattice is shown in Fig. 9 (b). While it looks different than the one in the square lattice, it has very similar features. It also shows large regions with bulk gaps, where the system is a IQH system characterized by an integer quantized Hall conductance. For all practical purposes, the system behaves similarly to the more conventional square lattice version. The main difference will be where we put the leads, as detailed in Sec. VI. In Figs. 9 (c) and (d), we contrast the spectra of the SG and the dual SG both from the triangular lattice, as well as from each other. The main observation is that the spectra of SG and dual SG are markedly different as a consequence of the two distinct cutting procedures.

IV. THE METHOD
We use two complementary methods to study the properties of the edge states: eigenstate spectroscopy and the non-equilibrium Green function technique. The former one is straightforward, so we do not comment on it any further. The latter, on the level it is used here, is equivalent to the Landauer-Büttiker approach. We implement it numerically and in order to achieve larger system sizes, we combine it with the recursive Green function method. The setup we study for square based geometries is shown in Fig. 10. It also applies to the triangular systems, although with slightly modified lead locations. Four leads are connected to the central scattering region (S) that hosts the fractal. The leads can in principle sit at different temperatures T a and at different chemical potentials µ a (a = 1, 2, 3, 4).
The main advantage of this setup is that it allows to study both the diagonal voltage drop V xx and the Hall or transverse voltage drop V xy , which provide complementary insights.
Following a standard Keldysh Green function treatment, one can derive the following expression for the cur- Figure 11. The Hofstadter model phase diagram, in terms of ρx,y for a full depth fractal of a generation 5 Sierpinski carpet. Picture from Ref. [22]. rent into lead a: where is the transmission function between leads a and b. The trace in Eq. (12) extends over all internal indices, such as for instance orbitals and the sites along the length of the interface. G are the retarded (advanced) Green functions describing propagation from lead a to lead b and vice versa (these are generically large matrices that have to be handled numerically). The broadening function Γ a (ω) describes the hybridization function of lead a with the scattering region S. It is given by where G a is the lead Green function, meaning the Green function of the decoupled lead at the interface. This quantity is also calculated numerically. Furthermore, f a denotes the Fermi-Dirac distribution within lead a, where β a = (k B T a ) −1 describes the inverse temperatures and µ a is the chemical potential of the respective lead a.
In the following analysis, we always consider zero temperature in the leads, implying that the distribution function is the step function. Consequently, the current into lead Figure 12. Zoom in on the Hall resistance for three consecutive generations (4, 4), (5,5) and (6,6). The large scale structure of the edge modes is already present in generation (4,4). The difference between the various panels is the amount of high frequency variations, as seen in the two cuts at fixed Φ/Φ0 in the lower panel. The color scheme is the same as in Fig. 11. Picture from Ref. [22].
a reads where in the last line we assume µ a(b) = µ + V a(b) and T a,b (ω) ≈ T a,b (µ). In this paper, we are mostly concerned with calculating the transverse resistance and, to a lesser extent, the longitudinal one. The pattern of currents that allows us to access both is given by I 1 , I 2 , I 3 , I 4 = I 0 , −I 0 , 0, 0. If we ground e.g. lead no. 4, i.e. V 4 = 0, we can determine all the potentials and the Hall resistivity is given by The main advantage of the above formulation in terms of Green functions is that we are able to express the transmission in terms of only a very limited number of elements of it. This allows to use a very efficient numerical technique, the recursive Green function method. The algorithm has the benefit of reducing computation time significantly. Instead of inverting the whole central scattering region S, one can do it iteratively, slice by slice.
However, it needs to be adapted to accommodate the presence of more than two leads. Technical details on the Recursive Green function formalism and slicing procedure are given in Appendices A and B for the interested reader.

V. EDGE STATES AT THE BOUNDARY OF THE SIERPINKSI CARPET FRACTAL
In this section, we discuss the physics of the edge states on the Sierpinski carpet of both the Hofstadter and the Chern insulator model. Our main finding is that the stability of the edge states of both two dimensional models differ significantly upon rendering the lattice fractal.
In a strictly two dimensional setting, the Hall conductivity σ xy is quantized according to σ xy = Ce 2 /h, where C is the Chern number of the underlying electronic structure. Alternatively, the transverse conductivity can be considered from the point of view of chiral edge modes in the system, following the bulk-boundary correspondence. It is important that the bulk itself is gapped. Then, the Hall conductivity is quantized as |σ xy | = ne 2 /h where n is the number of protected chiral edge modes. Since the Chern number is only properly defined and quantized for gapped systems in two dimensions (or more generally, in even dimensions), we will henceforth not talk about the Chern number but only about σ x,y (and ρ x,y ), since this quantity is always well defined.

A. The Hofstadter model on a Sierpinski carpet
In this part, we investigate the generalization of the Hofstadter model to a fractal geometry. Some of the results have already been published elsewhere [22], and we only give a very condensed version here. The main results are summarized in Fig. 11 and Fig. 12. Similar results have also been obtained in Ref. 26 using the Kubo-Bastin formula for Hall conductivity.  Fig. 6, the number of regions with perfect quantization is massively reduced.
It was found that the vital quantity that determined the maximum number of edge modes was not G or F , but the difference between them, the 'fractal distance' We will see that this quantity plays a crucial role also in the systems that we consider in this work. This becomes more apparent in Fig. 12, which shows a zoom-in into the upper panel of Fig. 11. From left to right, we increase the fractal generation (all at full depth, ∆ = 0), making sure that the results have converged. In the lower panel, we see a corresponding cut at fixed flux, again for different generations.
For some fluxes we find plateaus, whereas for others no quantization is visible. This suggests that generically the quantum Hall physics of two dimensions is unstable to modifying the dimension. From the point of view of Chern numbers and the periodic table, this is not unexpected [1-3]. It is not obvious how to define Chern numbers if the dimension is not two, and therefore one could expect that for fractals the quantization, in general, is gone.
The specific problem with fractals is that even if the embedding space is two-dimensional the system does not possess any finite period, which prevents the formation of a two-dimensional Brillouin zone. One can artificially impose periodic boundary conditions on any finite-size fractal. However, these "periodic" systems would have just as many bands as there are sites in the fractal, and this number would grow with system size, preventing the number of bands from being a well defined quantity in the thermodynamic limit.
Alternative methods are currently being developed to describe the topology in systems lacking translation invariance. These include calculations in real-space [43], perpendicular space [44] or using non-commutative geometry techniques [45][46][47]. The real-space calculations on Sierpinski fractals have been attempted in Refs. 21 and 22, but as far as we are aware, these do not guarantee quantization.
Since the quantization of the Hall conductivity is also related to the existence of chiral edge states (assuming there is a bulk gap), this begs the question of what happens to them.
A priori nothing forbids the existence of edge states in a fractal between one and two dimensions. In order to analyze this, we consider the modes that correspond to states at the Fermi level in one of the plateaus shown in Fig. 12. We choose a flux of Φ/Φ 0 = 0.375 and identify edge states at all depths of a fractal of generation 4, see Fig. 13 (upper row). Given that edge modes can survive on the exterior edges, it is natural to ask whether the same is true also on the interior edges, which are formed by the holes that are cut away from the fractal. Indeed, as shown in the middle panel of Fig. 13, one may identify states that are localized on the inner edges as well (note that we have added weak disorder in all plots to remove accidental degeneracies due to lattice symmetries). These states are also identified in Fig. 14 via a Hall resistance measurement involving leads placed on the inner edges instead of the outer ones. In the figures, the leads are placed in the central, F = 1, "hole" and ρ x,y is computed for ∆ = 0, 1, 2 and G = 3, 4, 5. We find, again, that ∆ determines whether or not edge modes are stable on the inner edge, just like on the outer ones. The key to understanding the stability and the associated quantization plateau is that for some energies the wave function is localized on the outermost boundary sites. This is a fine tuned situation and it renders the system less stable against adding disorder, as some of us showed in Ref. [22]. For reference, a number of typical bulk states are shown in the lower panel of Fig. 13 for the same cutting depths.

B. The Haldane Chern insulator on a Sierpinski carpet
We now turn our attention to the Haldane model defined in Eq. (7). The setup we consider is of the type shown in Fig. 10. Our starting point is the square lattice in the scattering region S, to which we then apply the cutting algorithm to generate the Sierpinski carpet. We first investigate the Hall resistance for signatures of quantization, in the sense of the IQH effect.
Starting from a two-dimensional system of size G = 5, we compute the Hall resistance for different cutting depths, i.e. for F -values between 0 and 5. The results are shown in Fig. 15.
A low number of cuts leaves the Hall resistance largely intact. However, starting from cutting depth F = 3, the quantization appears less stable and completely vanishes for F = 4 or F = 5. We will discuss later that this is related to the space accessible to a boundary mode.
This behavior is also observed for other system sizes, as can be seen in Fig. 16. We find that the relevant measure for the breakdown of the quantization is not the fractal  generation G itself but the "fractal distance" ∆ = G − F . If ∆ > 2, the quantization is intact and the Hall resistivity remains unchanged. For ∆ < 2, no well-defined region with a finite Hall current exists. This holds for regions with one edge mode, as well as regions with two edge modes in the non-fractal system. At ∆ = 2, there are still regions with a finite Hall resistivity, however they show features of instability.
An interpretation for ∆ can be found by considering the boundary of our system. Since ∆ refers to the difference in number of actual cuts made compared to the maximum number of cuts possible, it measures the width of the system boundary that is still intact. The relation between the number of sites that have not yet been touched by the cutting procedure, i.e. the boundary width b and ∆ is given by This relation holds independently of the size of the full system. A graphical representation for some values of ∆ is given in Fig. 17.
Since the quantized Hall current is carried by edge modes, details of the edge along which they are traveling should indeed have an impact on their stability. Understanding why for any boundary with a width of 9 sites or less (∆ ≤ 2) the Hall resistance becomes unstable requires a closer look at the edge states themselves.
We consider a cylinder geometry, like the one shown in Fig. 18, and investigate the spatial dependence of the wave function transverse to the edge direction. By performing a Fourier transform on only the x-component, we may study the localization of the edge modes in the finite y-direction. The resulting Hamiltonian is given by Δ=1 Δ=2 Δ=0 Figure 17. The boundary width b is dependent on the cutting depth (colored regions). Independent of the size of the fractal, for any given ∆, b remains the same.
W y x Figure 18. The cylinder setup with width W used for calculating the wave functions of edge modes. In the x-direction the system has periodic boundary conditions.
and W is the width of the ribbon shown in Fig. 18. In Fig. 19, we show the edge states on the non-fractal square lattice for a system size corresponding to generation G = 5 (F = 0). This serves as reference point for the following discussion. We show the absolute value of the wave functions as a function of the distance to the edge for the first 25 sites. There is no visible difference between the shown curve and that of G = 4, from which we conclude that for the system sizes considered, finite-size effects are not important. We now compare the typical extension to the length scale defined by the cutting procedure.
Vertical lines indicate at which points a cut would effectively terminate the intact boundary of a fractal system for a given ∆. For ∆ = 0 (∆ = 1) this implies that there are already missing sites in the second (fourth) row of atoms of the system. What this means in practice, is that a sizeable weight of the wave function will be localized around holes, which leads to situations where the different parts of the wave function are effectively counterpropagating, see Fig. 19 (d).
Such systems can therefore not sustain stable edge modes and it is not surprising that the quantization of the Hall conductivity breaks down in these cases, as shown in Fig. 16. If, however, the boundary is larger than ∆ = 2, the weight of the wave function is far enough from any holes in the sample. Therefore, for ∆ ≥ 2 the Hall resistance in Fig. 16 is hardly influenced by the presence of cuts. The case ∆ = 2 marks the threshold between both behaviors. There, the main peak of the wave function still gets supported by the lattice. However, already the second peak in Fig. 19 (a),(b) cannot completely fit into the boundary when the fractal gets cut out. This explains why, in the Hall resistance plots, we can still find an area with well defined edge modes, but in combination with additional substructures that come from the presence of another cut-out edge hybridizing with the outer modes.
Further insight can be gained by considering the stability of edge modes on the fractal under the influence of dis- order. We implement disorder via a random on-site potential , which is chosen in the range ∈ [− max , max ]. The results are shown in Fig. 20 for the parameter value M = 4. Since edge modes were already shown to be unstable for ∆ < 2, these cases will not be considered. The results in Fig. 20 were calculated for G = 5 fractals, but smaller systems show the same qualitative behavior. As in the two-dimensional case, i.e. F = 0, there exists a critical disorder at which edge modes cease to exist. This critical threshold also exists in the fractal system. However, the value for such a threshold becomes slightly smaller as ∆ decreases. At ∆ = 2, we see some features already for small disorder, but in general there still seems to be some stable Hall current. This can be explained by the finite size of the edge modes making them unstable on smaller boundary sizes, as was previously discussed.
We end this section with a comment. The reader might wonder if there is a different explanation for the breakdown of topology that does not rely on the mechanism of inner and outer counter-propagating edge modes being gaped out by the presence of the fractal cuts. For instance, the perturbation introduced by fractal formations could affect different single-particle states differently. One could then imagine that the periodicity that is introduced by a certain fractal depth happens to be commensurate with electron states that bear the largest Berry curvature Ω(k). Those states may then be strongly scattered by the newly added fractal perturbation and may open up trivial gap due to Weiss oscillations [48,49]. It is thus justified to wonder what the fate is of the edge modes if one instead of fractal cuts perform periodic cuts of the same size as the cuts at a certain "fractal-distance" ∆.
We have not pursued this direction in this work, but it makes an interesting follow up study. One may argue that if the hybridization of edge modes is not important, then topologically non-trivial band structures would persist even at "periodic" cuts of ∆ = 1 or ∆ = 0, for at   Fig. 21 a). Main panel : The Hall resistivity, ρx,y. Inset: Transmission from lead 1 to lead 3, T1→3. For Φ/Φ0 = 0.07, bands of Chern numbers up to 5 can be identified, with the corresponding Hall voltage being quantized at ρxy = 1/5h/e 2 . For Φ/Φ0 = 0.3, a large energy gap with one edge mode is observed in the energy range 0 < E < 3.5.
least some range of M . One could (but we have not) test this hypothesis by making periodic cuts of finer and finer size, and compare with the fractal cuts. The type of "periodic" deformations described above would also allow for a multi-band Brillouin zone to still exist, and if the number of bands is reasonably small, one could still compute Chen numbers in the infinite system.

VI. THE SIERPINSKI GASKET
We return to the Hofstadter model, but on an underlying triangular lattice and the Sierpinski gasket geometry.
Since we are interested in the Hall voltage, this begs  Fig. 21. In setup a), VH J has trouble picking up any signals of Hall quantization. In setup b), some quantization can be seen in VH , but it does depend on the lattice regularization. We conclude that the lead setup in Fig. 21 b) is preferable over setup a) for the purposes of detecting a quantized Hall response.
the question of how to define the Hall voltage. In this geometry, we adapt the "standard" Hall measurement to the triangular setup by putting two of the leads on the same side. The measurement then proceeds in the standard way, by running a current between "opposite" leads and measuring ρ x,y on the remaining two leads (Fig. 21).

Lead placement and shallow cuts
Here, we explore two different placements of the leads in order to stabilize the Hall voltage measurement, see Fig. 21. In setup (a), leads 1 and 3 are placed symmetrically over the pinching points that appear already in the first fractal generation G = 1. In setup (b), the leads 1 and 3 are placed asymmetrically next to the G = 1 pinching points, but centered on the pinching points for the second fractal generation G = 2.
As a reference, we study the non-fractal system first. The results are displayed in Fig. 22, where we show the Hall voltage at low field Φ/Φ 0 = 0.07 and high field Φ/Φ 0 = 0.3. In the low-field regime, Landau bands with up to Chern number 5 can be identified by the quantization of ρ x,y , whereas in the high-field regime, the spectrum is dominated by the large gap 0 < E < 3.5, with only a single edge mode. Note that in the non-fractal case the two different lead placements in Fig. 21 lead to equivalent results (not shown).
The SG differs substantially from the SC in that already at the shallowest cut F = 1, the edge contains pinching points which are only one lattice site wide. A direct consequence is that these can only allow for one edge mode to pass and this limits the maximal transport through the pinching point to be at most 1. This also means that we should only expect a quantized Hall voltage that is either ρ xy = 0, ±1.
In Fig. 23, we investigate whether there is a preferred way to place the leads to detect this one edge mode. We here focus on the case Φ/Φ 0 = 0.3 and compare the ρ xy for both the SG and dual SG regularizations with respect to the placement of the leads in Fig. 21. For this purpose, we only make a single cut F = 1 to try and detect the possibility of a single edge mode. For the "symmetric" placement, we see that ρ x,y shows strong high-frequency fluctuations precisely in the region where a quantized ρ x,y response could be expected. We speculate the these high frequency oscillations are due to the symmetrically placed leads acting as strong impurities and interfering with the path of the edge mode, as it tries to navigate the pinching point which acts like a point contact.
On the other hand, in the "asymmetric" lead placement (lower panel) a clear Hall plateau is observed for the SG system. The same can, unfortunately, not be said for the dual SG system in this energy range. We note however that in the range −1.5 < E < −1 both the SG and dualSG lattices show a quantized Hall response. This shows that at least for a shallow cutting, both lattices are able to support edge modes. We conclude that the "asymmetric" placement of the leads is preferred for detecting edge modes, and we will thus only use that one in what follows.
We note that a further consequence of the pinching points at F = 1 is that "bulk" and "edge" currents will need to pass though the same point. This will lead to the possibility of mixing between the transverse and longitudinal resistivity. Indeed, the oscillatory behavior of the Hall resistivity shown in Fig. 23 is similar to the behavior of Hall measurement in 2D electron gasses that occurs when partially filled Landau levels cross the Fermi level as the B-field is increased. In the current setup, longitudinal resistivity cannot really be measured, as this would require a 6 terminal setup. Our computational method can easily be modified to include also a 6 terminal setup, but it starts to become computationally expensive, and would introduce another layer of finite-size effects, so we choose not to do it in this work. Thus, we cannot conclusively say that the longitudinal and transverse resistances are not mixed.

Fractal cuts
Next, we perform cuts at maximal depth, such that the fractal depth is F = G − 1. The result is displayed in Figure 24. The Hall voltage for Φ/Φ0 = 0.3 and Φ/Φ0 = 0.07, and F = G − 1 for G = 5, 6, 7 with the lead setup shown in Fig. 21(b). In each panel, only a segment of the spectrum is shown. The value of ρx,y is shifted by 0.1 between the three generations for increased readability. In addition, values of ρx,y that deviate from h/e 2 are made successively whiter to suppress noise in ρx,y and to highlight the regions where ρx,y is quantized. Note that for all three system sizes roughly the same behavior of ρx,y is observed and that the different system sizes share the same regions with quantized ρx,y. It is thus reasonable to suspect that these plateaus will be present also for larger system sizes.   In the panels of Fig. 26, we can see that for all three system sizes, roughly the same behavior of ρ x,y is observed. Also, we note that the different system sizes share the same regions with quantized ρ x,y . It is thus tempting to conclude that the fractal is able to support topological states in the thermodynamic limit. Indeed, by diagonalizing the Hamiltonian without leads, we find clearly distinguishable edge modes in the full depth fractal. These are depicted in panels (a) and (b) of Fig. 25.

Interior transmission
Interestingly enough, from the point of view of the pinching points, there is no particular reason why the edge states of the Sierpinski gasket have to run only on the outside of the gasket. In fact, the interior of the gasket forms an edge that is just as valid as the exterior edge. Indeed, in Fig. 25 (c) and (d), one can see edge states at the interior of the gasket.
If we place leads around the inner holes of the main triangle, as depicted in Fig. 26, we can clearly see regions of quantized Hall resistivity. We note two things: (i) the regions of quantized resistivity are the same for the inner and the outer placement of the leads. (ii) here ρ x,y = −h/e 2 instead of ρ x,y = +h/e 2 , as in the outer placement.
From the second observation, we conclude that for the inner placement the edge modes are propagating as 1 → 2 → 4 → 3 → 1, whereas in the outer placement the circulation direction is 1 → 3 → 4 → 2 → 1, as indicated by the arrows in Fig. 25.
Combining the second observation with the first observation that the interior and exterior quantization happens in the same energy range, we draw the conclusion that the pinching points of the gasket are not enough to gap out the edge modes. Rather, it looks as if the pinching point is alternating between sustaining propagation on the exterior and interior edges, allowing these in practice to coexist for the purposes of transport.
The observation of coexisting transport can be understood from the perspective of a quantum Fabry-Perrot interferometer, as shown in Fig. 27a). In an interferometer of this kind, edge modes can circulate on each of the regions I-III without being destroyed by the quantum point contacts connecting them. At the same time, currents may tunnel between region I and II (or II and III) coherently. An analogous situation can be found in the case of the Sierpinski gasket with F = 1. The three triangles that are formed at F = 1 can alternatively be thought of as forming a Fabry-Perrot interferometer necklace, Fig. 27 (c), or as a pinched Corbino geometry. In both cases, the counter-propagating edge modes do not gap each other out completely allowing for coexistence and a quantized Hall resistivity response on the inside, as well as the outside boundary.

VII. CONCLUSION AND DISCUSSION
In this paper, we investigated the stability of the IQH effect when paradigmatic models are transferred onto lattices that implement fractals with dimensions between one and two, specifically the SC, SG, and the dual SG. On general grounds, starting from bulk considerations, one should not expect a robust quantization of the transverse conductance: The Chern number is only well defined in two spatial dimensions.
The question can also be investigated from the point of view of protected edge modes. They are known to exist in two-dimensional IQH systems by virtue of the bulkboundary correspondence. We study how their stability is compromised by rendering the lattices increasingly fractal. We use a spectroscopic method and the Green function method.
Our main finding is that under generic circumstances, the edge channels become unstable and no quantization can be expected. The main reason for the instability is rooted in the fact that the fractal introduces new edges on the inside of the sample. The associated 'inner edge states' are counterpropagating with respect to the outside edge states and can hybridize. This eventually gaps them out. The only exceptions we find are in situations where either the edge states are extremely localized or where the fractal does not extend all the way to the edges due to the cutting depth. We show that in those situations, the protection against disorder is also strongly reduced compared to the two-dimensional counterpart.
We stress, however, that even though a robust quantization of the transverse conductance is not expected, there are fine tuned situations where it still exists. In these situations, many of the topological features that are observed are also scale invariant, in that they only depend on the fractal distance, ∆ = G − F , and not on G or F individually.
For future, it would be interesting to investigate other classes of topological insulators on fractal geometries, and also the connection with Weiss oscillations [48,49] on periodic structures with the same length scale as the fractal cuts. means that using exact diagonalization, by directly solving Eq. (A1), in order to calculate the Green functions is not feasible for larger systems.
Fortunately, the transmission function Eq. (12) does not depend on knowing the full Green function matrix. Instead, it is enough to know the part of the Green function responsible for propagation from left to right together with the lead contributions.
In order to visualize the dimensionalities of different parts of the Green function, we can write it as a matrix of matriceŝ the interpretation of which is depicted in Fig. 28. Note, that while the amount of matrix elements in the whole system grows with the volume of the scattering area, the corner elements G L,R , G L,R that describe propagation from one lead to another only grow as the surface of the leads attached on the sides. The same holds for the other corner elements G R,R and G L,L , which will be required to compute the lead contributions Γ L and Γ R in Eq. (13). Therefore, if there is a way to compute them without solving the whole system, we will be able to compute the transmission function (12) for much larger system sizes. In the following, we will show that this is indeed achievable using a recursive approach, a flowchart of which is shown in Fig. 29.
Since we are interested in the transport from the left to the right lead, the starting and ending point of the recursion are always given by the leads themselves. Thus, the first Green function to be computed is always the surface Green function of the left lead G L , before it is attached to the scattering region S. For a semi-infinite lead described by e.g. a simple tight-binding chain, G L can even  Figure 29. A flowchart depicting the steps taken in the recursive Green function algorithm. Each step is described in the text in Sec. A.
be solved exactly, but in more complicated situations G L needs to be determined numerically. Each slice contains the sites to which the system from the previous step directly connects, and that were not considered before. In a simple square geometry, this implies that the cuts are equivalent to taking all sites that share the x-coordinate x = j + 1 if the previous slice ended at x = j. This remains true even if a fractal is cut into that square geometry. An example of a simple linear slicing is shown in the left panel of Fig. 30.
For each slice, interaction matricesV ← andV → are determined by taking the appropriate terms out of the Hamiltonian that connects both slices. For the Haldane model in (6), this implies that if we attach a slice that connects sites with x-coordinate j + 1 to the previous slice, where x = j, the interaction matrices are symbolically given byV where it is understood that only terms containing elements on both collumn j and j + 1 and included. Note that on a fractal geometry, some of these Hamiltonian elements may not exist due to missing lattice sites. This implies that the amount of sites that connect to a slice, as well as the size of a slice itself are not constant, and therefore the dimensions of the interaction terms and Green functions vary. If we let G (n) j,k denote the Green's function from slice j into slice k after n slices have been attached, and H n denote the Hamiltonian of the nth slice, the recursive algorithm reads We may take G 0,0 = G Lead to denote the initial semiinfinite lead. Note how G S and G T get updated as more slices are added to the system, whereas G (n) j,k refers to a specific setup.
In the last step of the algorithm n = N + 1, the right lead needs to be attached to the scattering area. In this step, the bare Green function of the new slice gets replaced by the lead Green function again, such that Note that the transmission (12) as well as the recursive algorithm (A3) require the lead Green's function G Lead . Since it is needed as a starting point for the recursive algorithm, it has to be found before applying the recursion method. However, computing it is closely related to finding the surface function G (n) n,n in the recursion Eq. (A3) with the addition that, for a semi-infinite system, adding another slice does not change the surface Green function. This implies that we need to solve the self-consistency equation where here, H S refers to the Hamiltonian on the surface of a semi-infinite lead. In this work, the leads are assumed to be simple tight-binding chains with nearestneighbor interaction, which in the case of Haldane model do not couple the two orbitals of the scattering area. Thus here, H S is a one-dimensional tight binding chain with a length equal to the width of the scattering area.

Appendix B: Recursive Greens function slicing for four leads
In this work, we are studying a scattering area with more than two leads. As a result, the slicing procedure that was discussed briefly in Appendix A needs to be adjusted.
In the treatment of the leads, the direction in which they extend infinitely was integrated out exactly, such that they can be coupled to the scattering area as an effective self-energy. On the inverse Green function level, this self-energy affects the sites that directly couple to the leads only. Thus, in order for the recursive method to work, we need to be careful when we add a slice of our scattering area that is connected to an external lead. To be precise, we need to slice our system such that each lead is attached to one, and only one, slice. In Fig. 30 we give an example for a slicing procedure in a four-lead setup, in direct comparison to the slicing of the same scattering area for a two-lead system. Furthermore, the simple fact that we consider currents between more than two leads, as shown in Fig. 10, means that the target sites for which we calculate the transmission Green function will not always sit in the very last slice. Thus, when a slice labeled by index n is attached, in addition to the computation of the surface Green's function G S = G