Effects of the Vertices on the Topological Bound States in a Quasicrystalline Topological Insulator

The experimental realization of twisted bilayer graphene strongly pushed the inspection of bilayer systems. In this context, it was recently shown that a two layer Haldane model with a thirty degree rotation angle between the layers represents a higher order topological insulator, with zero-dimensional states isolated in energy and localized at the physical vertices of the nanostructure. We show, within a numerical tight binding approach, that the energy of the zero dimensional states strongly depends on the geometrical structure of the vertices. In the most extreme cases, once a specific band gap is considered, these bound states can even disappear just by changing the vertex structure.

The physics is similar in the two cases: The systems have a gapped bulk, and bound states located at the physical boundaries with energy in such gaps. The bound states are predicted to enjoy, for topological reasons, enhanced robustness with respect to perturbations and are hence promising for quantum technological applications.
To make progress, new paradigms for topological protection have hence recently been explored. A very promising one is higher order topology [64][65][66][67][68]. A d-dimensional higher order topological insulator or superconductor is a system in which the topologically protected states live in n dimensions less than the gapped bulk, with 1 < n ≤ d. The index n is referred to as the order of the topological insulator. Very recently, three-dimensional gapped systems with one-dimensional topological states have been realized experimentally [69,70]. On the other hand, proposals for the highly relevant case of two-dimensional bulk and zero-dimensional topological bound states have been put forward, but not experimentally accessed.
Among such proposals, a noteworthy one has been envisioned by S. Spurrier and N. Cooper [71]. The main ingredient is a quasicrystalline two-dimensional lattice, generated by two layers of graphene-like lattices with a 30 degree relative rotation, and possessing 12-fold rotational symmetry. Each layer consists of a Haldane model and the two layers have opposite Chern number. Moreover, the edges of the bilayer structure they consider are of the armchair and bearded type. The result is that the system, thanks to the presence of the rotational symmetry, hosts topological zero-dimensional bound states, which are degenerate, isolated in energy from the others, and localized at the vertices of the structure. The system is hence a higher order topological insulator with d = 2 and n = 2. A subsequent work [72] showed that these properties also characterize samples with armchair-zigzag edges. However, in this case, the energy at which the gap opens drastically changes. While in the armchair-bearded case the gap opens at the charge neutrality point, in the armchair-zigzag case it is substantially moved away from it. This observation poses questions about the sensitivity to details of the model.
In this context, the aim of the present work is to address some questions related to the robustness of the properties of the topological bound states with respect to lattice perturbations at the corners of the structure: Do these corner modes always exist in a given gap? Which is their energy? How strong is the protection granted to their degeneracy?
We will show how the energy, degeneracy and even existence of the predicted corner modes astonishingly depend upon the precise disposition of the lattice sites around the sample vertices. These results seem to narrow the path for a hypothetical experimental realization of the system proposed, which would require site-wise precision in the realization of the sample in order to achieve the HOTI phase close to the Fermi level.
The paper is organized as follows: In Section 2 we introduce the model, both for how it concerns the lattice and for how it concerns the Hamiltonian. In Section 3 we present and discuss our results, answering one by one the questions we addressed above. In Section 4 we draw our conclusions.

Model
In order to proceed with the discussion of the results achieved, we must first introduce the model of HOTI on a quasicrystal, which will be the object of our analysis. The present section is organized in two parts: in the first one we introduce the system lattice, while in the latter we present the model Hamiltonian. In what follows, we will consider a square sample instead of a dodecagonal one, as is the case in [71]. Indeed, when dealing with the robustness of the bound states, the choice of a square sample provides major computational advantages, while leaving the physics qualitatively unchanged.

Lattice
Let us begin by describing the construction of a finite-size square sample of the quasicrystalline lattice under inspection: one starts with the two honeycomb layers superimposed with AA stacking (i.e., each atom of the top layer is right on top of an atom of the bottom layer) and rotates one of the two layers by thirty degrees with respect to the other, taking a rotation axis that goes through the center of two superimposed hexagons (see Figure 1). The resulting lattice is crystallographically equivalent to that of the graphene quasicrystal, which has been shown to possess local bulk C 12 rotational symmetry [73].
We then crop our lattice with a square shape centered in the rotation center, and oriented in such a way that the edges are not jagged. As already noted in [72], the edges of the resulting system-which possesses C 4 global rotational symmetry-can only be given by two possible pairs of the upper and lower layer edges: bearded-armchair or zigzag-armchair. However, even once the combination of the edges has been chosen, there are still two families of square lattices that can be obtained: these differ with each other for the disposition of the lattice sites around the vertices. The zoom of a vertex for each of the four possible configurations just discussed is shown in Figure 2. In [72], it was shown how the properties of the model for an HOTI on the (spinless) graphene quasicrystal first proposed in [71], depend on the shape of the edges. Here, we want to discuss how, even when the edge combination is fixed, the features of the HOTI phase can still depend on the specific configuration of the sites near the vertices. In order to approach this target, let us first introduce the model Hamiltonian.

