Implementation and Characterization of the Dice Lattice in the Electron Quantum Simulator

Abstract

Materials featuring touching points, localized states, and flat bands are of great interest in condensed matter and artificial systems due to their implications in topology, quantum geometry, superconductivity, and interactions. In this theoretical study, the experimental realization of the dice lattice with adjustable parameters is proposed by arranging carbon monoxide molecules on a two‐dimensional (2D) electron system at a (111) copper surface. First, a theoretical framework is developed to obtain the spectral properties within a nearly free electron approximation and then compare them with tight‐binding calculations. This investigation reveals that the high mobility of Shockley state electrons enables an accurate theoretical description of the artificial lattice using a next‐nearest‐neighbor tight‐binding model, resulting in the emergence of a touching point, a quasi‐flat band, and localized lattice site behavior in the local density of states. Additionally, theoretical results for a long‐wavelength low‐energy model that accounts for next‐nearest‐neighbor hopping terms are presented. Furthermore, the model's behavior under an external magnetic field is theoretically examined by employing Peierl's substitution, a commonly used technique in theoretical physics to incorporate magnetic fields into lattice models. The theoretical findings suggest that, owing to the exceptional electron mobility, the highly degenerate eigenenergy associated with the Aharonov‐Bohm caging mechanism may not manifest in the proposed experiment.

Introduction

In the past few years, there has been a growing interest in studying lattice systems characterized by spectra containing one or more quasi-flat bands.^[1] This interest stems from the discovery of superconductivity and other correlated phases in layered systems.^[2,3] The prevailing understanding suggests that correlated phenomena are observable in lattice systems with quasi-flat bands, primarily due to the careful comparison of energy scales. Specifically, correlation effects become visible as electronic kinetic energy diminishes to the same order of magnitude as, or even below, the electron interaction energy, thus highlighting an intricate interplay of energy scales at play.

However, the interest of the condensed matter community in the properties of lattice systems with flat bands roots back to the past 30 years, when the effects of frustration and destructive interference were investigated in the kagome,^[4] the dice,^[5] and the Lieb lattice.^[6] All three of these lattice systems are characterized by a perfect flat band that is a finite energy for the case of the kagome lattice and at zero energy for the latter two lattices. The bulk wave function of these flat bands can be expressed as a linear combination of the so-called compact localized states (CLSs).^[7–9] These states, corresponding to wave functions with support only on a finite number of lattice sites, can be regarded as a third possible type of basis alongside the Bloch and Wannier bases. However, it's important to note that CLSs do not necessarily constitute a complete and orthogonal set.^[9] The CLSs associated with a flat band exhibit high degeneracy and possess intriguing quantum geometry properties.^[10] Signatures of CLSs have been observed in different systems presenting flat bands, as the photonic diamond chain^[11] and photonic Lieb lattice,^[12] and in the electronic realm in the quasi-1D diamond-necklace chain.^[13]

Photonic and electronic platforms represent ideal quantum simulators of single particle physics.^[14,15] The concept of quantum simulators can be traced back to a suggestion by Richard Feynman.^[16] He proposed using a controllable quantum system to simulate another quantum system with unknown properties. Since then, many quantum simulation platforms have been developed.^[17] In this article, we will mainly consider a specific condensed-matter quantum simulator, also known as the electron quantum simulator.^[15] In this type of platform, the free electrons of a metallic surface of a noble metal are unambiguously constrained via atoms or molecules to move along a potential profile meticulously designed to resemble one of the specific lattices under investigation, the properties of which are yet to be determined. Within this technique, it has been possible to explore the physics of graphene,^[18,19] aperiodic structures as fractals,^[20] and the Pensore tiling,^[21] in addition to many other different lattice, as summarized in refs. [15, 22, 23].

Of particular interest for this manuscript are the proposal and the realization of the Lieb lattice in the electron quantum simulator.^[24–27] The reason boils down to the spectral properties of the Lieb lattice; this is a face-centered square lattice characterized by three energy bands. In the simplest next-neighbor tight-binding approximation, the energy spectrum of the Lieb lattice is energy symmetric, with one of the three bands completely flat and at zero energy. The three bands share an accidental threefold degeneracy at the M point of the first Brillouin zone (BZ). For this momentum, it is possible to present a long-wavelength low-energy approximation for the Hamiltonian of the Lieb lattice, resulting in an effective pseudo-spin one Hamiltonian.^[28] Additionally, the flat band of the Lieb lattice can be interpreted in terms of CLSs.^[7,9]

Here, we will investigate and propose a possible experimental implementation of another lattice with a flat band at zero energy: the dice lattice.^[5] We present a sketch of the dice lattice in Figure 1a. It is a triangular lattice with a base of three sites in the unit cell; as in the case of the Lieb lattice, the connectivity of the three sites is not homogeneous. Specifically, we have one site, dubbed H with connectivity six, and two inequivalent sites, named R and R' with connectivity three. Similar to the Lieb lattice, the energy spectrum of the dice lattice consists of three bands, two of which are energy symmetric, and the third one is completely flat and at zero energy. The three bands have an accidental degeneracy at the K point of the BZ. In this point, it is possible to describe the long-wavelength low-energy expansion of the system Hamiltonian in terms of an integer pseudo-spin theory.^[29] Some previous work proposed the experimental implementation the dice lattice in the framework of cold atoms in optical lattices.^[29,30]

[IMAGE OMITTED. SEE PDF]

A peculiar property of the dice lattice compared to the Lieb one is the evolution of the energy spectrum under the application of a perpendicular magnetic field. In the absence of disorder, the dice lattice exhibits extreme wavefunction localization, leading to a peculiar energy spectrum characterized by three perfectly flat bands: two with finite and opposite energies and a third at zero energy. This localization phenomenon, known as Aharonov-Bohm caging,^[31–33] arises from the interplay between lattice geometry and the applied magnetic field, contrasting with the ordinary Anderson localization occurring in a disordered system. A few experiments have shown the signature of the Aharonov-Bohm caging in the dice lattice.^[34–37]

