Making artificial $\textit{p}_{x,y}$-orbital honeycomb electron lattice on metal surface

TL;DR

This paper proposes a method to create and analyze a $p_{x,y}$-orbital honeycomb electron lattice on a Cu surface using STM, enabling exploration of exotic quantum properties in a controllable artificial system.

Contribution

It demonstrates a feasible approach to construct and identify $p$-orbital honeycomb lattices on metal surfaces, combining theoretical models and STM measurements.

Findings

Identification of $p$-orbital surface bands via STM LDOS patterns

Discovery of two types of edge states in the $p$-orbital honeycomb lattice

Proposal of potential exotic quantum phenomena such as ferromagnetism and QAH effect

Abstract

We theoretically demonstrate that the desired $p_{x,y}$-orbital honeycomb electron lattice can be readily realized by arranging CO molecules into a hexagonal lattice on Cu(111) surface with scanning tunneling microscopy (STM). The electronic structure of the Cu surface states in the presence of CO molecules is calculated with various methods, \textit{i.e.}~DFT simulations, muffin-tin potential model and tight-binding model. Our calculations indicate that, by measuring the LDOS pattern using STM, the $p$-orbital surface bands can be immediately identified in experiment. We also give an analytic interpretation of the $p$-orbital LDOS pattern with $k \cdot p$ method. Meanwhile, different from the case of graphene, the $p$-orbital honeycomb lattice has two kinds of edge states, which can also be directly observed in STM experiment. Our work points out a feasible way to construct a $p_{x,y}$-orbital honeycomb electron lattice in a real system, which may have exotic properties, such as Wigner crystal, ferromagnetism, $f$-wave superconductivity, quantum anomalous Hall (QAH) effect. Furthermore, we also propose a simple way to calculate and identify the modified Cu surface bands in the Cu/CO systems with the DFT simulations. Considering the recent works about $p$-orbital square lattice in similar systems [M. R. Slot, \textit{et al.} Nat. Phys. \textbf{13}, 672 (2017); Liang Ma, \textit{et al.} Phys. Rev. B \textbf{99}, 205403 (2019)], our work once again illustrates that the artificial electron lattice on metal surface is an ideal platform to study the orbital physics in a controllable way.

Equations

$$H_{Cu}=-\frac{{\hslash }^2}{2m^{\ast }}{\nabla }^2+U(r).$$

$$\begin{array}{r}\text{LDOS}(r,\epsilon )=\sum _{mk\sigma }{|{\varphi }_{mk\sigma }(r)|}^2\delta (\epsilon -{\epsilon }_{mk\sigma }).\end{array}$$