Model Hamiltonian
The HOTI model we analyze [71] is governed by the following tight-binding Hamiltonian where we denote i creating a fermion on the i-th lattice site of the top/bottom layer of honeycomb lattice. The Pauli matrices τ i act in the space of layer pseudospin, coupling the Fermionic operators of the two layers. The Hamiltonian H presents two in-plane hopping terms: the first one describes nearest neighbors ( ij ) hoppings with amplitude t, while the second one involves next-nearest neighbors ( ij ) hoppings of amplitude λ H . The latter is associated to a complex phase iν ij (ν ij = ±1), given by the usual Haldane coupling [74,75]. Also, an inter-layer hopping term is present, parameterized by t ⊥ ij . Its explicit form is borrowed from tight-binding models for real twisted bilayer graphene (TBG) samples [76,77].
with d ij being the vector that goes from the i site on one layer to the j site on the other one. We take the the parameters appearing in Equation (2) as in tight-binding models for real TBG samples as well [76,77]: δ = 0.184a is the decay length of the transfer integral, d 0 = 1.362a the distance between the two layers and t ⊥ = −0.178t a coupling coefficient. All lengths are expressed in units of the in-plane next-nearest neighbor distance a. Finally, with the aim of tuning the interlayer coupling strength, we add the parameter λ ⊥ in Equation (1) (the parameter t ⊥ is taken in such a way that if λ ⊥ = 1, then the value of the coupling for interlayer vertical hoppings matches the value used in realistic tight binding models for TBG [76,77]).
In order to achieve a slightly more intuitive perspective about the model, it is worth noting [71] that without the interlayer hopping term (λ ⊥ = 0), the Hamiltonian actually consists of two Haldane models [74] with opposite Chern number (±1 in the upper/lower layer respectively) and with a relative 30 • twist. This is almost the same as the Kane-Mele model [7,8], just with the spin mapped into layer pseudospin. Indeed, with λ ⊥ = 0 the low energy spectrum would be gapless, with the edges of the two planes hosting a pair of counterpropagating modes-just as in the Kane-Mele model. What the interlayer hopping does is to couple together these two Haldane models: more specifically, as was shown in [71], it gaps out the counter-propagating edge modes on the two layers, making it possible for the second order topological insulator (SOTI) phase to occur. Such a phase is characterized by the presence of bound states localized at the corners of the system (corner modes), degenerate among themselves and with energy laying inside the edge-gap.
The phenomenology just discussed is encoded in the effective low energy theory of the model derived in [71]. This theory is only based on the symmetries of the model and it leads to a few possible values for the energy of the corner modes inside the energy gap. For the present case of square sample (C 4 global symmetry), these are given by: Having introduced all the necessary ingredients, we can now proceed to presenting and discussing our results.

Results
We will now show the results obtained through the numerical diagonalization of Hamiltonian (1) on different square samples of the quasicrystalline lattice. For the construction and diagonalization of the model, the python package pybinding [78] has been used. The fixed parameters of the Hamiltonian were set as described above. We set λ ⊥ = 2, resulting in interlayer hoppings of double strength with respect to real TBG samples. This choice does not change the physics on a qualitative level, but it increases the gap width, reducing the hybridization of the corner modes for the same edge lengths. Also, we set λ H = 0.3t, so that the Haldane bulk gap is maximized [74].
The present section is organized with the aim of separately answering the questions addressed in the introduction about the stability of the corner modes with respect to local geometrical perturbations of the lattice. We start by addressing the issue of the existence of the corner modes in a certain gap.

Existence
Let us begin by presenting the low energy eigenvalues obtained for four different square samples, one for each of the categories in Figure 2. These are reported in Figure 3, where the plots of the eigenvalues are ordered the same way as the corresponding lattices, exemplified in Figure 2. Figure 3a,b reports the low-energy eigenvalues for two square samples with beardedarmchair edges and different shapes of the corners (see Figure 2a,b); in both cases we find four degenerate eigenvalues located inside the edge-gap. The existence of the topological bound states in the gap is hence robust with respect to variations in the corner shape. The lattices considered for the diagonalization present an edge with a length of L ≈ 321a and L ≈ 347a, respectively. The slight difference in the energy of the degenerate quadruplet located inside the edge-gap will be discussed later in the article.   [72]). In this case, the lattices considered for the diagonalization present an edge of length of L ≈ 280a and L ≈ 306a respectively. Just as before, these two lattices mainly differ for the disposition of the sites near the vertices, as can be understood by comparing Figure 2c,d. The results obtained by considering a lattice such as the one in Figure 2c or the one in Figure 2d are quite different. In the first case, as can be seen in Figure 3c, the quadruplet of degenerate eigenvalues is located right at the middle of the higher energy gap; on the other hand, in the latter case ( Figure 3d) the quadruplet is located close to the bottom of the lower energy gap. This result demonstrates that, given a gap, the existence of topological bound states crucially depends on the shape of the vertices.
Before proceeding, it is worth discussing the probability density related to the bound states. This is pictured in Figure 4, for the first in-gap state in the set up of Figure 2a. All other eigenstates behave in a qualitatively analogous way. As expected, the bound states are located at the vertices of the structure. Moreover, their localization length is related to the difference in energy separating them to the lowest closest continuum energy states.
From the analysis just reported, we learn that the existence of the eigenvalues associated to the corner modes in a certain gap is strongly dependent on the disposition of the lattice sites around the vertices in the armchair-zigzag case, while it is not in the bearded-armchair one.