In this article, we propose a realization of the dice lattice in the electron quantum simulator. Specifically, we will consider the case of decorating with CO molecules the surface (111) of Cu.^[15,18]. We propose in total three different realizations differing for the amount of CO molecules implemented and the area of the unit cell. We evaluate the spectral properties of these configurations and show how they map into the simplest next-nearest-neighboring tight-binding approximation. The advantage of implementing the dice lattice compared to the Lieb lattice roots is down to the compatibility between the space group of the dice lattice with one of the Cu(111) substrates. It is important to note that the set of parameters we used for our implementation are in agreement with the experimental parameters.^{[18,20,22,25,38]}

The paper is organized as follows. In Section 2.1, we revise the properties of the dice lattice spectrum in the next-nearest-neighboring tight-binding approximation. In Section 2.2, we present our three proposals for the experimental implementation and find the corresponding parameters to map into next-nearest-neighboring tight-binding models. In Section 2.3, we revise the long-wavelength low-energy approximation for the dice lattice Hamiltonian expanded around the K point of the BZ, including terms describing the next-nearest-neighboring hopping terms. Finally, in Sections 2.4, we revise the properties of the dice lattice in a perpendicular magnetic field in the tight-binding approximation. In Sections 3, we summarize the work and give some possible outlook. In Appendix A, we show what is happening when changing the system parameters compared to the standard one used in the experiments.^{[18,20,22,25,38]}

Results

In the following, we describe the general properties of the dice lattice within the simplest tight-binding model without a magnetic field. Successively, we present a possible theoretical implementation with an electron quantum simulator platform based on CO molecules on Cu(111). The numerical implementation is based on the nearly-free electron model, in which a muffin-tin potential models the CO molecules. We compare the results with a next-nearest-neighbor tight-binding model. Finally, we present the results for the spectrum when a perpendicular magnetic field is applied.

Tight-Binding Model

The dice lattice is a regular tessellation of the Euclidean plane with regular π/3 rhombi, resulting in a bipartite structure with hexagonal symmetry. Each unit cell contains three sites; in the following, we name them a hub, indicated with the letter H, and two rims, indicated with R and R', respectively.^[31] We present a sketch of the lattice in Figure 1a. The following translational vectors characterize the unit cell: 1 $$$\begin{align} \bm {a}_1= \frac{a}{2}{\left(3,\sqrt {3}\right)}, \quad \bm {a}_2= \frac{a}{2}{\left(-3,\sqrt {3}\right)} \end{align}$$$where a is the distance between the lattice sites. In the following, we consider a rhombic unit cell for simplicity, while the Wigner-Seitz cell is hexagonal. We describe the dice lattice with a simple s-orbital tight-binding model, including next-neighbor (solid lines in Figure 1a) and next-nearest-neighboring (dashed lines in Figure 1a) hopping terms, t and t′, respectively. Neglecting the next-nearest-neighboring hopping terms t′, the hub sites are sixfold connected, while the rim sites are threefold connected. Thus, including t′ terms will destroy this property, and all the lattice sites will have the same connectivity. This is important since the flat band originates in unequal connectivity among the sites.