$$U_{atom}(\text{r})=\{\begin{array}{r}-U_0,|\text{r}|\le \frac{L}{2}\\ 0,|\text{r}|>\frac{L}{2}\end{array}$$

$$\varphi (r,\theta )=\text{R}(r)\text{Y}(\theta ),$$

$$\begin{array}{r}\frac{d^2\text{Y}(\theta )}{d{\theta }^2}=-n^2\text{Y}(\theta )n=0,1,2,\cdots \end{array}$$

$$\{\begin{array}{r}{r}^2\frac{d^2\text{R}(r)}{d{r}^2}+r\frac{d\text{R}(r)}{dr}+({\lambda }_1{r}^2-n^2)\text{R}(r)=0,r\le \frac{L}{2}\\ r^2\frac{d^2\text{R}(r)}{d{r}^2}+r\frac{d\text{R}(r)}{dr}-({\lambda }_2{r}^2+n^2)\text{R}(r)=0,r>\frac{L}{2}\end{array}$$

$$\{\begin{array}{rl}& c_1{J}_n(\sqrt{{\lambda }_1}{r}_0)=c_2{K}_n(\sqrt{{\lambda }_2}{r}_0)\\ & c_1{J}_n^{{}^{\prime }}(\sqrt{{\lambda }_1}{r}_0)\sqrt{{\lambda }_1}=c_2{K}_n^{{}^{\prime }}(\sqrt{{\lambda }_2}{r}_0)\sqrt{{\lambda }_2}\\ & |c_1{|}^2{\int }_0^{r_0}r{J_n}^2(\sqrt{{\lambda }_1}r)dr+|c_2{|}^2{\int }_{r_0}^{\infty }r{K_n}^2(\sqrt{{\lambda }_2}r)dr=1\end{array}$$

$$\begin{array}{r}{\varphi }_{p_x}(r,\theta )=\frac{1}{\sqrt{\pi }}\text{R}(r)\cos \theta \\ {\varphi }_{p_y}(r,\theta )=\frac{1}{\sqrt{\pi }}\text{R}(r)\sin \theta \end{array}$$

$$\text{R}(r)=\{\begin{array}{r}{c}_1{\text{J}}_1(\sqrt{{\lambda }_1}r),r\le \frac{L}{2}\\ c_2{\text{K}}_1(\sqrt{{\lambda }_2}r),r>\frac{L}{2}\end{array}$$

$$H_{\sigma }=t_{\parallel }\sum _{\vec{r}\in A,i}\{p_{\vec{r},i}^{†}{p}_{\vec{r}+{\hat{l}}_i,i}+h.c.\},$$

$$H_{\pi }=-t_{\perp }\sum _{\vec{r}\in A,i}\{p_{\vec{r},i}^{\prime †}{p}_{\vec{r}+{\hat{l}}_i,i}^{\prime }+h.c.\}.$$

$$\displaystyle H(k)=t_{\parallel}\left(\begin{array}[]{cc}0&H_{\parallel}\\ H^{+}_{\parallel}&0\\ \end{array}\right)-t_{\perp}\left(\begin{array}[]{cc}0&H_{\perp}\\ H^{+}_{\perp}&0\\ \end{array}\right),$$

$$\displaystyle H_{\parallel}=\left(\begin{array}[]{cc}\frac{3}{4}(e^{i\textbf{k}\textbf{l}_{1}}+e^{i\textbf{k}\textbf{l}_{2}})&\frac{\sqrt{3}}{4}(e^{i\textbf{k}\textbf{l}_{1}}-e^{i\textbf{k}\textbf{l}_{2}})\\ \frac{\sqrt{3}}{4}(e^{i\textbf{k}\textbf{l}_{1}}-e^{i\textbf{k}\textbf{l}_{2}})&\frac{1}{4}(e^{i\textbf{k}\textbf{l}_{1}}+e^{i\textbf{k}\textbf{l}_{2}})+e^{i\textbf{k}\textbf{l}_{3}}\\ \end{array}\right),$$

$$\displaystyle H_{\perp}=\left(\begin{array}[]{cc}\frac{1}{4}(e^{i\textbf{k}\textbf{l}_{1}}+e^{i\textbf{k}\textbf{l}_{2}})+e^{i\textbf{k}\textbf{l}_{3}}&\frac{\sqrt{3}}{4}(e^{i\textbf{k}\textbf{l}_{2}}-e^{i\textbf{k}\textbf{l}_{1}})\\ \frac{\sqrt{3}}{4}(e^{i\textbf{k}\textbf{l}_{2}}-e^{i\textbf{k}\textbf{l}_{1}})&\frac{3}{4}(e^{i\textbf{k}\textbf{l}_{1}}+e^{i\textbf{k}\textbf{l}_{2}})\\ \end{array}\right).$$

$${\psi }_{p_{x,y}}^{A,B}(\text{k},\text{r})=\frac{1}{\sqrt{N}}\sum _{{\text{R}}_{A,B}}{\varphi }_{p_{x,y}}^{A,B}(\text{r}-{\text{R}}_{A,B})e^{i\text{k}{\text{R}}_{A,B}},$$

$${\Psi }_{m\text{k}}(\text{r})=$$

$$+$$

$$\text{LDOS}(\epsilon ,\text{r})=\sum _{m\text{k}\sigma }|{\Psi }_{m\text{k}\sigma }(\text{r}){|}^2\delta (\epsilon -{\epsilon }_{m\text{k}\sigma }).$$

$${\Psi }_{m\text{k}}(\text{r})=$$

$$+$$

$${\Psi }_{1\text{k}}(\text{r})=$$

$$+$$

$${\Psi }_{2\text{k}}(\text{r})=$$

$$+$$

$${\Psi }_{3\text{k}}(\text{r})=$$

$$+$$

$${\Psi }_{4\text{k}}(\text{r})=$$

$$+$$

$$\underset{\Delta E\to 0}{\lim }\frac{{\rho }_{sf}^{3D}\times \Delta h}{\Delta E}={\rho }_{sf}^{2D},$$

Full text

I introduction

Orbital is an independent degree of freedom of electrons in solid in addition to their charge and spin. In many transition metal oxides, orbital physics plays an important role, and the interplay between the orbital, spin and charge degree of freedom can induce many exotic phenomena such as colossal magnetoresistance Millis et al. (1996); Ramirez (1997); Ling et al. (2000); Mathieu et al. (2004), superconductivity Bednorz and Müller (1986); Luke et al. (1998); Ohtomo and Hwang (2004); Kamihara et al. (2006), and metal-insulator transition Van Aken et al. (2003); Chen and Zou (2007); Kong et al. (2008). Though the orbital physics is important, the study about the influence of orbital degree of freedom in real materials is still a big challenge, because orbital is always coupled with other degrees of freedom such as charge, spin, or crystal field. Artificial lattice systems, e.g. cold atoms in optical lattice Müller et al. (2007), photonic lattice Milićević et al. (2017), offer ideal platforms to investigate the orbital physics due to their unprecedented controllability. In last decade, great efforts have been made to simulate the orbital physics in optical lattices Wu (2009); Li and Liu (2016). In experiment, bosons in p-band optical lattice have been realized and intensively studied Wirth et al. (2011); Ölschläger et al. (2013); Kock et al. (2015).

Interestingly, some exotic orbital related quantum states, which do not exist in real materials, can be realized in artificial lattice systems. An intriguing example is the $\textit{p}_{x,y}$-orbital honeycomb lattice, which was proposed by Wu and coworkers in 2007 Wu et al. (2007); Wu and Das Sarma (2008). Different from graphene, in optical lattice, the energy of $s$ orbital is rather largely separated away from the $p$ orbital, so that a p-orbital honeycomb lattice can be constructed without $sp$ hybridization. The p-orbital honeycomb lattice has a unique band structure. It has four $p$ bands, where two of them are flat bands and the other two give rise to a Dirac cone at the $K$ points in Brillouin zone. Due to the quenched kinetic energy of flat bands, strong correlation effects, such as Wigner crystal Wu et al. (2007), ferromagnetism Zhang et al. (2010), may appear. Meanwhile, the band structure is able to be tuned into topological nontrivial state Zhang et al. (2019); Zhu et al. (2019), e.g. the QAH effect Wu (2008); Zhang et al. (2011, 2014); Milićević et al. (2019); Chen and Xie (2019)and QSH effect Li et al. (2018); Canonico et al. (2019). It also enables unconventional $f$-wave Cooper pairing with a conventional attractive interaction Lee et al. (2010). Beyond the bulk properties, this $p_{x,y}$-orbital honeycomb lattice has two kinds of edge states, which are distinct from that of graphene. One is zero energy edge states, which have a similar origin to the conventional edge states in graphene. The other is novel dispersive edge states. Both of them have been observed in photonic lattice in a recent experiment Milićević et al. (2017).

However, fermions in p-orbital honeycomb lattice, which have been intensively studied in theory, have not been reported in experiment. Very recently, it was also theoretically proposed that the $p_{x,y}$-orbital honeycomb lattice can be founded in some special materials Zhou et al. (2016); Barreteau et al. (2017); Song et al. (2018); Jiang et al. (2019); Zhou et al. (2019); Hu et al. (2020), but these proposals have not been confirmed in experiment so far. To the best of our knowledge, experimental realization of this interesting $p$-orbital honeycomb lattice has only been reported in photonic lattice of coupled micropillars Milićević et al. (2017); Jacqmin et al. (2014). Recently, we propose that $p$-orbital bands can be realized in the artificial electron lattice on metal surface Ma et al. (2019), where the metal surface electrons are transformed into an electron lattice by periodically arranged adatoms. And the $p$-band features on square and Lieb electron lattice have already been observed in a recent scanning tunneling microscopy (STM) experiment Slot et al. (2017); Ma et al. (2019).

In this work, we illustrate that the interesting $p_{x,y}$-orbital honeycomb electron lattice can be readily constructed on metal surface with the same technique by designing a proper adatom pattern. And it actually has been realized in experiment Gomes et al. (2012); Wang et al. (2014), only further measurement is needed to give a confirmation. The scheme is illustrated in Fig. 1. Here, the CO molecules are periodically deposited on Cu(111) surface with the help of the STM tip. This Cu/CO system has been successfully used to construct the artificial honeycomb Park and Louie (2009); Gomes et al. (2012); Wang et al. (2014), square Slot et al. (2017), Lieb electron lattices Slot et al. (2017) on metal surface. In Fig. 1 (a), CO molecules are represented by black balls, where each applies a repulsive potential on metal surface. The repulsive potential can be approximately described as a muffin-tin potential, which is denoted as the blue circles in Fig. 1 (a). As pointed out in Ref. Park and Louie (2009); Gomes et al. (2012); Wang et al. (2014), if the CO molecules are arranged into a hexagonal lattice (Fig. 1), the Cu surface states can be forced into a honeycomb electron lattice because of the repulsive potential applied by the CO molecules. Former works Park and Louie (2009); Gomes et al. (2012); Wang et al. (2014) only consider the lowest two energy bands, which correspond to the $s$ orbital and has a graphene like band structure. We would like to point out that, in the same system, the higher bands are from the artificial $p$ orbitals, which is just the desired $p_{x,y}$-orbital honeycomb electron lattice [Fig. 1 (a),(c)]. As shown in Fig. 1 (b) and (d), the $p_{x,y}$ orbitals will form $\sigma$ bonds parallel to the bond direction, and $\pi$ bonds perpendicular to the bond direction as well Wu et al. (2007); Wu and Das Sarma (2008). The LDOS of the $p$ orbitals can be directly observed by STM. We use density functional theory (DFT) simulations, muffin-tin potential model, tight-binding (TB) model and $k\cdot p$ method to interpret this $p$-orbital picture. We also illustrate that the two kinds of edge states of this $p$-orbital honeycomb lattice could be directly observed in this artificial electron lattice system.

To calculate the modified Cu(111) surface bands in the presence of CO adatoms, the muffin-tin model is the most common method, while the DFT simulations do not work well so far Ropo et al. (2014a); Paavilainen et al. (2016). The reason is that, in DFT simulations, we have to consider a large supercell, including a thick Cu slab and CO adatoms, in which the surface bands and the Cu bulk bands are heavily folded and thoroughly mixed together. Thus, the major obstacle for the DFT simulations is to distinguish the modified surface bands from massive bulk bands. Here, we illustrate that the modified Cu surface bands can be identified by proper layer resolved band structure projections, so that now the desired Cu surface bands in the Cu/CO systems can be directly calculated using the DFT simulations. It should be a general method to calculate the adatom modified metal surface states via DFT simulations.

The outline of this paper is as follows: in Sec. II, we give the models and methods used in the calculations; in Sec. III, we show the numerical results and the corresponding discussions; finally, a short summary is given in Sec IV.

II model and methods

Here, we introduce the theoretical methods, i.e. muffin-tin potential model, tight-binding model and DFT simulations, which we used here to calculate the modified Cu surface bands in the Cu/CO systems.

II.1 Muffin-tin potential model

As mentioned above, the Cu surface states in Cu/CO system can be intuitively described by a muffin-tin potential model. We illustrate this muffin-tin potential model in Fig. 2 (a). Here, CO molecules can be approximately viewed as a repulsive muffin-tin potential $U(r)$ (black discs), which is $U_{\textrm{CO}}>0$ inside the black discs and zero elsewhere. Thus, the surface electrons are confined into the regions in between CO molecules, which actually gives rise to a honeycomb lattice. The sites of the honeycomb lattice are represented by the orange discs in Fig. 2 (a). The Hamiltonian of the muffin-tin model is

[TABLE]

This Hamiltonian can be solved with plane wave method to get the energy bands of Cu surface states in the presence of CO molecules Park and Louie (2009); Li et al. (2016); Qiu et al. (2016). Here, $m^{*}=0.38m_{0}$ is the effective mass of Cu surface states, where $m_{0}$ is the electron mass Gomes et al. (2012). The parameters of the muffin-tin potential $U(r)$ for the Cu/CO system are $U_{\textrm{CO}}=$ 0.8 eV and the diameter $d=$ 0.5 nm, which are got by fitting the DFT results, as will be shown later.

The muffin-tin potential of a single CO adatom has a fixed value. However, we can get a stronger effective potential applied on Cu surface via replacing the single CO adatom with a CO cluster Slot et al. (2019). For example, several kinds of CO clusters with different sizes are plotted in Fig. 2 (b). We use cluster-$l$ to label the CO cluster where the side length of hexagon is $l$ times the length of the Cu-Cu distance. Note that the larger the cluster size is, the stronger influence of the effective potential is.

The LDOS can be got by

[TABLE]

We will show that the low energy bands got by the muffin-tin potential model can be well interpreted as the $s$ and $p$ bands of a honeycomb lattice.

II.2 Orbitals in artificial atom

We then discuss the concept of orbitals in the artificial electron lattice. It should be noted that each site (orange disc) of this honeycomb lattice in Fig. 2 (a) actually can be viewed as an artificial atom. In each artificial atom, the electrons are confined in the region enclosed by the adjacent CO molecules. Considering this potential confinement, the eigenstates of an isolated artificial atom are discrete, which are very similar to the orbitals of a real atom. Then, through hopping between these artificial orbitals, we get the energy bands of this artificial electron lattice.

However, this orbital picture (and TB model) does not always work. This is because that the confinement potential is finite. If the electron kinetic energy is large enough, surface electrons can not be confined. Thus, the surface states should be described by the nearly free electron model instead of the TB model. In this situation, the concept of artificial orbitals is invalid.

We can roughly use a cylindrical potential well $U_{atom}$ to describe this two dimensional artificial atom. $H_{atom}=-\frac{\hbar^{2}}{2m^{*}}{\nabla}^{2}+U_{atom}(\textbf{r})$, where

[TABLE]

$U_{0}$ and L are the value and diameter of the potential well, respectively. One artificial atom corresponds to one cylindrical potential well, which is represented by one orange disc in Fig. 2 (a). The electrons trapped in $U_{atom}$ form the orbitals, and the hopping between orbitals gives rise to TB bands.

The orbitals can be got by calculating the bound states in the potential well. The eigenfunction of $H_{atom}$ is

[TABLE]

where $r$ is the radial coordinate and $\theta$ is the polar angle. We have

[TABLE]

where, $\textrm{Y}_{n}(\theta)=a_{n}\cos(n\theta)+b_{n}\sin(n\theta)$ for $n>0$, and $Y_{0}(\theta)=1/\sqrt{2\pi}$. Here, $a_{n}$ and $b_{n}$ are the coefficients, which can take on any value (but the wavefunction should be normalized). The equations of radial part are

[TABLE]

where $\lambda_{1}=\frac{2m(\varepsilon+U_{0})}{\hbar^{2}}$, $\lambda_{2}=-\frac{2m\varepsilon}{\hbar^{2}}$. At the boundary $r_{0}=\frac{L}{2}$, the continuity of the wavefunction and its derivative should be satisfied,

[TABLE]

Here, the last equation is the requirement of wave function normalization, $c_{1}$ and $c_{2}$ are two coefficients to be determined. $J_{n}$ ($K_{n}$) is the Bessel (Hankel) function. Together with this boundary condition, we can get the eigenvalues and eigenfunctions of the bound states in this potential well.

For example, in the case of $U_{0}=0.9$ eV and $L=1.8$ nm, there is one bound state with $n=0$, i.e. the $s$ orbital. With $n=1$, the allowed $p$ orbitals are

[TABLE]

where

[TABLE]

The coefficients are $c_{1}\approx 2.20$, $c_{2}\approx 1.16$, and the energy of the $p_{x,y}$ orbitals is $\varepsilon_{p_{x}}=\varepsilon_{p_{y}}=\varepsilon\approx-0.074$ eV. However, no solution exists for $n>1$, which implies that only $s$ and $p$ orbitals are valid in this situation. Note that, the larger $U_{0}$ is, the more artificial orbitals are allowed.

The value of $U_{atom}$ for the Cu/CO system can be estimated by fitting the bands calculated from the muffin-tin model. For the case in the experiment Gomes et al. (2012), a reasonable $U_{atom}$ is: $U_{0}=0.9$ eV, $L=1.8$ nm, which is just the above example. So, we argue that only $s$ and $p$ orbitals are valid in such Cu/CO systems.

We illustrate the artificial orbitals in the honeycomb lattice in Fig. 1 (a). Red discs represent the $s$ orbitals, and green lobes denote the $p_{x,y}$ orbitals. Different from graphene, the $s$ orbitals are far away from the $p$ orbitals in energy, so that there is no $sp$ hybridization here. The $s$ orbitals give rise to an artificial graphene, which has already been observed in experiment Gomes et al. (2012); Wang et al. (2014). In the following, we focus on the $p$ orbitals. Note that in two dimensional electron system, there are only $p_{x,y}$ orbitals, and no $p_{z}$ orbital. Thus, we actually get a $p_{x,y}$-orbital honeycomb lattice, which is illustrated in Fig. 1 (a) and (c). The honeycomb lattice has two sublattices, A and B. So, we have four $p$ bands above the two $s$ bands in this artificial honeycomb lattice.

Note that the cylindrical potential well is a rough approximation of the confinement potential applied by CO molecules, which is just used to illustrate the concept of artificial orbitals. As shown in Fig. 2 (a), the regions where the surface electrons are confined are actually anisotropic and have three-fold rotational symmetry. Such anisotropy of the confinement potential is not included in the cylindrical (TB) model. This does not obviously influence the $s$ bands, but will give rise to a difference between the cylindrical (TB) band structure of $p$ orbitals and that got from muffin-tin model.

II.3 Tight-binding model

The tight-binding model of $p_{x,y}$-orbital honeycomb lattice is given in Ref. Wu and Das Sarma (2008). Here, we give a short introduction to this TB model. First, we define three vectors $\vec{l}_{1,2,3}$ [see in Fig. 1 (c)], which point from A site to three nearest neighbor B sites, i.e., the directions of bonds. $a_{0}$ is the lattice constant [see in Fig. 1 (a)]. Considering the lattice geometry, $p_{x,y}$ orbitals are projected to directions parallel to $\vec{l}_{i=1,2,3}$ forming the projected orbital $p_{i}=(p_{x}\vec{e}_{x}+p_{y}\vec{e}_{y})\cdot\frac{\sqrt{3}\vec{l}_{i}}{a_{0}}$, where $\vec{e}_{x}$ ($\vec{e}_{y}$) is the unit vector along x (y) direction Wu and Das Sarma (2008). We also define $p^{\prime}_{i=1,2,3}$ as the projected $p$ orbitals perpendicular to $\vec{l}_{i=1,2,3}$. Note that here only two of $p_{i}$ or $p^{\prime}_{i}$ orbitals are linearly independent. Bonds exist between these projected orbitals. There are two kinds of bonds between the $p$ orbitals on neighboring sites, i.e., $\sigma$ and $\pi$ bond. As shown in Fig. 1 (d), $\sigma$ bond corresponds to the hopping between $p_{i}$ orbitals (“head to tail”) , and $\pi$ bond corresponds to the hopping between $p^{\prime}_{i}$ orbitals (“shoulder by shoulder”). The Hamiltonian is $H=H_{\sigma}+H_{\pi}$, where

[TABLE]

[TABLE]

Here, $t_{\parallel}$ ($t_{\perp}$) is the hopping of $\sigma$ ($\pi$) bond. The Hamiltonian can be viewed as the summation of all $\sigma$ and $\pi$ bonds.

It is convenient to use the basis $[p^{A}_{x},p^{A}_{y},p^{B}_{x},p^{B}_{y}]$, the corresponding Hamiltonian matrix is,

[TABLE]

where

[TABLE]

[TABLE]

Same expression can be obtained if we apply the Slater-Koster formalism Slater and Koster (1954) to $p_{x,y}$ orbitals. Diagonalizing $H(k)$, we can get the TB energy bands and wavefunctions. Using the $p$-orbital basis in Eq. (8), we can further calculate the LDOS based on the TB model. The details to calculate the LDOS are given in Sec. III.3. Meanwhile, we can also get the edge states by calculating the bands of the ribbon structure based on this TB model Milićević et al. (2017).

II.4 DFT simulations

In order to give more information about the experiment, DFT simulations can be done where more details of the Cu/CO system are included Ropo et al. (2014b); Paavilainen et al. (2016); Ma et al. (2019). As will be shown later, our DFT simulations are also in good agreement with the muffin-tin potential model, which means that the muffin-tin potential model captures the essential features of the experimental observations.

The DFT simulations use the Vienna Ab-initio Simulation Package (VASP), where projector-augmented wave method and a plane wave basis set are used Kresse and Furthmüller (1996). We select the Perdew-Burke-Ernzerhof (PBE) version of the generalized gradient approximation Perdew et al. (1996). The plane wave cutoff energy is $400$ eV and the $k$-point mesh is $9\times 9\times 1$. We use a Cu slab to mimic the Cu(111) surface, see in Fig. 1. The lattice constant $a_{0}$ is the distance between two adjacent CO molecules. Here, $a_{0}$ is eight times the length of the Cu-Cu distance (about 2.04 nm), namely an $8\times 8$ supercell is used. To avoid the interaction between nearest cells we set the inter-cell vacuum space to be $22{\textrm{\AA}}$. CO molecules are adsorbed on the top of Cu atoms, and the Cu-C distance and C-O distance are set to be $1.85{\textrm{\AA}}$ and $1.155{\textrm{\AA}}$, respectively. We only optimize the positions of C, O atoms and the Cu atoms closest to C atoms. The convergence criteria for force acting on each atom is $0.02\textrm{eV}/{\textrm{\AA}}$. This has been shown to be accurate enough to capture the physics discussed here Ma et al. (2019).

We first calculate the energy bands of the Cu/CO systems, and then identify the modified surface bands via proper layer-resolved band projections. We also simulate the LDOS of the surface states and the corresponding STM image with the Tersoff-Hamann scheme Tersoff and Hamann (1983). Finally, we calculate the charge difference caused by CO adsorption using the VASPKIT code Wang et al. (2019), in order to study the effect of charge transfer between CO molecules and Cu slab.

III Results and Discussions

III.1 Energy bands

Various theoretical methods, i.e. DFT simulations, muffin-tin model and TB model, are used to calculate the modified Cu surface bands in the Cu/CO systems. The muffin-tin model and TB model can give an intuitive picture about the $p$ bands, which clearly interprets the results of DFT simulations. Meanwhile, the parameters of the muffin-tin model and TB model can be got by fitting the results of DFT simulations.

III.1.1 DFT simulations

We first discuss how to calculate the Cu surface bands via DFT simulations. Here, we use a Cu slab to mimic the Cu surface. As mentioned above, the challenge here is to distinguish the surface bands from the bulk bands of the Cu slab. In the following, we use three steps to illustrate the Cu surface bands in the Cu/CO system.

At the first step, we illustrate how to identify the Shockley surface bands of a Cu slab in the DFT simulations. In Fig. 3 (a), we plot the energy bands of a 51-layer-thick pure Cu slab, in which the surface bands are well known in former literature Butti et al. (2005); Berland et al. (2012); Courths et al. (2001); Kevan (1983); Jeandupeux et al. (1999); Tamai et al. (2013). As we see here, in such Cu slab, the surface bands are mixed up with the bulk bands. In order to identify the surface bands, we project all the bands of the Cu slab into the surface layers by layer resolved band structure projections. The most important finding is that the Shockley surface states are mainly composed of the $p_{z}$ orbitals of surface layers, i.e., Cu $p_{z}$ orbitals of the top two layers, see the red dotted lines in Fig. 3 (a). Note that the red dotted lines here are just the surface bands reported in former literatures Butti et al. (2005); Berland et al. (2012); Courths et al. (2001); Kevan (1983); Jeandupeux et al. (1999); Tamai et al. (2013), and these states are strongly localized at the surface layers (See Fig. 11 (h) and (i) in Appendix A). The effective mass of the surface bands is about $m^{*}=0.45m_{0}$ (the blue solid line is the fitted $k^{2}$ dispersion). Thus, we can use this property to identify the surface bands in the DFT simulations.

A key fact here is that the Cu slab has two surfaces, which gives rise to two degenerate surface bands in the DFT simulations. If the thickness of Cu slab is largely reduced, the two degenerate surface bands will be split due to the finite size effect. In Fig. 3 (b), we plot the bands of an 8-layer-thick pure Cu slab. With the same method, we get two split surface bands instead (red dotted lines), due to the coupling between the top and bottom surface states. The band split is about 0.2 eV. Such surface band split is unreal and does not exist in reality, since the sample of Cu(111) surface is thick enough. In principle, it can be avoided by using a thick Cu slab in the calculations. However, when we consider the CO adatoms, a very large supercell is needed for the Cu slab, so that the thickness of the Cu slab in the DFT simulations should be limited. For example, if we use an $8\times 8$ supercell ($a_{0}=2.04$ nm) and an 8-layer-thick Cu(111) slab, we have 512 Cu atoms in one supercell which has already reached the limit of our calculation resource. Fortunately, except for the band split, such coupling only slightly changes the effective mass of the surface bands in this situation ($m^{*}=0.5m_{0}$ and blue solid line is the fitted $k^{2}$ dispersion). It can still be viewed as the free electron band with $k^{2}$ dispersion. Therefore, in order to compare with experiment, a reasonable and feasible way is to choose an 8-layer-thick Cu slab (maximumly allowable in our DFT simulations), and artificially select one of the two calculated surface bands.

At the second step, we show the folded surface bands where a supercell of the Cu slab is considered. The surface band of Cu(111) surface has a $k^{2}$ dispersion near the $\Gamma$ point, which can be approximately viewed as a free electron gas. But if the CO molecules are periodically adsorbed, a supercell of the Cu slab has to be used as illustrated in Fig. 1 (a). In this situation, the surface bands are folded into a small BZ defined by the supercell. For example, Fig. 3 (c) gives the folded surface bands (red dots) of the 8-layer-thick Cu slab with an $8\times 8$ supercell. Note that here we only consider a pure Cu(111) surface and CO molecules have not been included yet. Since the 8-layer-thick Cu slab has two surface bands, we have two sets of folded surface bands here, which have an energy split about $0.2$ eV. We select one set of the folded surface bands (low energy one), see the red solid lines in Fig. 3 (c). Such folded surface band calculated by DFT simulations can be well understood by the muffin-tin model in Eq. (1). Here, we set $U_{\textrm{CO}}=0$ in the muffin-tin model, since CO molecules have not been adsorbed yet. Thus, the muffin-tin model gives the folded free electron bands here. In Fig. 3 (d), the red dots are the replotted folded surface bands got by DFT simulations, i.e., the red solid lines in Fig. 3 (c), and the blue solid lines are the folded surface bands got by the muffin-tin model ($m^{*}=0.50m_{0}$ and $U_{\textrm{CO}}=0$). We see that they coincide well with each other, especially in the low energy region. We argue that the small discrepancy is because that the assumption of $k^{2}$ dispersion is invalid for the DFT surface bands at high energy, as shown in Fig. 3 (a) and (b).

At the last step, we calculate the modified Cu surface bands in the presence of CO adatoms via DFT simulations. In Fig. 3 (e), we show the bands of Cu/CO system, where an $8\times 8$ supercell is used and CO adatoms are included as shown in Fig. 1 (a). Here, the red dots are the calculated surface bands via the layer projection method. Similarly, we select one of the two sets of surface bands, which are replotted in Fig. 3 (f). The blue lines are the surface bands calculated by muffin-tin model with $U_{\textrm{CO}}=0.8$ eV. The consistency of the DFT bands and that got from the muffin-tin model clearly indicates that the muffin-tin model gives a rather well description about the surface bands of the Cu/CO system. And here we also get the value of muffin-tin potential parameters, i.e. $U_{\textrm{CO}}=0.8$ eV and $d=0.5$ nm, by fitting the DFT bands. It is clear to see the influence of the CO adatoms by comparing Fig. 3 (f) to (d). In Fig. 3 (f), the CO potential lifts the band degeneracy along some high symmetry lines, and a Dirac point at the lowest two bands emerges at $K$ point Park and Louie (2009); Nádvorník et al. (2012), denoted by the blue circle. Actually, the lowest two bands here are just the $s$ bands of the honeycomb lattice, which give rise to an artificial graphene. Above the lowest two bands, there is another Dirac point at $K$ point, also denoted by a blue circle. As will be shown later, these high energy bands above the $s$ bands, as well as the Dirac point, just result from the $p_{x,y}$ orbitals of the honeycomb lattice.

III.1.2 Muffin-tin and TB models

Through DFT simulations, we have shown that the modified Cu surface bands in Cu/CO systems can be well described by the muffin-tin model. Here, based on the muffin-tin and TB models, we then illustrate that such Cu/CO systems do host a $p_{x,y}$ honeycomb electron lattice.

The surface bands calculated by muffin-tin model are plotted in Fig. 4 (a). Here, to demonstrate, we use an ideal set of parameters for the muffin-tin model, $U_{\textrm{CO}}=9$ eV and $d=0.5$ nm. As discussed in former literature Ma et al. (2019); Slot et al. (2019), the lowest two bands (cyan lines) are from the $s$ orbitals, which are similar to the bands of graphene. This is just the reason why this artificial structure is named as “artificial graphene” Park and Louie (2009). We are interested in the upper four bands (blue lines). Actually, the four bands here are from the $p_{x,y}$ orbitals on the honeycomb lattice. To prove this point, we use the TB model in Eq. (12) to fit the $p$-orbital energy bands in Fig. 4 (b). The red dashed lines are the results of TB model, and we also replot the $p$ bands from muffin-tin potential model (blue lines) as a comparison. We see that the $p$ bands got from the TB model coincide well with that from the muffin-tin potential model in the low energy region. The fitted parameters of the TB model are : $t_{\parallel}=0.49$ eV, $t_{\perp}=0.15$ eV and $\varepsilon_{p}=1.2$ eV. At high energy, there is an obvious discrepancy. There are two reasons. One is due to the anisotropy of the confinement potential as mentioned in the last section, which has been ignored in the TB model. It does not obviously influence the $s$ bands, but can not be completely ignored for the $p$ bands. The other is due to hybridization between the $p$ orbitals and the states with higher energy. Since the confine potential is finite, the surface states at higher energy here are more like that of nearly free electrons, for which the orbital picture and TB model do not work.

To make the $p$-orbital picture clearer, we consider an extreme case with $U_{\textrm{CO}}=9$ eV, $d=1.6$ nm. In this case, confinement potential is extremely large. Thus, more orbitals appear and these orbitals are well separated in energy. We plot the bands in this extreme case in Fig. 4 (c), (d). The lowest two bands are from $s$ orbitals. The upper four bands are just from the $p_{x,y}$ orbitals. Now, we see that there are two nearly flat bands and a Dirac point at the K point, which are the characteristics of the bands of $p_{x,y}$-orbital honeycomb lattice Wu et al. (2007). We also give a fitting with TB model, as shown in Fig. 4 (d). Now, there is no hybridization with higher bands, and the discrepancy of the band structure at high energy should only result from the influence of the anisotropy of the confinement potential.

Interestingly, in Fig. 4 (c), the two bands above the $p$ bands are similar to the lowest two $s$ bands (black lines). Actually, they correspond to a new $s$ orbital of the artificial atom, which appears due to the very large repulsive (confinement) potential in this extreme case. So, we get another two graphene-like $s$ bands at high energy. However, they have distinct LDOS patterns, which will be shown later.

For the Cu/CO system, we get $U_{\textrm{CO}}=0.8$ eV by fitting the DFT bands. However, the muffin-tin model tells that such repulsive potential is not strong enough to completely separate the $p$ bands from the $s$ bands in energy, and may not give a clear signal of $p$ band in experiment. In Fig. 5 (a), we plot the muffin-tin bands calculated with $U_{\textrm{CO}}=0.8$ eV (green solid lines) and $U_{\textrm{CO}}=9$ eV (red dashed lines), respectively. We see that, with a weak repulsive potential, e.g. $U_{\textrm{CO}}=0.8$ eV, the lowest $p$ band is badly bent, and obviously deviates from the typical shape of $p$ band. It implies that, with a single CO adatom, the $p$ band features in experiment are not quite evident. A better way is to use a CO cluster instead of a single CO adatom, as shown in Fig. 2 (b), which has been successfully utilized in experiment Slot et al. (2019). The reason is that the CO cluster is essentially equivalent to a stronger effective potential in muffin-tin model. In Fig. 5 (a), we also plot the muffin-tin bands calculated with the cluster-1 of CO, see the blue solid lines. It obviously demonstrates that, with a CO cluster, we can get the same surface band structures as that from a strong repulsive potential (i.e., single adatom with $U_{\textrm{CO}}=9$ eV). So, in order to give a more clear $p$ band characteristic, unless specified otherwise, we use the cluster-1 of CO in the muffin-tin model instead of a single CO adatom. This should be a more appropriate experimental design to observe the $p$ orbital bands in Cu/CO system. We emphasize that the physical picture of $p$ orbital remains the same no matter whether a single CO adatom or a CO cluster is used.

In Fig. 5 (b), we use TB model to fit the muffin-tin bands of the cluster-1 case. The red (green) dashed lines are the TB bands of $p$ ($s$) orbitals.

III.2 LDOS pattern: numerical results

Now, we discuss the LDOS pattern in this artificial honeycomb lattice, which can be directly measured using STM. The observed LDOS pattern can be used to identify the $p$ bands in experiment.

Fig. 6 (a-e) are the LDOS calculated with the TB model at $E-E_{f}=$ -0.05, 0.36, 0.66, 1.59, 1.79 eV, respectively. The corresponding TB bands have been given in Fig. 5 (b), and the energy of the LDOS are denoted by red arrows. In Fig. 6 (a), the LDOS corresponds to an energy at $s$ band ($E-E_{f}$= -0.05 eV). It clearly shows a honeycomb lattice pattern Park and Louie (2009); Gomes et al. (2012), in which the disc at each site represents an onsite $s$ orbital. Fig. 6 (b-c) are the LDOS of the bottom $p$-orbital flat band. Note that the flat band here is not exactly flat, so that the LDOS at different energy of this band will show distinct features. For example, Fig. 6 (b) gives the LDOS at energy around the K point, while Fig. 6 (c) is around the $\Gamma$ point. Near the $\Gamma$ point, the pattern shows a bonding state feature, in which the electron density mainly distributes in between the two lattice sites and is tiny near the center of the artificial atom. In contrast, near the K point the distribution becomes more complex and nonlocal, as shown in Fig. 6 (b). As for the middle two $p$ bands, we find two different kinds of LDOS patterns. At lower energy, the LDOS shows a bonding state feature quite like that in Fig. 6 (c). At higher energy, the LDOS has an anti-bonding state feature, which has a node in the middle of the bond, as shown in Fig. 6 (d). The LDOS of the top flat band has a stronger anti-bonding state feature, as shown in Fig. 6 (e), in which the electron density inside of the bond is pushed to the opposite directions outside of the bond. We will give a further analysis of the above LDOS patterns with $k\cdot p$ method later in Sec. III.3.

Fig. 6 (f-j) are the LDOS calculated with the muffin-tin potential model at $E-E_{f}=$ -0.05, 0.36, 0.66, 0.84, 0.98 eV, respectively. The corresponding band structures are also given in Fig. 5 (b), and the energy of the LDOS are denoted by blue arrows. Note that the $s$ bands and the bottom $p$-orbital flat band coincide well with that from the TB model. At higher energy, the $p$ bands have a discrepancy with that from the TB model, as shown in Fig. 5 (b). So it is reasonable to expect that the muffin-tin LDOS patterns have similar shapes as the TB LDOS patterns at low energy, while they are different at high energy. Fig. 6 (f) shows the LDOS of the $s$ bands and Fig. 6 (g-h) show the LDOS of the bottom $p$-orbital flat band near the K point and $\Gamma$ point. They are qualitatively in agreement with the TB results. When we move to the middle $p$ bands, only the anti-bonding features can be observed, as shown in Fig. 6 (i-j). We think that the muffin-tin potential model is more realistic. We expect that the LDOS patterns observed in experiment should be more like that from the muffin-tin potential model.

The LDOS in Fig. 6 are for the Cu/CO system. As a comparison, if the muffin-tin potential is extremely large (e.g. $U_{\textrm{CO}}=$ 9 eV, $d=$ 1.6 nm), it is possible to get a standard band structure of $p_{x,y}$-orbital honeycomb lattice, see Fig. 4 (c-d). We then plot the LDOS of this extreme case in Fig. 7. The corresponding band structure is shown in Fig. 4 (c-d), and the energy of the LDOS are denoted by blue arrows. Fig. 7 (a) is the LDOS of the bottom $p$-orbital flat band. The flat band now is exactly flat, so that the LDOS of this band only shows a bonding state feature. As for the middle two $p$ bands, the LDOS in Fig. 7 (b-c) show bonding and anti-bonding state features, respectively, which are similar to the TB LDOS in Fig. 6 (c-d). The LDOS of the top flat band has a stronger anti-bonding state feature, much obvious than that in Fig. 7 (c).

An interesting case is the new $s$ bands above the $p$ bands in Fig. 4 (c). Though a similar graphene like band structure is found [black lines in Fig. 4 (c)], the new $s$ bands have a different LDOS pattern, as shown in Fig. 7 (d). It is because that the $s$ bands correspond to a new $s$ orbital with a node in the radial part of the wavefunction, which can be obtained by solving Eq. (5) and Eq. (6).

The LDOS can also be calculated with the DFT simulations. However, due to the finite thickness of Cu slab (8-layer-thick Cu slab) in the DFT simulations, the calculated LDOS contains the contributions from bulk states. It implies that the LDOS calculated with the DFT simulations are not the pure $p$ band LDOS as required. The muffin-tin model should be a better way. And we give the results and discussions about the LDOS from DFT simulations in Appendix A.

Finally, we argue that the above $p$ band LDOS can also be observed in other similar artificial honeycomb lattice systems, e.g., nano-patterned semiconductor heterostructure Nádvorník et al. (2012); Polini et al. (2013); Jacqmin et al. (2014); Klembt et al. (2018).

III.3 LDOS pattern: $k\cdot p$ method

After diagonalizing the $p_{xy}$-orbital tight-binding Hamiltonian Eq. (12), we have coefficients $c^{A}_{p_{x}}(m\textbf{k})$, $c^{A}_{p_{y}}(m\textbf{k})$, $c^{B}_{p_{x}}(m\textbf{k})$ and $c^{B}_{p_{y}}(m\textbf{k})$ for bloch wavefunctions $\{\psi^{A}_{p_{x}}(\textbf{k},\textbf{r}),\psi^{A}_{p_{y}}(\textbf{k},\textbf{r}),\psi^{B}_{p_{x}}(\textbf{k},\textbf{r}),\psi^{B}_{p_{y}}(\textbf{k},\textbf{r})\}$, which are the basis of TB Hamiltonian. Here, $m$ is the band index. The basis are

[TABLE]

where $\phi^{A,B}_{p_{x,y}}(r-R_{A,B})$ are just the wavefunctions of $p$ orbitals given in Eq. (8). Then, the eigenfunction of the $m$-th $p$ band is

[TABLE]

So, the TB LDOS can be obtained with

[TABLE]

Here, we show how to understand the bonding and antibonding features of the LDOS shown in Fig. 6 by analyzing the wavefunctions. The wavefunction in Eq. (16) can be divided into the components from A site and B site

[TABLE]

where the expression in the bracket of the first (second) line is from A (B) site. Near the $\Gamma$ point, the coefficients $c^{A,B}_{p_{x,y}}(m\textbf{k})$ can be got by $k\cdot p$ method Wu and Das Sarma (2008). For example, the wavefunction of the bottom $p$-orbital flat band ($m=1$) near the $\Gamma$ point can be written as

[TABLE]

We can see that the $p$ orbitals at A and B sites differ in a minus sign. This implies that the two $p$ orbitals are antiparallel (head to head), so that the LDOS in between A and B sites is amplified, i.e., the bonding state feature, as shown in Fig. 6 (c) and (h).

The wavefunctions of other three $p$ bands near the $\Gamma$ point can be written as

[TABLE]

[TABLE]

[TABLE]

Similarly, the second $p$ band of Eq. (20) also has a bonding state feature. While for the upper two $p$ bands ($m=3$ and $m=4$), the $p$ orbitals at A site and B site are parallel (head to tail), as shown in Eq. (21) and Eq. (22). Thus the LDOS in between A and B sites is canceled, i.e., the antibonding state feature, as shown in Fig. 6 (d) and (e).

But for LDOS near the K point, the wavefunctions of A and B sites do not have such special relation as that near the $\Gamma$ point Wu and Das Sarma (2008). Thus, we cannot observe a clear bonding or antibonding state feature, as shown in Fig. 6 (b) and (g).

III.4 Edge states

It is pointed out that the $p_{x,y}$-orbital honeycomb lattice has two kinds of edge states, which have been observed in photonic lattice system Milićević et al. (2017). Meanwhile, a recent STM experiment has detected the $s$-orbital edge states of the artificial graphene on Cu surface Wang et al. (2014). Thus, it should be ready to measure the edge states of the $p_{x,y}$-orbital honeycomb lattice in this Cu/CO system. We first calculate the bands of a zigzag ribbon of the $p_{x,y}$-orbital honeycomb lattice with the TB model. The width of the zigzag ribbon is about $34$ nm. The calculated bands are given in Fig. 8. Here, the green lines are the zero energy edge states, and the red lines are the dispersive edge states. These results are all in agreement with former work Milićević et al. (2017). Note that, in order to compare with Fig. 5 (b), we set $\varepsilon_{p_{x}}=\varepsilon_{p_{y}}=1.248$ eV in the TB model for the ribbon calculation.

In Fig. 9, we plot the LDOS of the zigzag ribbon near a zigzag edge, whose energy are denoted by red arrows in Fig. 8. Fig. 9 (a) and (c) are the LDOS plotted at $E-E_{f}=$ 0.55 eV and $E-E_{f}=$ 0.95 eV, respectively. These LDOS patterns correspond to bulk states, as shown in Fig. 8. So, the electron density distributes all over the bulk. At $E-E_{f}=$ 0.67 eV, the states are mainly the dispersive edge states (red lines in Fig. 8), but there are also some bulk states. So in Fig. 9 (b), we see that the intensity of LDOS is extremely large near the edge, which results from the dispersive edge states. Meanwhile, some bulk states can also be found in Fig. 9 (b). In Fig. 9 (d), we plot the LDOS at $E-E_{f}=$ 1.248 eV. As shown in Fig. 8, there are only the zero energy edge states at $E-E_{f}=$ 1.248 eV. Thus, we see that LDOS becomes nearly zero away from the zigzag edge. Our calculations indicate that we can distinguish different kinds of edge states by measuring the LDOS with STM.

III.5 Position of Fermi level

In STM experiment, to access the $p$ orbitals, the Fermi level should be close to $p$ bands. In the Cu/CO system, the position of the Fermi level relative to the surface bands depends on CO-CO (cluster-cluster) spacing $a_{0}$ Li et al. (2016); Qiu et al. (2016); Ma et al. (2019). In each unit cell of this artificial honeycomb lattice, the number of surface electrons is $S_{0}\textrm{N}_{e}$, where $S_{0}=\sqrt{3}a_{0}^{2}/2$ is the area of unit cell, $\textrm{N}_{e}\approx 0.72\,\textrm{nm}^{-2}$ is the electron density of Cu surface states. The surface bands can be got by the muffin-tin model. We know that the total number of surface electrons is just the summation of all the surface states below Fermi level. With this relation, we can calculate the position of Fermi level relative to the surface bands, as a function of the CO-CO spacing.

Based on the principle above, we give a simple relation to estimate the required CO-CO spacing, with which one can let the Fermi level cross the $m$-th surface band. We know that a filled band corresponds to 2 electrons per unit cell (including spin up and spin down). Thus, if the Fermi level crosses the $m$-th surface band, we have $S_{0}\textrm{N}_{e}\leq 2m<S_{0}\textrm{N}_{e}+2$. To let the Fermi level cross the $p$ bands, the relation above give a rough estimation about $a_{0}$, $2.5\,\textrm{nm}<a_{0}<4.4\,\textrm{nm}$, which is our suggested value of the further experiment. Note that, with different value of $a_{0}$, the surface band shape is similar, and only the band width is changed. This method to predict Fermi level was used in our previous work Qiu et al. (2016); Ma et al. (2019), and is in good agreement with the experimental observation Slot et al. (2017).

Here, we point out that the DFT method we used here can not correctly predict the position of the Fermi level, though the modified surface bands can be well achieved. For pure Cu(111) surface, the STM and ARPES experiments indicate that the Fermi level relative to the bottom of surface band $E_{f}-E_{0}$ is about 0.40$\sim$0.45 eV Kevan (1983); Jeandupeux et al. (1999); Tamai et al. (2013). However, DFT simulations with different methods give different values of Fermi level. For example, a full-potential linearized augmented plane-wave (FLAPW) method gives $E_{f}-E_{0}=$ 0.526 eV Butti et al. (2005); DACAPO code using the ultrasoft pseudopotential plane-wave method gives $E_{f}-E_{0}=$ 0.443 eV Berland et al. (2012); WIEN97 code using the relativistic self-consistent FLAPW method gives $E_{f}-E_{0}=$ 0.40 eV Courths et al. (2001); VASP package we used here based on projector-augmented wave method gives $E_{f}-E_{0}=$ 0.55 eV [Fig. 3 (a)]. From these studies, we can see that the calculated $E_{f}-E_{0}$ has an error range about 150 meV, which is comparable to the width of $s$ bands [about 430 meV, see Fig. 3 (f)]. Thus, we argue that it is quite hard to predict whether the Fermi level crosses the $s$ or $p$ bands in such Cu/CO system by DFT simulations. Furthermore, when a thin slab is used in the DFT simulations, there is an artificial shift of surface bands (about 150 meV) resulted from the finite size effect, see Fig. 3 (b). It makes the prediction of $E_{f}$ even worse.

In above, we assume that the number of Cu (111) surface electrons is fixed even if CO molecules are adsorbed. It is reasonable because that we can not expect one hundred of CO molecules can significantly change the Fermi level of Cu bulk. To further illustrate this point, we calculate the CO adsorption induced charge difference per supercell with different CO spacing. The results are plotted in Fig. 10 (a). Here, $z$-direction is perpendicular to the Cu surface, and the position of Cu, C and O atoms are denoted in the figure. We see that the charge difference only occurs at the top two Cu surface layers. The integral of the charge difference along the $z$-direction gives a charge transfer of 0.04 $e^{-}$ per supercell from Cu slab to the CO molecules for all the CO spacing used in the calculation. This value is very tiny, considering the large number of bulk states in each supercell of the Cu slab. Thus, it indicates that the charge transfer here can not influence the absolute position of bulk Fermi level. We also calculate the DOS for various cases, with or without CO adatoms, in Fig. 10 (b). There is nearly no change of DOS no matter whether the CO molecules are absorbed or not, and the Fermi level is always at the same position of DOS. All these calculation results suggest that the influence of charge transfer on Fermi level can be safely ignored.

IV Summary

In summary, we illustrate that the desired $p_{x,y}$-orbital honeycomb electron lattice can be realized on Cu surface by depositing CO molecules into a hexagonal lattice with STM tip. By DFT simulations, we get the modified surface bands and the effective muffin-tin parameters, which gives a direct support to the muffin-tin model. Based on these model parameters, we argue that using a CO cluster instead of a single CO adatom may be a more suitable experimental design to observe the $p$ band features in further experiment. Then, we use the muffin-tin and TB model to illustrate that such Cu/CO systems do give rise to a $p_{x,y}$-orbital honeycomb electron lattice. We further calculate the LDOS patterns of the $p$ bands and then give an analytic interpretation with $k\cdot p$ method. The LDOS of $p$ bands have some unique patterns, which can be used to identify the $p$ bands in further STM experiment. Finally, we point out that the two kinds of edge states, i.e. zero energy and dispersive edge states, can be readily observed in this $p$-orbital honeycomb electron lattice. We expect that such two kinds of edge states can be confirmed in further STM experiment.

A further suggestion about experiment is the lattice constant $a_{0}$, namely the distance between adjacent CO atoms. In order to detect the $p$ bands, the Fermi level should be across the $p$ bands. In the Cu/CO systems, the position of Fermi level depends on $a_{0}$. Our analysis suggests that an appropriate region is $2.5\,\textrm{nm}<a_{0}<4.4\,\textrm{nm}$.

Note added: During the revision of our previous manuscript Qiu et al. (2019), we note that Thomas S. Gardenier and collaborators also study the $p$ orbitals in artificial honeycomb lattice in an STM experiment Gardenier et al. (2020).

Acknowledgement

This work was supported by the National Key Research and Development Program of China (Grant No. 2017YFA0403501), the National Natural Science Foundation of China (Grants No. 11874160, 11534001, 21873033, 11274129, 21873033) and the program for HUST academic frontier youth team. Wen-Xuan Qiu also thanks Li-Lan Qin for the helpful support on the DFT simulations.

Appendix A DFT simulations of LDOS

By DFT simulations, we also calculate the LDOS around the Cu surface (Cu/CO system with a $8\times 8$ supercell). Fig. 11 (a-e) are the LDOS got by DFT simulations with bias voltage in regions [-0.70 V, -0.69 V], [-0.29 V,-0.28 V], [0.01 V, 0.02 V], [0.19 V, 0.20 V] and [0.33 V, 0.34 V], respectively. Note that the energy intervals above all belong to the energy region of the $s$ and $p$ bands in Fig. 3 (f).

Here, we first would like to discuss the relation between the LDOS calculated with the DFT simulations and that with muffin-tin potential model. The muffin-tin potential model deals with a 2D electron system, so that the unit of LDOS ${\rho}_{sf}^{2D}$ got by Eq. (2) is $\textrm{eV}^{-1}\textrm{nm}^{-2}$. As for the DFT simulations, we use a 3D Cu slab to mimic the surface states. The corresponding LDOS is actually the 3D partial charge density at Cu surface in a specified energy interval $\Delta E$, namely, ${\rho}_{sf}^{3D}$, so that the unit is e/$\textrm{\AA}^{3}$. We have a simple relation

[TABLE]

where, $\Delta h$ is the thickness of the surface states and is a constant. When we calculate the LDOS at Cu surface with DFT simulations, we choose very small $\Delta E$, in which the 2D surface states dominate if the slab is thick enough, so ${\rho}_{sf}^{3D}\propto{\rho}_{sf}^{2D}$. In this case, the obtained DFT results indeed reflect the LDOS of surface state electrons and are qualitatively in agreement with that from muffin-tin potential model.

In principle, near the surface, the LDOS in the energy region of surface bands should only result from the Cu surface states. However, in the DFT simulations, if the thickness of the Cu slab is finite, the contribution of the bulk bands can not be completely ignored. To illustrate this point, we project all the bands of a pure Cu slab into the surface layers (top and bottom two), in order to see their contributions to the surface LDOS. Note that, in Fig. 3, we project the bands to a special orbital of top two surface layers so as to identify the surface bands. But here in Fig. 11 we need to project to all the orbitals of surface layers, since that LDOS contains all the surface orbitals. In Fig. 11 (a), we plot the results of a 51-layer-thick Cu slab, where the surface bands in this case are in Fig. 3 (a). We see that, in the energy region of the surface bands, no other bands contribute to the surface layers, except the surface bands. Fig. 11 (b) gives the same projection for a 8-layer-thick Cu slab, where the corresponding surface bands are plotted in Fig. 3 (b). We see that, except for the surface bands, other bulk bands also have distribution in the surface layers. It implies that the surface LDOS got from the DFT simulations contains some contributions from bulk bands, which is different from the real situation. Therefore, compared with the DFT simulations, the muffin-tin model can give a more reliable description about the $p$ band LDOS patterns of the Cu/CO systems.

Fortunately, despite the interference from the bulk bands, the LDOS got by DFT simulations in Fig. 11 (a-e) are still qualitatively in agreement with that from the muffin-tin model, see in Fig. 6 (f-j). When we increase the energy from the bottom of surface band, all the typical LDOS of the $p$ bands in Fig. 3 (f) are found in the DFT simulations.

Finally, we show that, even in the 8-layer-thick Cu slab, the surface bands identified in Fig. 3 are indeed localized at the surface layers. Our way is also the layer resolved band projection. Here, we project all the bands of a Cu slab to the bulk part of the Cu slab. For example, in Fig. 11 (h), we consider the case of a 51-layer-thick Cu slab, where we project the bands into the inner 47 layers. We see that the surface bands here have nearly no contribution to the bulk layers, see the blue box. It indicates that in this case the surface bands are strongly localized at the Cu surface. The situation is similar in the case of 8-layer-thick Cu slab, as shown in Fig. 11 (i).