Energy
By looking at the spectra in Figure 3, it is evident how the energy of the corner modes is strongly influenced by the disposition of the sites at the vertices of the system. In order to be a little more quantitative, we can try to compare the actual results reported above with the predictions of the low energy theory developed in Ref. [71]. If we call m the half-width of the gap and we refer with E mg to the midgap energy and with E cm to the energy of the corner modes, we can compute the normalized energy of each quadruplet with respect to the corresponding midgap as: The values of such normalized energy for the four spectra in Figure 3 are reported in Table 1. Table 1. Values of the normalized energy ε cm , defined in Equation (4), for the spectra in Figure 3.
Recall that the prediction of the low energy theory for the energy of the corner modes was from Equation (3), ε th cm = cos pπ 4 , with p an integer. Comparing with the values in Table 1, we see that the cases (a) and (c) substantially meet the theoretical value for p = 3 [71] and p = 2 respectively. The other two configurations instead, lead to values of the normalized energy which are not compatible with the predictions of the low energy theory. With reference to Figure 2, we can observe that the cases (b) and (d) correspond to samples with lattice sites located just at the vertices. This leads us to the conclusion that the presence of lattice sites just at the corner may act as a sort of on site chemical potential term, causing a shift of the corner mode energies.
Indeed, this statement is corroborated by the results achieved considering the same lattice of Figure 3b, but with the lattice sites located right at the vertices removed, so that the resulting corners are as in Figure 5a. By diagonalizing the Hamiltonian in Equation (1) for this specific lattice configuration, one finds the low energy spectrum reported in Figure 5b. Its comparison with the previous one in Figure 3b shows that the removal of the corner sites modifies the energy of the corner states, which is now visibly higher.  (1) for a square lattice with an edge of length L ≈ 347a, with the corner sites removed, as depicted in Panel (a). For the diagonalization we set λ ⊥ = 2 and λ H = 0.3t as before and all other parameters as specified in the main text.

Degeneracy
In order to assess the degree of protection of the degeneracy of the bound states, we consider a square sample with bearded-armchair edges, with an edge of length L ≈ 347a, and we remove the lattice sites located at only one of its corners. In other words, the resulting lattice is just like the ones considered in order to obtain the spectra in Figures 3b and 5b, though with three of the corners as in Figure 2b and with the remaining one as in Figure 5a. Strictly speaking, the system considered would now have lost the original exact C 4 symmetry, which is broken just because of the different disposition of the sites at one of the vertices. On the other hand, the C 4 symmetry is still present concerning the bulk and edges of the system. This should be enough to grant the validity of the low energy theory [71], at least for what it concerns the existence of the bound states. The spectrum resulting from the numerical diagonalization is shown in Figure 6.   (1) for a square lattices with an edge of length L ≈ 347a, with the sites located at just one of the four corners removed, as in Figure 5a. For the diagonalization we set λ ⊥ = 2 and λ H = 0.3t as before and all other parameters as specified in the main text. In Panel (b) a zoom of the spectrum around the gap.
One can clearly see that the original quadruplet of Figure 3b splits into a triplet plus a singlet. By carefully examining the precise energies of the eigenvalues, one finds that the triplet is at the energy of the quadruplet in Figure 3b (let us call it ε), while the singlet is at the energy of the quadruplet in Figure 5b (let us call it ε ). Moreover, one can also plot the probability density of the eigenstates associated with the in-gap eigenvalues (not shown here), finding that corner-modes corresponding to the three degenerate eigenvalues with energy ε are localized on the three "untouched" vertices (Figure 2b), while the one with energy ε is localized on the cropped vertex (Figure 5a). This analysis drives us to two main conclusions: first, the degeneracy of the corner modes is protected just as long as the C 4 symmetry is exact on the whole system (corners included). If we geometrically perturb the system at any vertex causing a local C 4 symmetry breaking, we destroy the degeneracy between the corner modes, though they still survive at different energies inside the gap-if originally present. Second, as we already noted in the previous subsections, the energy of the corner modes inside the edge gap crucially depends on the disposition of the lattice sites at the corners.

Conclusions
Our article deals with the quasicrystalline higher order topological insulator predicted in Equation (1). We address three questions related to the properties of the topologically protected zero-dimensional bound states characterizing the model. Such questions are related to the robustness with respect to changes in the structure of the vertices and are the following ones: Do the topological bound states always exist in a given gap? Is their energy fixed? Is their degeneracy protected? We show that the existence of the bound states in a given energy gap is independent of the structure of the vertices if the edges are of the armchair-bearded type, while it is not in the armchair-zigzag case. On the other hand, the energy of the bound states crucially depends on the structure of the vertices in all cases, and the degeneracy is consequently not protected. We hence conclude that if the system is intended to be used as a higher order topological insulator, a precise control of the lattice sites at the vertices is needed.