The Bloch Hamiltonian that includes both hopping terms reads: 2 $$$\begin{align} \hat{h}(\bm {k})= \def\eqcellsep{&}\begin{pmatrix} 0 && f_{\text{RH}}(\bm {k})\, t && f_{\text{RR'}}(\bm {k})\, t^{\prime } \\[6pt] f_{\text{RH}}^\ast (\bm {k})\, t && -\Delta && f_{\text{HR'}}(\bm {k})\, t\\[6pt] f_{\text{RR'}}^\ast (\bm {k})\, t^{\prime }&& f_{\text{HR'}}^\ast (\bm {k})\, t&& 0 \end{pmatrix} \end{align}$$$expressed in the following base Ψ = (ψ_R, ψ_H, ψ_R')^T. Here, we have introduced the following auxiliaries functions: 3 $$$\begin{align} f_{\text{RH}}(\bm {k})&= 1+\text{e}^{-\text{i}\, \bm {k}\cdot \bm {a}_2} +\text{e}^{-\text{i}\, (\bm {k}\cdot \bm {a}_1+\bm {k}\cdot \bm {a}_2)}\end{align}$$$4 $$$\begin{align} f_{\text{HR'}}(\bm {k})&=1+\text{e}^{-\text{i}\, \bm {k}\cdot \bm {a}_1} +\text{e}^{-\text{i}\, (\bm {k}\cdot \bm {a}_1+\bm {k}\cdot \bm {a}_2)}\end{align}$$$5 $$$\begin{align} f_{\text{RR'}}(\bm {k})&=1+\text{e}^{-\text{i}\, \bm {k}\cdot \bm {a}_1 } +\text{e}^{\text{i}\, \bm {k}\cdot \bm {a}_2} \end{align}$$$where $$\bm {k}$$ is the crystal momentum vector defined within the BZ, see Figure 1b. In Equation (2), we have introduced an energy term Δ describing the difference between the hub and rim on-site energy. It is not possible to find a compact expression for the energy bands; however, for t′ = 0, the energy spectrum is given by 6 $$$\begin{align} \mathcal {E}_{\text{c}}&=0\end{align}$$$7 $$$\begin{align} \mathcal {E}_\alpha (\bm {k})&=\alpha t\sqrt {6 + {\left(\frac{\Delta }{4 t}\right)}^2+4 \Lambda (\bm {k})}-\frac{\Delta }{2} \end{align}$$$where we have introduced the function 8 $$$\begin{align} \begin{aligned} \Lambda (\bm {k})&=\cos [\bm {k}\cdot \bm {a}_1]+\cos [\bm {k}\cdot \bm {a}_2]+\cos [\bm {k}\cdot (\bm {a}_1+\bm {a}_2)] \\ & = 2\cos {\left(\frac{3 k_x a}{2}\right)}\cos {\left(\frac{\sqrt {3}k_y a}{2}\right)}+\cos (\sqrt {3}k_y a)\, \end{aligned} \end{align}$$$

We start presenting the spectral properties of the dice lattice in the simple case of t′ = Δ = 0. In this limit, of the three bands composing the spectrum, two are dispersive $$\mathcal {E}_\pm$$, and one is completely flat at zero energy $$\mathcal {E}_{\text{c}}$$ — Figures 1c. The two dispersive bands are identical to the one of the honeycomb lattice.^[39] The flat band $$\mathcal {E}_{\text{c}}$$ arises because of the bipartite nature of the dice lattice^[8,32] and can be interpreted in terms of CLSs as well.^[40] The three bands are degenerate in the K and K′ points of the BZ. This threefold degeneracy is accidental since, in two dimensions, there are no irreducible representations of order three.^[41,42] Performing a low-energy long-wavelength approximation of the energy spectrum around this K or K′, it is possible to obtain an effective pseudo-spin one Hamiltonian model for the system,^[28–30,43] see Section 2.3 for additional details.

By inspecting Equation (7), we notice that the onsite energy Δ shifts the two dispersive bands and partially removes the degeneracy at the K(K′)-points by opening a gap between a dispersive band and the flat one — Figure 2a. Changing the sign of Δ, the energy gap opens between the flat band and the other dispersive one. Thus, we conclude that this term breaks the threefold degeneracy at the K(K′)-points. Furthermore, a finite value of the next-nearest-neighboring hopping t′ induces a finite dispersion in the flat band — Figure 2b, but it does not open an energy gap at the threefold degeneracy point in K. The energy symmetry is broken when Δ and t′ are finite — see Figure 1d. We note on passing that the effects of t′ are similar to the effects of evaluating the energy spectrum by considering an overlap integral.^[44]

[IMAGE OMITTED. SEE PDF]

In Figure 2b,d, we present the density of states (DOS) corresponding to the spectra depicted in Figure 2a,c, respectively. We define the DOS as9 $$$\begin{equation} \rho (\epsilon)=\sum _n \int _{\text{BZ}} \frac{d\bm {k}}{2 \pi ^2}\, \delta (\epsilon -\epsilon _{n \bm {k}}) \end{equation}$$$where $$\epsilon _{n \bm {k}}$$ is the energy of the band n for the 2D momentum $$\bm {k}$$ within the BZ and we have taken into account the spin degeneracy. Notably, it features characteristic Van Hove singularities, indicative of the saddle points at the M point in the energy spectrum, see Figure 2b, whereas, when t′ ≠ 0, it features a double peak one close to the M point and another one close in energy, see Figure 2d. Additionally, in both cases, it exhibits a distinct peak associated with the flat or quasi-flat band $$\mathcal {E}_{\text{c}}$$.

Finally, we present in Figure 2e the bandwidth evolution of the $$\mathcal {E}_{\text{c}}$$ band as a function of Δ, and t′. We have defined the bandwidth as difference between the maximum and the minumum value of $$\mathcal {E}_{\text{c}}(\bm {k})$$, with $$\bm {k}$$ in the BZ. Here, we can observe that the bandwidth grows linearly with the next-nearest-neighboring hopping t′ up to values of this of the order of t′ ≈ 0.7t, then it starts to change slope, especially for large values of the onsite energy Δ. Additionally, we note that for t′ < t/2, the bandwidth is independent of the onsite energy Δ.

Implementation with CO Molecules on Cu(111)

In this section, we explore diverse implementations of the dice lattice within the electron quantum simulator. These are realized by decorating the Cu(111) surface with CO molecules.^[15,22,23] The Cu(111) surface is characterized by a triangular lattice structure featuring an interatomic distance of a_Cu = 0.256 nm. Specifically, we introduce three distinct configurations for the dice lattice, as illustrated in Figure 3, each distinguished by the number of CO molecules utilized to construct the unit cell. In all proposed configurations, CO molecules play a crucial role in designing the unit cell's lattice sites with varying connectivity and in tuning the hopping probabilities between the rim sites. This tunneling is possible since electrons can tunnel beneath the CO potential barriers, enabling a nonzero next-nearest-neighbor hopping t′.^[22] It is noteworthy that the sizes of the three sites within the lattice are not uniform, leading to different onsite energies for the hub and the rim sites.

[IMAGE OMITTED. SEE PDF]

In the following, we summarize the key steps for obtaining the spectral properties for these three dice lattice implementations with CO molecules. In general, we want to find the spectral properties of the following Hamiltonian 10 $$$\begin{equation} H=\frac{\bm {p}^2}{2m^*}+V_{\text{CO}}(\bm {R}) \end{equation}$$$where m* is the effective mass of the free electrons on Cu(111),^[22] the periodic potential $$V_{\text{CO}}(\bm {R})$$ is obtained by considering the set of all CO molecules. Without loss of generality, we can assume that each CO molecule at the position $$\bm {r}_0$$ is a cylindrical repulsive potential for the Cu(111) surface electrons:^[18] 11 $$$\begin{equation} V(\bm {r}-\bm {r}_0) = V_0 {\begin{cases} 1 & |\bm {r}- \bm {r}_0| \leqslant \frac{D}{2} \\[6pt] 0 & |\bm {r}- \bm {r}_0|&gt; \frac{D}{2} \end{cases}} \end{equation}$$$The parameters of this potential are chosen in agreement with the experiments as V₀ = 0.9 eV and a cylinder diameter of D = 0.6 nm.^[18,22,25] In Appendix A, we show what happens by changing the shape of the potential describing the CO molecules and the behavior by changing the system parameters away from the optimal ones. The complete set of CO molecules can be described by a periodic potential $$V_{\text{CO}}(\bm {R})$$ that can be expanded in terms of Fourier components as 12 $$$\begin{equation} V_{\text{CO}}(\bm {R})=\sum _{\bm {G}} V_{\bm {G}}\, \text{e}^{\text{i}\bm {G}\cdot \bm {R}} \end{equation}$$$here $$\bm {G}$$ are vectors of the reciprocal lattice. Because of the shape of each cylindrical potential in Equation (11), the Fourier components $$V_{\bm {G}}$$ are 13 $$$\begin{align} V_{\bm {G}}& = \int _{\text{UC}} d\bm {r}\ V(\bm {r}-\bm {r}_0)\text{e}^{-\text{i} \bm {G}\cdot \bm {r}} \nonumber \\ & = (2\pi r_0 V_0) \frac{\mathcal {J}_1(|\bm {G}|D)}{|\bm {G}|}\, \end{align}$$$where $$\mathcal {J}_1(x)$$ are the Bessel functions of the first order, and the integral is performed over the unit cell (UC). Summing up over all the Fourier components $$V_{\bm {G}}$$, we can rewrite Equation (12) as: 14 $$$\begin{align} V_{\text{CO}}(\bm {R}) & = \frac{2\pi r_0 V_0}{A} \sum _{\bm {G}}\frac{\mathcal {J}_1(|\bm {G}|r_0)}{|\bm {G}|} \text{e}^{\text{i} \bm {G}\cdot \bm {R}} \nonumber\\ & = \frac{2\pi r_0 V_0}{A}\!\!\!\! \sum _{\lbrace \ell,m\rbrace \in \mathbb {Z}}\!\!\frac{\mathcal {J}_1(|\ell \bm {b}_1+m\bm {b}_2|D)}{|\ell \bm {b}_1+m\bm {b}_2|} \text{e}^{\text{i} (\ell \bm {b}_1+m\bm {b}_2)\cdot \bm {R}} \end{align}$$$

At this point, assuming that the system wave function is also periodic,^[45] we can transform the Schrödinger equation in a set of coupled linear equations 15 $$$\begin{align} {\left[\frac{\hbar ^2 {\left|\bm {q}-\bm {G}\right|}^2}{2m}-\mathcal {E}_n(\bm {q})\right]}&c^n_{\bm {q-G}} \\ &\hspace{-28.45274pt}+ \sum _{\bm {G^{\prime }}} V_{\bm {G}^{\prime }-\bm {G}} \, c^n_{\bm {q-G}}=0 \nonumber\end{align}$$$where $$\bm {q}$$ is the crystalline momentum defined in the BZ, and n the band index that can vary between one and the total number of reciprocal lattice vectors $$\bm {G}$$ in Equation (12). Solving this set of n linear equations results in the band structure $$\mathcal {E}_n(\bm {q})$$ for the periodic potential in Equation (12) and the Bloch wave function given by 16 $$$\begin{equation} \psi _{n\bm {q}}(\bm {r})=\text{e}^{\text{i}\bm {q}\cdot \bm {r}} \sum _{\bm {G}} c^n_{\bm {q}-\bm {G}} \,\text{e}^{\text{i}\bm {G}\cdot \bm {r}} \end{equation}$$$In this expression, the term denoted by the summation defines the periodic part $$u_{n\bm {q}}(\bm {r})$$ of the Bloch wave function $$\psi _{n\bm {q}}(\bm {r})$$. Finally, we can use the Bloch wave function in Equation (16) to calculate the local density of states (LDOS):^[46] 17 $$$\begin{equation} \rho (\bm {r},\epsilon)=\sum _n \int _{\text{BZ}} \frac{d\bm {k}}{2\, \pi ^2}\,\vert \psi _{n \bm {k}}(\bm {r})\vert ^2\, \delta (\epsilon -\epsilon _{n \bm {k}}) \end{equation}$$$which can be measured in the STM experiments.^[25] We note on passing that integrating $$\rho (\bm {r},\epsilon)$$ over the unit cell results in the total DOS as in Equation (9).

As stated before, we consider the three configurations of the dice lattice we presented in Figure 3. These differ for the total number of CO molecules, 9 for the case of CI Figure 3a, 24 for the case of CII Figure 3b, and finally, 15 for the case of CIII Figure 3c. A summary of the properties of the three configurations is given in Table 1. We present in Figure 4a,b,c the band structure for the three possible realizations of the dice lattice, CI, CII, and CIII, respectively. Together with the band structure obtained by solving Equation (15), we present the optimal fitting with the tight-binding model from Equation (2) at the K-point; it has to be noted that we include a homogeneous shift of the spectrum ϵ.

Table 1 Summary of the properties of the three configurations in Figure 3.

	ϵ/t	Δ/t	t′/t	#CO unit cell	□ (nm2)
CI	−1.18182	1.71499	−0.87715	9	8
CII	−5.67708	0.128125	−0.428646	24	13
CIII	−5.51934	0.0431492	−0.40663	15	18

[IMAGE OMITTED. SEE PDF]

In Figure 4d,e, f, we present the LDOS $$\rho (\bm {r},E)$$ at the center of the rim and hub sites for the three configurations, this quantity would be accessible via STM measurements. These results help us to understand the role played by the next-nearest-neighboring model (t′ ≠ 0). As introduced before, the perfectly flat band case t′ = 0 can be characterized by a linear combination of CLSs around the hub sites.^[5,8,40,47] Surprisingly, for t′ ≠ 0, we can observe that the LDOS is still localized around the rim sites for energies corresponding to the $$\mathcal {E}_{\text{c}}$$ band — Figure 2b,d. The survival of imperfect CLSs is more evident from the plots of the LDOS in space Figure 4g,h,i, these density plots for $$\rho (\bm {r},E)$$ are obtained for an energy corresponding crossing between the upper and the central band at the K point. Here, it is possible to recognize that the LDOS still has a shape resembling the CLSs around the hub sites, with values that are much larger on the rim sites compared to the hub sites.^[8]

Low-Energy Theory

In this section, we will present a low-energy, long-wavelength approximation for the Hamiltonian (2) around the K/K′ points. In the limit of t′ = 0, the approximated Hamiltonian reads 18 $$$\begin{equation} \mathcal {H}_0 = v_{\text{F}}\, \bm {\Sigma }\cdot \bm {p} \end{equation}$$$where, v_F = 3ta/2 is the Fermi velocity, and $$\bm {p}=-\text{i}\hbar (\partial _x,\partial _y,0)$$ is the momentum operator in the lattice xy-plane. We can observe that Equation (18) bears a remarkable resemblance to the Hamiltonian governing the behavior of electrons in graphene. In Equation (18), there is a significant difference as compared to the graphene one.^[39] Here, the pseudo-spin vector Σ = (Σ_x, Σ_y, Σ_z) has a total spin S = 1, which reflects that the dice lattice has three inequivalent sites per unit cell as compared to the two sites in graphene. These 3 × 3 matrices can be expressed as a linear combination of Gell-Mann matrices^[48] λ_i and read as follow: 19 $$$\begin{align} \Sigma _x&=\frac{1}{\sqrt {2}}(\lambda _1+\lambda _6) = \frac{1}{\sqrt {2}}\def\eqcellsep{&}\begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 0 \end{pmatrix}\end{align}$$$20 $$$\begin{align} \Sigma _y&=\frac{1}{\sqrt {2}}(\lambda _2+\lambda _7) =\frac{1}{\sqrt {2}}\def\eqcellsep{&}\begin{pmatrix} 0 & -\text{i}& 0 \\ \text{i}& 0 & -\text{i}\\ 0 & \text{i}& 0 \end{pmatrix}\!\! \end{align}$$$21 $$$\begin{align} \Sigma _z&=(\lambda _3+\sqrt {3}\lambda _8)=\def\eqcellsep{&}\begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -1 \end{pmatrix} \end{align}$$$which satisfy angular momentum commutation relations [S_i, S_j] = iϵ_ijkS_k, with ϵ_ijk the Levi-Civita tensor.^[29,30,49] We can introduce a chiral symmetry operator defined as 22 $$$\begin{align} \Gamma =\def\eqcellsep{&}\begin{pmatrix} 1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \end{align}$$$that is unitary $$\Gamma ^\dag \Gamma =\Gamma ^2=\mathbb {I}_3$$ where $$\mathbb {I}_3$$ is the identity matrix of dimension 3, and anticommutes with the Hamiltonian (18): $$\lbrace \Gamma,\mathcal {H}_0\rbrace =0$$. This chiral symmetry operator governs the energy symmetry of the system, see Figure 1c.^[38,50–52]

We can complement the matrices (19) with two additional Gell-Mann matrices, λ₄ and λ₅, so to describe the corresponding low-energy next-nearest-neighboring hopping term t′. This Hamiltonian term reads: 23 $$$\begin{align} \mathcal {H}_\Xi = \Xi (\lambda _4 p_x+\lambda _5 p_y) \end{align}$$$where the coupling constant is Ξ = 3t′a/2. Importantly, this additional term breaks the chiral symmetry of the system: $$\lbrace \Gamma,\mathcal {H}_\Xi \rbrace \ne 0$$ and thus the energy symmetry, see Figure 2c. For Ξ = 0, we can write a simple expression for the energy spectrum: 24 $$$\begin{align} E^{\text{lwa}}_\alpha =\alpha v_{\text{F}} \sqrt {k_x^2+k_y^2} \end{align}$$$with α ∈ {0, ±}, whereas, for Ξ ≠ 0, it is not possible to present a compact expression for the energy spectrum. In Figure 5a,b we present the spectrum associated to Equations (18) and (23) for the dice lattice when Ξ is zero or finite, respectively. In Figure 5a, we can identify the threefold degeneracy with the flat band in the (0,0) point corresponding to the K point where we have performed the expansion of the Hamiltonian (2). In Figures 5b, for finite Ξ, we observe that the linear dispersion around the (0,0) point is maintained but the flat band acquires a dispersion, and there is breaking of the energy symmetry.

[IMAGE OMITTED. SEE PDF]

Magnetic Field Properties: Tight-Binding Approach

In the following, we show how to obtain the electronic energy spectrum in the presence of an external magnetic field in the tight-binding approximation. As customary, we perform the Peierls substitution in the hopping amplitude:^[53] 25 $$$\begin{align} t_{(\bm {R},\ell)\rightarrow (\bm {R}^{\prime },\ell ^{\prime })} \rightarrow t_{(\bm {R},\ell)\rightarrow (\bm {R}^{\prime },\ell ^{\prime })}\, \text{e}^{-\text{i}\, \frac{\hbar }{e}\, \Phi _{(\bm {R},\ell)\rightarrow (\bm {R}^{\prime },\ell ^{\prime })}} \end{align}$$$with 26 $$$\begin{align} \Phi _{(\bm {R},\ell)\rightarrow (\bm {R}^{\prime },\ell ^{\prime })} = \int _{\bm {R}_\ell }^{\bm {R}^{\prime }_{\ell ^{\prime }}} d\bm {r}\cdot \bm {A}(\bm {r}) \end{align}$$$where ℓ, ℓ′ = {R', H, R} and $$\bm {R}_{\ell } = n_1\, \bm {a}_1 +n_2\, \bm {a}_2 + \bm {\delta }_\ell$$ is the position of the ℓ-th site in the unit cell, with $$\bm {\delta }_\ell$$ the displacement within the unit cell.

We choose the Landau gauge $$\bm {A}=(0,B\, x,0)$$. Then, considering, for example, the $$\hat{y}$$-axis parallel to $$\bm {a}_2$$, the lattice periodicity is partially conserved along the y-direction. Within this choice of gauge, it is still possible to consider a change of basis 27 $$$\begin{align} \mathinner {|{n_1,n_2,l}\rangle } = \frac{1}{\sqrt {N_2}}\sum _{k_2} \text{e}^{-\text{i}\, k_2\, n_2}\, \mathinner {|{n_1,k_2,l}\rangle } \end{align}$$$in which the tight-binding Hamiltonian reads 28 $$$\begin{multline} \hat{H} = \sum _{k_2}\sum _{n_1}^{\text{PBC}} \Bigl [ \bm {c}_{n_1,k_2}^\dagger {\left(\mathcal {M} + \text{e}^{\text{i}\, k_2}\, \mathcal {P}_2 \right)} \bm {c}_{n_1,k_2} \\ +\bm {c}_{n_1,k_2}^\dagger {\left(\mathcal {P}_1 + \text{e}^{-\text{i}\, k_2}\, \mathcal {P}_3 \right)} \bm {c}_{n_1+1,k_2} + \text{h.c.}\Bigr]\end{multline}$$$where $$\bm {c}_{n_1,k_2}^\dagger =(\hat{c}^\dagger _{n_1,k_2,\text{R}},\, \hat{c}^\dagger _{n_1,k_2,\text{H}},\, \hat{c}^\dagger _{n_1,k_2,\text{R'}})$$ are the creation operators, the superscript in the n₁-summation refers to the periodic boundary conditions and 29 $$$\begin{align} \mathcal {M} &= \frac{t}{4} \def\eqcellsep{&}\begin{pmatrix} 0 & \text{e}^{-\frac{1}{2}\, \text{i}\, \pi f (6 \, n_1-1)} & 0\\ \text{e}^{\frac{1}{2}\, \text{i}\, \pi \, f\, (6\, n_1-1)} & -\frac{\Delta }{t} & \text{e}^{-\frac{1}{2}\, \text{i}\, \pi \, f\, (6\, n_1+1)}\\ 0 & \text{e}^{\frac{1}{2}\, \text{i}\, \pi \, f\, (6\, n_1+1)} & 0 \end{pmatrix} \end{align}$$$30 $$$\begin{align} \mathcal {P}_1 & = \def\eqcellsep{&}\begin{pmatrix} 0 & 0 & 0 \\ t & 0 & 0 \\ t^{\prime }\, \text{e}^{\frac{3}{2}\, \text{i}\, \pi \, f\, (2\, n_1+1)} & t & 0 \end{pmatrix}\end{align}$$$31 $$$\begin{align} \mathcal {P}_2 &= \def\eqcellsep{&}\begin{pmatrix} 0 & t\, \text{e}^{\frac{1}{2}\, \text{i}\, \pi \, f\, (6\, n_1-1)} & t^{\prime } \\ 0 & 0 & t\, \text{e}^{\frac{1}{2}\, \text{i}\, \pi \, f\, (6\, n_1+1)} \\ 0 & 0 & 0 \end{pmatrix}\end{align}$$$32 $$$\begin{align} \mathcal {P}_3 &= \def\eqcellsep{&}\begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ t^{\prime }\, \text{e}^{-\frac{3}{2}\, \text{i}\, \pi \, f\, (2\, n_1+1)} & 0 & 0 \end{pmatrix} \end{align}$$$contain the magnetic hopping terms. We have omitted, for simplicity in Equation (29), a constant shift ϵ corresponding to the effective shift of the tight-binding model fitted over the three experimental proposals.

Typically, for the dice lattice, the magnetic flux is measured through an elementary rhombus (one-third of the unit cell surface) and normalized with respect to the ϕ₀ = h/e magnetic flux quantum, that is $$f= B\, \sqrt {3}\, a^2/(2\, \phi _0)$$. The calculation of the spectrum as a function of f gives rise to the Hofstadter's butterflies^[54] presented in Figures 6 and 7. In particular, Figure 6a represents the ideal nearest neighbor t = 1, t′ = 0 case,^[31] whereas in Figure 6b, we present the case of the triangular lattice corresponding to the t = 1, t′ = −t case.^[55] Interestingly, for the case of the dice lattice at f = 1/2, all the states collapse in only three energy eigenvalues. This is a consequence of the Aharonov-Bohm caging mechanism,^[31] for which the number of sites visited by an initially localized wave packet is limited due to destructive interferences, resulting in an extreme localization phenomenon.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

As we have seen in Section 2.2, due to the high mobility of the electrons of the Shockley state, we have finite values of the t′ parameter; consequently, the response in a magnetic field of our experimental proposal would substantially differ from the idea case in Figure 6a. Considering the fitting parameters in Table 1, we obtain the Hofstadter butterflies represented in Figure 7. We can clearly see how these fractal structures are an interpolation^[56] of the pure Hofstadter butterfly for the case of the dice lattice in Figure 6a and the one for the triangular lattice in Figure 6b.

Conclusion and Outlook

In this study, we explored the realization of a dice lattice on a Cu(111) surface decorated with CO molecules. Our experimental proposal demonstrates that the key spectral features of the dice lattice are maintained in some of the proposed implementations. We have verified this by fitting the electronic spectral properties using a next-nearest-neighbor tight-binding model. Some of the proposed implementations reveal a touching point and a quasi-flat band. We note on passing that one of the largest experiments performed with the electron quantum simulator for the realization of the Pensore tiling^[21] included around 460 CO molecules. Extrapolating on this number would mean we could achieve around 20 unit cells for the dice lattice in the CII configuration. This number of unit cells is sufficient to already observe bulk properties in finite-size systems.

The finite bandwidth of the central quasi-flat band in our model is analogous to the Lieb lattice case, theoretically predicted^[24] and experimentally implemented.^[25] However, the actual value of the bandwidth of our proposal is of the order lower than 0.1 eV in the two gapless configurations compared to the value of almost 0.2 eV of the Lieb lattice. This proves that the experimental realization of the dice lattice can result in a quasi-flat band with a smaller bandwidth. By studying the LDOS, we have shown a clear signature of the quasi-flat band in terms of a site-specific localization resembling CLSs; this spatial arrangement can be observed experimentally with a scanning tunneling microscope technique.^[25]

As a final perspective for our results, we envision two important developments: as for the case of the honeycomb lattice, an opportune modulation of the lattice periodicity can result in a homogeneous shift of the Fermi energy.^[18] This property can use used to create an effective pn-junction where the super-Klein tunneling^[43] could be observed. Second, to enhance the possibility of observing the effects of electron correlation in the flat band, we plan to investigate in future studies the implementation of the dice lattice in semiconducting electron simulator platforms based on InAs and InSb decorated with In or Cs adatoms.^[57,58] Finally, regarding the external magnetic field, an interesting outlook would be to study the strain-induced pseudomagnetic field in the case of next-nearest-neighbor term^[59] since intensities of the order of 60 hundred Tesla have been achieved for the case the honeycomb lattice.^[18]

Appendix

Appendix A - Model and Parameters

Here, we discuss an additional shape for the muffin-thin potential and the dependence of the parameters for the case of the cylindrical one.

It turns out that the muffin-tin model with cylindrical potentials approximates very well the experimental results by choosing the parameters a potential barrier V₀ = 0.9 eV and a diameter D = 0.6 nm.^[18,22,25] An alternative description of the CO molecules is a 2D isotropic Gaussian potential: A1 $$$\begin{gather} V_{\text{Gauss}}(\bm {r}-\bm {r}_0) = V_{0}^{(G)}\, \text{e}^{-{\left(\frac{\bm {r}-\bm {r}_0}{\sqrt {2}\, \alpha }\right)}^2} \end{gather}$$$with $$V_0^{(G)}=1.8\, \text{eV}$$ and α = 0.15 nm, such that A2 $$$\begin{align} \int _{\text{UC}} d\bm {r}\, V_{\text{Gauss}}^n(\bm {r}) \sim \int _{{\bf R}^2} d\bm {r}\, V_{\text{Gauss}}^n(\bm {r}) = \int _{\text{UC}} d\bm {r}\, V^n(\bm {r}) \end{align}$$$where n = 1, 2 and V is cylindrical potential. In particular, we can assume that the Gaussian vanishes fast enough within the unit cell and we require that the first and second momenta of the two potentials coincide.

The Fourier components $$V_{\bm {G}}$$ are in the Gaussian case A3 $$$\begin{align} V_{\bm {G}}& = \int _{\text{UC}} d\bm {r}\ V_{\text{Gauss}}(\bm {r}-\bm {r}_0)\text{e}^{-\text{i} \bm {G}\cdot \bm {r}} \nonumber \\ & = (2\pi \alpha ^2 V_0) \text{e}^{-\vert \bm {G}\vert ^2 \alpha ^2/2}\, \end{align}$$$and we found that the kind of potential does not affect significantly the nearly free electron model (in Figure A1, we report the case with 24 CO molecules).

[IMAGE OMITTED. SEE PDF]

Finally, by increasing the muffin-tin cylindrical potential strength V₀, we may expect higher values of the on-site energies and higher ratio t/t′ between the nearest and next-nearest-neighbor hopping terms, which is a measure of the central band flatness (note that V₀ is tunable only by changing the CO molecule with one with different repulsive potential^[24]). This is shown in Figure A2 for the 24 CO molecules case, where the potential has increased by about 10%. This results in a lower gap Δ and in an increase of the t/t′ ratio by about 6%.

[IMAGE OMITTED. SEE PDF]

Acknowledgements

The authors acknowledged the useful discussion with Miguel Angel Cazalilla, Miguel Angel Jimenez Herrera, Rian Ligthart, and Ingmar Swart. D.B. acknowledged the support from the Spanish MICINN-AEI through Project No. PID2020-120614GB-I00 (ENACT), the Transnational Common Laboratory Quantum − ChemPhys, the Department of Education of the Basque Government through the project PIBA_2023_1_0007 (STRAINER), and the financial support received from the IKUR Strategy under the collaboration agreement between the Ikerbasque Foundation and DIPC on behalf of the Department of Education of the Basque Government and the Gipuzkoa Provincial Council within the QUAN-000021-01 project.

Conflict of Interest

The authors declare no conflict of interest.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

References

P. Törmä, S. Peotta, B. A. Bernevig, Nat. Rev. Phys. 2022, 4, 528.

E. Y. Andrei, A. H. MacDonald, Nat. Mater. 2020, 19, 1265.

D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. N. Basov, A. N. Pasupathy, A. Rubio, Nat. Phys. 2021, 17, 155.

R. J. Baxter, J. Math. Phys. 1970, 11, 784.

B. Sutherland, Phys. Rev. B 1986, 34, 5208.

E. H. Lieb, Phys. Rev. Lett. 1989, 62, 1201.

H. Aoki, M. Ando, H. Matsumura, Phys. Rev. B 1996, 54, [eLocator: R17296].

D. L. Bergman, C. Wu, L. Balents, Phys. Rev. B 2008, 78, [eLocator: 125104].

J.‐W. Rhim, B.‐J. Yang, Phy. Rev. B 2019, 99, [eLocator: 045107].

A. Graf, F. Piéchon, Phys. Rev. B 2021, 104, [eLocator: 085114].

S. Mukherjee, M. D. Liberto, P. Öhberg, R. R. Thomson, N. Goldman, Phys. Rev. Lett. 2018, 121, [eLocator: 075502].

S. Mukherjee, A. Spracklen, D. Choudhury, N. Goldman, P. Öhberg, E. Andersson, R. R. Thomson, Phys. Rev. Lett. 2015, 114, [eLocator: 245504].

S. N. Kempkes, P. Capiod, S. Ismaili, J. Mulkens, L. Eek, I. Swart, C. M. Smith, Quantum Frontiers 2023, 2, 1.

A. Aspuru‐Guzik, P. Walther, Nat. Phys. 2012, 8, 285.

A. A. Khajetoorians, D. Wegner, A. F. Otte, I. Swart, Nat. Rev. Phys. 2019, 1, 703.

R. P. Feynman, Int. J. Theor. Phys. 1982, 21, 467.

E. Altman, et al.PRX Quantum 2021, 2, [eLocator: 017003].

K. K. Gomes, W. Mar, W. Ko, F. Guinea, H. C. Manoharan, Nature 2012, 483, 306.

S. Wang, L. Z. Tan, W. Wang, S. G. Louie, N. Lin, Phys. Rev. Lett. 2014, 113, [eLocator: 196803].

S. N. Kempkes, M. R. Slot, S. E. Freeney, S. J. M. Zevenhuizen, D. Vanmaekelbergh, I. Swart, C. M. Smith, Nat. Phys. 2019, 15, 127.

L. C. Collins, T. G. Witte, R. Silverman, D. B. Green, K. K. Gomes, Nat. Commun. 2017, 8, [eLocator: 15961].

S. E. Freeney, M. R. Slot, T. S. Gardenier, I. Swart, D. Vanmaekelbergh, ACS Nanoscience Au 2022, 2, 198.

I. Piquero‐Zulaica, J. Lobo‐Checa, Z. M. A. El‐Fattah, J. E. Ortega, F. Klappenberger, W. Auwärter, J. V. Barth, Rev. Mod. Phys. 2022, 94, [eLocator: 045008].

W.‐X. Qiu, S. Li, J.‐H. Gao, Y. Zhou, F.‐C. Zhang, Phys. Rev. B 2016, 94, [eLocator: 241409].

M. R. Slot, T. S. Gardenier, P. H. Jacobse, G. C. Van Miert, S. N. Kempkes, S. J. Zevenhuizen, C. M. Smith, D. Vanmaekelbergh, I. Swart, Nat. Phys. 2017, 13, 672.

R. Drost, T. Ojanen, A. Harju, P. Liljeroth, Nat. Phys. 2017, 13, 668.

D. Bercioux, S. Otte, Nat. Phys. 2017, 13, 628.

N. Goldman, D. F. Urban, D. Bercioux, Phys. Rev. A 2011, 83, [eLocator: 063601].

D. Bercioux, D. F. Urban, H. Grabert, W. Häusler, Phys. Rev. A 2009, 80, [eLocator: 063603].

D. Bercioux, N. Goldman, D. F. Urban, Phys. Rev. A 2011, 83, [eLocator: 023609].

J. Vidal, R. Mosseri, B. Douçot, Phys. Rev. Lett. 1998, 81, 5888.

J. Vidal, B. Douçot, R. Mosseri, P. Butaud, Phys. Rev. Lett. 2000, 85, 3906.

J. Vidal, P. Butaud, B. Douçot, R. Mosseri, Phys. Rev. B 2001, 64, [eLocator: 155306].

C. C. Abilio, P. Butaud, T. Fournier, B. Pannetier, J. Vidal, S. Tedesco, B. Dalzotto, Phys. Rev. Lett. 1999, 83, 5102.

C. C. Naud, G. G. Faini, D. D. Mailly, Phys. Rev. Lett. 2001, 86, 5104.

H. Li, Z. Dong, S. Longhi, Q. Liang, D. Xie, B. Yan, Phys. Rev. Lett. 2022, 129, [eLocator: 220403].

J. G. Martinez, C. S. Chiu, B. M. Smitham, A. A. Houck, Sci. Adv. 2023, 9, [DOI: https://dx.doi.org/10.1126/sciadv.adj7195].

S. N. Kempkes, M. R. Slot, J. J. van den Broeke, P. Capiod, W. A. Benalcazar, D. Vanmaekelbergh, D. Bercioux, I. Swart, C. Morais Smith, Nat. Mater. 2019, 18, 1292.

A. H. C. Neto, A. H. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, A. K. Geim, Rev. Mod. Phys. 2009, 81, 109.

S. M. Zhang, L. Jin, Phys. Rev. B 2020, 102, [eLocator: 054301].

B. Bradlyn, J. Cano, Z. Wang, M. G. Vergniory, C. Felser, R. J. Cava, B. A. Bernevig, Science 2016, 353, [eLocator: aaf5037].

B. Bradlyn, L. Elcoro, J. Cano, M. G. Vergniory, Z. Wang, C. Felser, M. I. Aroyo, B. A. Bernevig, Nature 2017, 547, 298.

D. F. Urban, D. Bercioux, M. Wimmer, W. Häusler, Phys. Rev. B 2011, 84, [eLocator: 115136].

M. O. Goerbig, Rev. Mod. Phys. 2011, 83, 1193.

N. W. Ashcroft, N. D. Mermin, Solid state physics, Holt, Rinehart and Winston, New York, NY 1976.

H. J. Monkhorst, J. D. Pack, Phys. Rev. B 1976, 13, 5188.

A. Graf, F. Piéchon, Phys. Rev. B 2021, 104, [eLocator: 195128].

M. Gell‐Mann, Phys. Rev. 1962, 125, 1067.

D. Bercioux, O. Dutta, E. Rico, Ann. Phys. (Berl.) 2017, 529, [eLocator: 1600262].

Here the chiral symmetry is simply anticommuting with the system Hamiltonian since the flat band is at zero energy. In other lattice systems with three sites in the unit cell, such as the kagome lattice, the flat band is at finite energy. Thus they require the introduction of a generalized chiral symmetry.^{[38, 51, 52]}

X. Ni, M. Weiner, A. Alu, A. B. Khanikaev, Nat. Mater. 2018, 61, 1.

M. A. J. Herrera, S. N. Kempkes, M. B. de Paz, A. García‐Etxarri, I. Swart, C. M. Smith, D. Bercioux, Phys. Rev. B 2022, 105, [eLocator: 085411].

R. Peierls, Z. Phys. 1933, 80, 763.

D. R. Hofstadter, Phys. Rev. B 1976, 14, 2239.

L. Du, Q. Chen, A. D. Barr, A. R. Barr, G. A. Fiete, Phys. Rev. B 2018, 98, [eLocator: 245145].

M. F. Barnsley, Constr. approx. 1986, 2, 303.

V. D. Pham, Y. Pan, S. C. Erwin, F. v. Oppen, K. Kanisawa, S. Fölsch, Physical Review B 2022, 105, [eLocator: 125418].

E. Sierda, X. Huang, D. I. Badrtdinov, B. Kiraly, E. J. Knol, G. C. Groenenboom, M. I. Katsnelson, M. Rösner, D. Wegner, A. A. Khajetoorians, Science 2023, 380, 1048.

J. Sun, T. Liu, Y. Du, H. Guo, Phys. Rev. B 2022, 106, [eLocator: 155417].