Edge-Corrected Mean-Field Hubbard Model: Principle and Applications in 2D Materials
- 1Guangdong Provincial Key Laboratory of Micro/Nano Optomechatronics Engineering, Institute of Nanosurface Science and Engineering, Shenzhen University, Shenzhen, China
- 2Key Laboratory of Low-Dimensional Materials and Application Technologies (Ministry of Education) and School of Materials, Science and Engineering, Xiangtan University, Hunan, China
This work reviews the current progress of tight-binding methods and the recent edge-modified mean-field Hubbard model. Undercoordinated atoms (atoms not fully coordinated) exist at a high rate in nanomaterials with their impact overlooked. A quantum theory was proposed to calculate electronic structure of nanomaterials by incorporating bond order-length-strength (BOLS) correlation to mean-field Hubbard model, i.e., BOLS-HM. Consistency between the BOLS-HM calculation and density functional theory (DFT) calculation on 2D materials verified that (i) bond contractions and potential well depression occur at the edge of graphene, phosphorene, and antimonene nanoribbons; (ii) the physical origin of the band gap opening of graphene, phosphorene, and antimonene nanoribbons lays in the enhancement of edge potentials and hopping integrals due to the shorter and stronger bonds between undercoordinated atoms; (iii) the band gap of 2D material nanoribbons expand as the width decreases due to the increasing under-coordination effects of edges which modulates the conductive behaviors; and (iv) non-bond electrons at the edges and atomic vacancies of 2D material accompanied with the broken bond contribute to the Dirac-Fermi polaron (DFP) with a local magnetic moment.
2D Materials demonstrate extraordinary properties compared with bulk material and attract overwhelming attentions. For example, band gap (EG) occurs of graphene nanoribbons(GNRs) while bulk graphite is a good conductor; EG expands monotonically with the inverse of ribbon width of GNR  and, bare or Hydrogen ended ; unexpected magnetism [3–6] and quantum hall effect [7, 8] were detected at the edge and defect sites of GNRs while bulk graphite is diamagnetism; few-layered phosphorene [9, 10] shows high carrier motilities of 1,000 cm2V−1s−1, and a high current on/off ratio of up to 105 at room temperature; Antimony(Sb), non-hygroscopic, gray metal with a layered structure similar to that of BP . Antimonene is stable at high temperature as high as 1,000 K  and becomes semiconducting when it is a one atomic layer .
Once the size of nanomaterials decreases, the ratio of under-coordinated atoms located at surface, edges, and defects increase compared to the total atomic number. Under-coordination effects induce tunability to the properties of nanomaterials, in contrast with the constant properties in bulky species. For example, the quantities such as the Young's modulus, melting point, dielectric constant, and the extensibility of a solid can change with the solid size [14–16]. The binding energy of core electronic levels also generally shift to energies that are lower (larger in absolute value) than those of the bulk as size decreases or when the atom locates around under-coordinated sites such as surface, edge, and nano-islands .
The involvement of the broken bonds and the non-bonding states make the materials at the nanoscale much complicated and hardly to be understood. Some experimental and theoretical studies have been conducted on size dependent properties of noble metal nanostructures, electronic structure of metal nanoparticles and GNRs, and the water anomalies in literature. However, fundamental progress in theory is still lagging far behind the experimental and theoretical exploitations. Recently, bond order-length-strength (BOLS) correlation derived non-bonding electron polarization theory [16, 18, 19] were proposed to address those problems and obtained great success [20, 21]. Under the light shed by BOLS theory, we will explore the mysteries brought by broken bond and non-bond of nanomaterial.
Multi-scale computational modeling of materials is becoming a reliable tool to underpin scientific investigations and to complement traditional theoretical and experimental approaches [22–24]. At the atomic scale, the ab initio approaches were developed to model the material purely based on Quantum mechanics law without fitting from experimental data, such as Hartree-Forck theory and density functional theory (DFT) [25, 26]. At the large scale above ~1,000 atoms, the classical force field of MD is fit for simulations. Tight-binding (TB) approach, although developed earlier than DFT , is developing fast in recent decades and widely used in the investigation of system of 100–1,000 atoms as reviewed in articles [28–30]. However, under-coordinated modifications at edge, defect and surface sites to conventional TB method are usually neglected in low-dimensional systems or investigated case by case [31, 32].
Edge modification can be presented in the edge Hamiltonian matrix element as a function of distance between edge atoms by fitting distance-dependent TB parameters . Edge effects were also modeled for a specific system, like graphene, from the geometrical perturbation , curvature-strain [35, 36], and enhancement of hopping integral at edge  and so on.
By extending the BOLS mechanics, a quantum theory of BOLS-corrected Hubbard model (BOLS-HM) was proposed with applications for 2D materials graphene, phosphorene, and antimonene nanoribbons of electronic, lattice vibronic, catalytic, and magnetic properties of the material under consideration.
In this themed report, we firstly reviewed the current progresses of TB calculations (conventional TB, density functional TB, spin-polarized Hubbard model, and edge states), and proposed our BOLS-HM model. Then, the applications of BOLS-HM model on the electronic structure and magnetism of GNRs, EG expansion of phosphorene, and electronic properties of antimonene were reviewed and discussed. Finally, we summarized the themed report.
The starting point for any discussion of the tight-binding method for electronic and atomic structure calculations must be Slater and Koster . TB method is developing fast in recent decades and widely used in the investigation of large system as reviewed in articles [28–30] and books [38, 39]. TB method is a parameterized and semi-empirical calculation method which can deal with much larger system with typically two to three orders of magnitude faster than ab initio methods do . TB parameters and codes have been developed by various groups, such as Harrison , Papaconstantopoulos with NRL-TB code , Bowler with DensEL code , and Seifert with DFTB code [30, 42].
Table 1 summarized the TB methods of Slater and Koster , Harrison , Papaconstantopoulos with NRL-TB code , Bowler with DensEL code , and Seifert with DFTB code [30, 42]. The comparison is considered from basis set, applicable system, parameterization of Hamiltonian matrix, electron-electron interaction, and self-consistent charges. It can be seen from the table that in early years (before 1990's) TB parameters were mainly fit from experiments but now are mainly from DFT results. There still lacks a TB method designed to nanomaterial and under-coordination system.
The main advantage of BOLS-HM is that it is designed for nanomaterial and under-coordination systems, since BOLS correlation theory focuses on the bond length, energy, elastic, electronic, optic, dielectric, and other properties of under-coordination system. Besides, BOLS-HM not only improves the Hamiltonian matrix element, it can also predict the physical and chemical properties of nanostructures such as the splitting between the bonding and anti-bonding band of GNR and of valence d band of noble metals, induced by the enhancement of Hamiltonian integrals and by the localized electron-electron repulsions.
TB is firstly based on two assumptions: Born-Oppenheimer approximation [38, 39] and Single electron approximation . Without considering electron-electron interaction, the simplest TB Hamiltonian of one electron in a unit cell is written as:
where the first term is the kinetic energy T of electron and V(r) is the periodic potential from ion-cores, which can be divided into atomic potential Vatom(r) and crystal potential Vcry(r).
As Slater and Koster stated in the paper :
“In fact, if we were going to use n atomic orbitals per unit cell, we could make any n linear combinations of the original orbitals, form Bloch sums of these modified orbitals, and solve a secular problem using the modified Bloch sums, and in every case come out with the same answer in the end. The advantage in one choice of atomic orbitals over another is convenience in calculating the matrix components or solving the secular equation…”
Thus, choosing a proper basis set, although cannot change the final answer, can indeed make the calculation much simpler and of more efficiency. Bloch  provided the formal mechanism for dealing with periodic systems, such as crystals, by means of the Bloch sum in form of LCAO.
Solution of energy states of Equation (1) is transferred to an eigenvalue problem of: , with overlap integral matrix and Hamiltonian matrix element:
where, tijvv′ is the off-diagonal element between vth orbital at ith atom and v′th orbital at the equivalent nearest jth atoms. However, the atomic orbitals, ϕv,i(r − Ri), are not ideal for the purposes of analysis, as the orbitals on different atomic sites are not orthogonal to one another.
Löwdin  provided a scheme for creating an orthogonal basis set, Löwdin functions, which are defined as
where S is the overlap matrix.
While Bloch functions are periodic and non-localized, Wannier functions  construct orthogonal “atomic” wave functions which can take non-periodic terms into account. Wannier function can be expressed as Wannier :
where ψv,k(r) is the Bloch sum.
Unlike other basis set which are expanded in k-space, Wannier function is a real-space basis set and can be a powerful tool in the study of the electronic and dielectric properties of materials .
Besides the above, many other excellent basis sets have also been proposed and widely applied in the electronic structure modeling of materials, such as linear muffin-tin orbital (LMTO) proposed by Slater  and used in LMTO-TB calculation [47, 48] and Gaussian-type orbital (GTO) proposed by Boys . To choose a proper basis set for a system is an important part in saving time and improving the efficiency and accuracy in the electronic structure calculations.
Density-Functional TB Method
As discussed, the single-electron Hamiltonian only considers the potential energy from ion cores in crystal and neglects the interaction among electrons.
The total energy of a system of M electrons in the field of N nuclei at positions Ri may be written within DFT as a functional of a charge density n(r) :
where the first sum includes kinetic energy T and potential energy V(r) from crystal ion-cores and potential energy from other electrons Ve, the second term is the exchange-correlation (XC) contribution, and the last term covers the ion-ion core repulsion, Erep.
For 0th order density-functional tight-binding (DFTB), charge-density n(r) which should be solved by self-consistent field (SCF) method in DFT is now approximated by the input charge-density n0(r). And the total energy becomes,
This non-self-consistent-charge (non-SCC) DFTB approach has been successfully applied to various problems in different systems and materials . Moreover, second-order SCC-DFTB  has also been included n(r) = n0(r) + δn(r) in the corresponding Hamiltonian terms in Equation (6), and also been widely used and accepted . However, DFTB  is purely parameterized by DFT results without considering experimental results and complicated to be realized.
Spin-Polarized Hubbard Model
In incompletely filled electron shells with narrow energy bands, the correlations between electrons are too strong to be neglected as done in conventional TB. The Coulomb repulsion between electrons at the same atomic site will change the band structure significantly. The Hubbard model, named after Hubbard , includes the onsite repulsion which stems from the Coulomb repulsion between electrons at each atomic site in Hamiltonian. Using the real-space Wannier functions av,i(r − Ri) as basis set, Hamiltonian is expressed in second quantization form as :
where spin sign σ = ±1/2 and are the creation and annihilation operator in Wannier representation, indicating to create and to annihilate an electron at the lattice vector Ri. The first and second terms are the same as the . The last term is the inter-electron potential energy, which is usually a multi-center integral. Hubbard  indicated that the single center integral U =< ii, σ|1/r|ii,−σ > is about 10 eV, much greater than other two-center, three-center integrals. Thus, the simplified Hamiltonian becomes:
where represents the particle-number operator at site Ri and spin σ. Equation (8) is called Hubbard model with last term Uni↑ni↓.
Hubbard model plays a significant role in the conductor-insulator transition of metal and magnetism of narrow band . For example, Hubbard-Mott insulators  are a class of materials that should conduct electricity under conventional band theories, but are insulators when measured (particularly at low temperatures). This effect is due to electron-electron interactions in narrow bands which are not considered in conventional band theory. The Coulomb repulsions of semi-localized electrons in narrow bands are strong enough to split the band into two subbands and generate the Mott band gap . This situation is very similar to nanomaterials with limited quantity of electrons and semi-localized electrons at edge and surfaces.
Shockley and Tamm Edge States
At the surface, the termination of a crystal obviously causes deviation from perfect periodicity. Surface states that are calculated in the framework of near-free electron approximation are called Shockley states  and those from a tight-binding model are called Tamm states . However, there is no real physical distinction between the two terms, only the mathematical approach in describing surface states is different.
A simplified model of the crystal potential in one dimension is used: the periodic crystal potential jumps abruptly to the vacuum level Vg at surface, as shown in Figure 1.
The one-dimensional surface wave function at E state can be expressed as:
The Shockley states and the Tamm states are suitable to describe periodical systems. However, the simplification of the stepped potential prevents the method applied in nanomaterials or imperfect surface where the change of multi-well potential at surface should be considered.
BOLS-Corrected Hubbard Model
The BOLS is an extension of the “atomic coordination number (CN)—atomic size” correlation mechanism of Goldschmidt, Pauling, and Feibelman to energy domain. The shorter and stronger bonds between undercoordinated atoms provide significant modification to atomic cohesive energy (Ec) and associated properties when the fraction of the undercoordinated atoms is increased. The consequence of the broken bond follows the BOLS correlation [14, 16]:
where m is the bond nature indicator; subscript i and B denote the ith under-coordinated atom and bulk value, respectively. zzB = z/zB is the reduce coordination with zB = 12 being the bulk standard. First equation indicates that: as the CN z decreases, bond length (di) becomes shorter compared with bulk value (dB); the curve of z-d dependence is illustrated in Figure 2A, compared with measurements of Au particles, carbon nanotubes, Pt, Ir, Ti, Zr, and Zn chains . The second equation indicates the single bond energy increases by as length decreases. The atomic cohesive energy Ec will change as the sum over all the bonds of the atom.
Figure 2. (A) The BOLS correlation mechanism (solid line) formulated from the atomic “CN-radius.” The formulation has been further verified by the measurement (scattered symbols) from Au particles, carbon nanotubes, Pt, Ir, Ti, Zr, and Zn chains . (B) Schematic illustration of the broken-bond induced potential trapping at the terminating edges up to two inter-atomic layers [58, 59]. Reprinted with permission of Sun .
Figure 2B illustrates the BOLS correlation and the quantum trapping of potential and energy [14, 16]. Instead of single potential well of Quantum confinement, multi-well quantum trapping happens to nanostructures. At the edge and surface of nanostructures (particle, rod, chain, etc.), the remaining bonds of the under-coordinated atoms to contract spontaneously with an association of bond strength gain up to two inter-atomic layers. The localized strain in turn causes potential well depression with a consequence of localized densification of charge, energy and mass, which is called as Quantum trapping of the charge and energy.
Any detectable quantity localized nearby the broken bonds can be formulated by the above parameters of bond strain (d), bond energy (E), and atomic cohesive energy (EC) change as shown in Equation (10). The BOLS theory has enabled unification of the unusual performance of detectable quantities Q(z), such as mechanical, thermal, acoustic, chemical, electronic, dielectric, ferroelectric, optic, and magnetic properties and the transport dynamics of electrons, phonons, and photons at the skins of nanostructures of various shapes, as summarized in Sun  and listed in Table 2.
Table 2. Formulation of the measurable quantities on the bonding parameters .
Q of a nanostructure can be expressed as a function q(m, z, d, EB), as sampled in Table 2. Properties of a nanosolid with size (K) depends on shape (τ), and bond nature (m) [16, 60]. BOLS correlation has already been verified of bond identities in the metal nanoparticles, metal surfaces, carbon nanotubes, graphene nanoribbons, etc. [61, 62].
Applying variational principle to the TB Hamiltonian, Hamiltonian matrix element of secular equations becomes:
where Vcry is the ion-cores' potential subtracting the ith atomic potential.
Writing ϕiv(r − Ri) as |iv >, as far as the v and v′ orbitals are concerned, the Hamiltonian matrix element becomes:
where the diagonal element Δi consists of the eigen-energy εv = < iv|Hatom|iv > of the vth level plus the onsite integral α. is the hopping exchange integral of the nearest neighbor jth atom. f (r) is the periodic factor, . Since the atomic orbitals of a same atom are orthogonal to each other, ; also, due to the localization of atomic Hamiltonian, .
CN Imperfection Induced Quantum Trapping
According to the BOLS correlation theory, CN imperfection at edges of low-dimensional nanostructures or at surface skin of 0-D nanoparticles causes the remaining bonds of the under-coordinated atoms to contract spontaneously with an association of bond strength gain, which in turn causes potential well depression with a consequence of localized densification of charge, energy, and mass. BOLS considered that the shorter and stronger bonds between undercoordinated atoms provide significant modification to the Hamiltonian matrix elements.
Since BOLS considered the effective CN of three layers from surface or edges, the modification to potential multi-well at three-layer edges is adopted. As illustrated in Figure 2A, the broken-bond induces potential trapping at the terminating edges up to three atomic layers. Compared with the Tamm model approximating original periodic potential jumping to vacuum level, as shown in Figure 2B, BOLS consideration of the quantum entrapment of potential multi-well at “surface skin” better describes the real nanomaterial observed experimentally .
Enhancement of Onsite and Hopping Integrals
In the present BOLS-HM, we took the following relations of effective coordination number (z), bond length (d), single bond energy (E), potential (V), and Hamiltonian integrals α an into consideration, the modification of under-coordinated site i to the bulk B can be expressed as,
Cz is the bond contraction coefficient and m is a constant of bond nature indicator. As the effective CN becomes imperfect at the edge, the edge bond will be shortened by Cz and strengthened by ; and the potential well of the neighbor atom will become closer and deepened in energy space proportional to bond energy. BOLS-HM considers the relative ratio of a quantity at under-coordinated atom and in bulk. Since the main contribution of Vcry to the ratio of integral at ith atomic site is the atomic potential from the neighbor jth atom Vij, the α(zi) and βij(zi) are both considered as proportional to the bond energy and potential well depression compared to those of bulk.
Determining the bond nature indicator m is accomplished by the consistency between predictions of BOLS correlation theory and experimental results on the various quantities as shown in Table 2. For example, based on the Young's modulus [63, 64], the melting temperature [65, 66] of carbon nanotube and the energy shift of C1s level of grapheme , we obtained the bond nature parameter m = 2.56 of carbon. Figure 3 compares the relation between hopping integrals and the C-C bond length derived by BOLS-HM and obtained by TB fitting from DFT results . It should be noted that DFT treats low dimensional systems with vacuum slabs; while BOLS-HM develops the relations from the consistency between BOLS correlation theory and the experimental results.
Figure 3. Correlation between the βij and the bond length d derived by BOLS–TB compared with previous TB fitting from the DFT results .
Onsite Repulsion Effects
Since limited amount of electrons exist in nanomaterial compared with tremendous amount of electrons in bulk, the electron-electron repulsion at each atom, especially at edge sites, will induce considerable changes to electronic structure. As discussed, the Hubbard model  includes the onsite repulsion which stems from the Coulomb repulsion between electrons at each atomic site in Hamiltonian by introducing the Hubbard term Uni↑ni↓.
BOLS-HM considers the relation of di and Ui as:
where ith atom locates at edge sites and jth atom at the interior with coordination number of zj. Hubbard U is the onsite repulsion energy which should contain the information of the relative ratio of the electron density, .
Application of BOLS-HM in 2D Materials
Electronic Structure of Graphene Nanoribbons
The principle for the EG expansion with width is not clear, apart from the EG disparity between experiments and calculations. Either the carrier confinement [68–70] or the edge energy pinning couldn't be changed by EG expansion and it also couldn't be observed  when calculations for the edge and for the interior were identical which were computed by hopping integrals employed in band structure. At the edge, the hopping integral is extremely greater comparing with the integral in the GNR interior . It's obvious that the EG is determined intrinsically by the crystal potential and hence the Hamiltonian matrix, but the density and energy of carriers in GNR is of a great importance to the transport dynamics though from our proposition.
Hamiltonian of GNRs
As far as the π and π* orbitals of the C 2pz electrons are concerned, the Hamiltonian matrix elements are
where the onsite energy Δi is the sum of the exchange integral and the eigen-energy of the pz electron, ε2pz. V is the crystal potential and Vi is the ith intra-atomic potential; α is the onsite exchange integral of negative value; βij is the hopping integral and . Considered structures of armchair-edge GNR (AGNR) and reconstructed zigzag-edge GNR (RecGNR) are shown in Figure 4. Since periodic direction is along x axis, wave vector k is parallel to x.
Figure 4. Illustration of structure of (A) AGNR and (B) RecGNR edges. The reconstructed edge consists of 7- and 5-atom rings. Periodic direction is along x axis and ribbon width in y axis. The hopping integrals are indicated as β, β 1, and β 2. N indicates the counting of ribbon width; the numbers also mark the atomic position. Reprinted with permission from Zhang et al. .
BOLS-HM Parameter for Graphene
The following relations of effective coordination number (z), bond coefficient (Cz), single bond energy (E), potential (V), and hopping integral were be taken into consideration in the present BOLS-HM.
Abbreviation of the bond contraction coefficient is Cz. CzH corresponds to the C-C bond length coefficient terminated with H. Subscripts B corresponds to bulk. From the measured melting point [65, 66] of carbon nanotubes, elastic modulus [63, 64], and the C1s electron binding energy shift , we have obtained m = 2.56.
In conducting the edge-modified TB calculations, the βij are enlarged by C(z)−2.56 for under-coordinated atoms, as shown in Figure 4. The β1 and β2 (reduced from βij, β for the bulk) at the edge are proportional to the bond energy. Therefore, β1 > β2 > β = 2.4 eV. The H-end edge bond contraction coefficient, CH, is defined as the bond length ratio between GNR edge and interior, and taken as 96% according to DFT relaxation results . Hence, the Hamiltonian as an example, of AGNR is
The first Brillion zone is . The direction of k is along the x axis. is the lattice constant of the unit cell of AGNR.
Considering Hamiltonian matrix, state vector matrix at a specific k point applies the equation:
where Σ is a diagonal matrix with diagonal elements of etc. Ck is an orthogonal matrix with . Through diagonalizing the Hamiltonian Matrix, by using QR (orthogonal and upper triangular matrix decomposition) iteration , Arnoldi Iteration , or Jacobi method , the eigen-energies and eigen-vectors can be obtained at each k point. High symmetry k-points are usually calculated first such as Γ point at the center of Brillouin zone and M point at center of an edge of simple cubic.
EG Opening GNRs
In Figure 5, band structures of bare and H-end RecGNR-8 are compared clearly. Both BOLS-HM and DFT calculation results are plot. By the method of BOLS-HM and DFT, a very small EG (~0.1 eV) is generated near Γ point. EG indeed is dominated by the quantum entrapment of electrons at edges. The bond which is dangling results in an additional state at the Fermi level. The electron entrapment could pin and polarize the non-bonding electron locally . Dangling bonds generates localized energy state near EF (blue line), but does not determine EG [2, 75].
Figure 5. The energy dispersion of bare RecGNR-8 calculated by (A) TB; (B) BOLS-HM; (C) DFT; The energy dispersion of H-end RecGNR-8 calculated by (D) TB; (E) BOLS-HM; and (F) DFT. Reprinted with permission from Zhang et al. .
Figure 6 compares size-dependent EG (N) calculated by conventional TB and by BOLS-HM scheme. AGNR with or without hydrogenation were both considered by BOLS-HM. EG is always zero when N = 3p+2 (p is a natural number) calculated by TB. EG opens with or without hydrogenation calculated by BOLS-HM, accord with the EG observations of DFT calculation [68, 76] and STM results . Hydrogenation does not determine the EG opening, but the hopping integral enhancement does.
Figure 6. EG vs. ribbon width (N) of the in bare and H-end AGNRs calculated by TB and BOLS-HM. Reprinted with permission from Zhang et al. .
Spin and Magnetism of GNR
There was a puzzle about the mechanism of the edge-discriminative generation of the (Dirac-Fermi Polarons) DFPs . The DFPs have been observed as a sharp density of states at Fermi energy (EF) by scanning tunneling microscopy/spectroscopy [79, 80]. However, it is not likely to form DFPs at a clean graphite surface, the interior of GNR, or the armchair edge. Below the surface of graphene, the charge density fluctuations attribute to charge donations of impurities . It has been suggested that the interlayer interaction should generate the DFPs . However, the mechanism only applies to graphene of two or more layers [83–85]. Clarifying the principle of the DFPs generating  is a burning desire.
Spin-Polarized BOLS-HM Parameters of Graphene
In the TB Hubbard model, the spin-polarized Hamiltonian can be expressed as:
Where σ is the spin with sign. fij is the periodic factor. Subscript D denotes dangling σ-bond. The last term is the Hubbard electron repulsion.
We establish the Hamiltonian matrix for atomic vacancy, zigzag, armchair, reconstructed zigzag in the BOLS-HM calculations. Owing to the truth that orbitals have different symmetries in the z-direction , the hopping integral tD between the dangling σ–electron orbital and the pz orbital of the nearest neighbors is about −0.5 eV . Spin-polarized electronic structures are solved by self-consistent field scheme.
Local Spin Density of States
Considering spin-polarized Hamiltonian, spin-polarized electronic structures can be obtained by diagonalizing the Hamiltonian matrix using the mean-field method as described. Local spin density of states (LSDOS) of spin σ at ith atom, giσ(E), can be calculated from the spin-polarized electronic structure.
LDOS at ith atom [gi(E)] is calculated according to the atomic contribution (civk) to the LCAO eigen-vector cvk and expressed as:
Figure 7 illustrates the calculation results of LDOS of AGNR-20. Density distribution at each atomic position and at difference energy can be obtained. Arrows in left panel indicates the quantum trapping of bonding electrons and polarization of anti-bonding states at edge sites obtained by BOLS-HM.
Figure 7. LDOS contour of AGNR-20 calculated by BOLS-HM (left) compared with conventional TB (right). C1 and C2 are the edge sites.
The BOLS-HM derived LDOS of a hydrogenated vacAGNR (H-vacAGNR) with ribbon width N = 20 is shown in Figure 8. A sharp resonant peak at vacancy site near EF can be observed, which is consistent with that identified experimentally from graphite surface vacancy .
Figure 8. The structure of vacAGNR-20 (up) and the BOLS-HM derived LDOS of the defective AGNR (down). Reprinted with permission from Zhang et al. .
Magnetism of Dirac-Fermi Polarons
As Figure 9 shows, the local spin density of states (LSDOS) with and without hydrogenation are compared intuitively, as well as the spin electron density in real space. Zigzag edges and atomic defects demonstrates a spin-density as high as 0.7~1.0 electron/Å3 of dangling σ-bond. Surprisingly, after hydrogenation, zigzag edge still has a low magnetism, due to the contribution of pz electrons trapped at edges. Meanwhile, we found armchair edge exhibiting an extremely low spin-density.
Figure 9. DFT-derived spin electron density and BOLS-HM derived difference LSDOS for (A) ZGNR and (B) H-ZGNR. Zigzag edge manifest a strong local magnetism caused by both Dz and trapped edge pz electron. Hydrogenation annihilates Dz but leaves a weak magnetism. Reprinted with permission from Zhang et al. .
The sp2 hybridization of GNR edges are shown in Figure 10 to help to understanding the results of magnetism of DFPs. At armchair edge in Figure 10A, each edge atom loses one of its neighbors and leaves a dangling electron. Due to the short distance (d) between two atoms along the edge, triple bonds are formed between two atoms and the sp2 hybrid may also change. At zigzag edge in Figure 10B, atoms along the edge have a longer distance () and leave the dangling σ bond which may provide edge Dirac states.
EG Expansion of Black Phosphorene
Few-layer phosphorene [9, 10] shows high carrier motilities, up to 1,000 cm2V−1s−1, and a high current on/off ratio of up to 105 at room temperature. These values make this material a promising candidate for field-effect transistors [91–93]. Modulating the band gap of phosphorene nanoribbon (PNR) [94, 95] is crucial for the semiconducting applications. Because of its anisotropy, the electronic properties are affected considerably by the orientation of the PNRs [96, 97], leading to promising applications. However, the physical origin behind the EG modulation is still unclear.
According to BOLS-HM, the perturbation to the crystal potential can be expressed as follows Vcry (N) = Vcry (∞) (1+ΔN). Vcry (N) and Vcry (∞) denote the crystal potential of nanomaterial size (N) and the bulk, respectively. ΔN is a size-dependent perturbation to the crystal potential. The relationship can also be described as Froyen and Harrison :
γi is the surface-to-volume ratio, and τ is the dimensionality with τ = 1, 2, 3 for a nanoparticle, nanoribbon, or nanofilm, respectively. The bond index m is intrinsic to a specific material. The expressions indicate the size-dependent Vcry and EG are a function of shape and surface-to-volume ratio of a nanomaterial.
Figure 11 compares the DFT-derived EG vs. the ribbon width(N) of zigzag and armchair PNRs with theoretical prediction of BOLS-HM . The m of the PNRs in Equation (20) was optimized to be 4.60 based on the measured data. Using the optimized m-values, we calculated, plotted, and compared the theoretical curves. The necessary values of the reference bandgap energy EbZ = 0.82 eV and EbA = 0.43 eV were obtained. They are very close to the bulk value of two typical edge materials. A further refinement of the derived m, EbZ, and EbA values was achieved by carefully matching the DFT calculations to the whole nanoribbon length ranges.
Figure 11. BOLS reproduction of the DFT calculated size-dependency of (A) zigzag PNRs and (B) armchair PNRs. Reprint with permission from Liu et al. .
Electronic Properties of Antimonene
Encouragingly, the 2D material has been found in group V elements. Antimony(Sb), non-hygroscopic, gray metal with a layered structure similar to that of BP . Antimonene is stable at high temperature as high as 1,000 K  and becomes semiconducting when it is a one atomic layer . Because there are a large number of dangling bonds and broken bonds in the edge of antimonene nanoribbon (SbNR) , the bond parameters of the surface atoms are different from the internal atoms. Like other 2D material, the EG modulation is crucial for its potential applications while the physical origin is not conclusive yet.
Though antimony belongs to rhombohedral system, hexagonal coordinate system is convenient to describe its structure. We identify the SbNR structures by the number of N along the ribbon orientations. Vacuum slabs of 10 Å were inserted in the width direction and the z axis. DFT calculation was performed by CASTEP using PBE functional. Figure 12 compares the DFT-derived EG vs. the ribbon width(N) of zigzag and armchair PNRs with theoretical prediction of BOLS-HM. The reproduction of these quantities confirms the importance of the size of nanoribbons, which supports the proposals that the origin of these novel properties is mainly due to the change of the bond length and strength in different sizes of SbNRs.
Figure 12. DFT-derived band gap dependence of width of the zigzag SbNRs. The red lines correspond to the BOLS theoretical prediction.
In this work, a quantum theory was proposed to calculate the under-coordinated effects on the electronic structure of materials by incorporating BOLS correlation theory to mean-field Hubbard model. The whole process of BOLS-HM dealing with the chemical bonds can be summarized as:
1. Input the effective coordination number z, the bond length at the under-coordinated sites to calculate the increase of bond energy, and the mortification of the inter-atomic matrix elements according to BOLS correlation.
2. Construct the Bloch sums ψv,k(r) of wave functions for low-dimentional nanomaterial.
3. Construct the Hamiltonian matrix with strong correlation Hubbard model of the nanosystem.
4. Solve the spin-polarized electronic structure self-consistently to get eigenenergy and eigenvectors for each k point.
5. Based on the eigenenergy and eigenvector results, properties such as Mulliken charge, Mayer bond order, local density of states can thus be obtained.
By virtue of BOLS-HM, we clarified that (i) bond contractions and potential well depression occur at the edge of graphene, phosphorene, and antimonene nanoribbons; (ii) the physical origin of the band gap opening of graphene, phosphorene, and antimonene nanoribbons lays in the enhancement of edge potentials and hopping integrals due to the shorter and stronger bonds between undercoordinated atoms; (iii) the band gap of 2D material nanoribbons expand as the width decreases which modulates the conductive behaviors; and (iv) non-bond electrons at the edges and atomic vacancies of 2D material accompanied with the broken bond contribute to the DFP with a local magnetic moment. Although, this work shows the BOLS-HM applications on three kinds of 2D materials, the potential applications are not limited to those materials but applied on other species such as Mo/WS2.
XZ coordinates the project. XZ and TW proposed the theoretical model. TW, WC, and DP wrote the manuscript. SW did the DFT calculations.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Financial supports from National Natural Science Foundation (No. 51605306) of China, Guangdong Innovation Youth Fund (No. 2015KQNCX144), and Natural Science Foundation of Guangdong (No. 2016A030310060), and of SZU (No. 827000131), Shenzhen foundation fund (No. JCYJ20160427105015701) are gratefully acknowledged.
3. Lehtinen PO, Foster AS, Ma Y, Krasheninnikov AV, Nieminen RM. Irradiation-induced magnetism in graphite: a density functional study. Phys Rev Lett. (2004) 93:187202. doi: 10.1103/PhysRevLett.93.187202
9. Lu SB, Miao LL, Guo ZN, Qi X, Zhao CJ, Zhang H, et al. Broadband nonlinear optical response in multi-layer black phosphorus: an emerging infrared and mid-infrared optical material. Opt Express (2015) 23:11183–94. doi: 10.1364/OE.23.011183
10. Guo Z, Zhang H, Lu S, Wang Z, Tang S, Shao J, et al. From black phosphorus to phosphorene: basic solvent exfoliation, evolution of raman scattering, and applications to ultrafast photonics. Adv Funct Mater. (2016) 25:6996–7002. doi: 10.1002/adfm.201502902
11. Cheung CH, Fuh HR, Hsu MC, Lin YC, Chang CR. Spin orbit coupling gap and indirect gap in strain-tuned topological insulator-antimonene. Nanoscale Res. Lett. (2016) 11:459. doi: 10.1186/s11671-016-1666-4
13. Zhang S, Yan Z, Li Y, Chen Z, Zeng H. Atomically thin arsenene and antimonene: semimetal-semiconductor and indirect-direct band-gap transitions. Angew Chem. (2015) 54:3112–5. doi: 10.1002/anie.201411246
17. Yang DQ, Sacher E. Strongly enhanced interaction between evaporated Pt nanoparticles and functionalized multiwalled carbon nanotubes via plasma surface modifications: effects of physical and chemical defects. J Phys Chem C (2008) 112:4075–82. doi: 10.1021/jp076531s
18. Sun CQ, Fu SY, Nie YG. Dominance of broken bonds and unpaired nonbonding pi-electrons in the band gap expansion and edge states generation in graphene nanoribbons. J Phys Chem C (2008) 112:18927–34. doi: 10.1021/jp807580t
19. Sun CQ, Sun Y, Ni Y, Zhang X, Pan J, Wang XH, et al. Coulomb repulsion at the nanometer-sized contact: a force driving superhydrophobicity, superfluidity, superlubricity, and supersolidity. J Phys Chem C (2009) 113:20009–19. doi: 10.1021/jp907726b
20. Sun CQ, Nie Y, Pan J, Zhang X, Ma SZ, Wang Y, et al. Zone-selective photoelectronic measurements of the local bonding and electronic dynamics associated with the monolayer skin and point defects of graphite. RSC Adv. (2012) 2:2377–83. doi: 10.1039/c2ra00512c
33. Porezag D, Frauenheim T, Köhler T, Seifert G, Kaschner R. Construction of tight-binding-like potentials on the basis of density-functional theory: application to carbon. Phys Rev B (1995) 51:12947–57. doi: 10.1103/PhysRevB.51.12947
37. Wang ZF, Li Q, Zheng H, Ren H, Su H, Shi QW, et al. Tuning the electronic structure of graphene nanoribbons through chemical edge modification: a theoretical study. Phys Rev B Condens Matter (2007) 75:113406. doi: 10.1103/PhysRevB.75.113406
41. Bowler D. DensEl: A. Tight Binding Code. Available online at: http://www.ucl.ac.uk/~ucapdrb/DensEl.html
42. Seifert G. The DFTB Website. Available online at: http://www.dftb.org
47. Andersen OK, Pawlowska Z, Jepsen O. Illustration of the linear-muffin-tin-orbital tight-binding representation: compact orbitals and charge density in Si. Phys Rev B (1986) 34:5253. doi: 10.1103/PhysRevB.34.5253
49. Boys SF. Electronic wave functions. I. A general method of calculation for the stationary states of any molecular system. Proc R Soc Lond A Math Phys Sci. (1950) 200:542–54. doi: 10.1098/rspa.1950.0036
50. Elstner M, Porezag D, Jungnickel G, Elsner J, Haugk M, Frauenheim T, et al. Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties. Phys Rev B (1998) 58:7260–8. doi: 10.1103/PhysRevB.58.7260
58. Huang WJ, Sun R, Tao J, Menard LD, Nuzzo RG, Zuo JM. Coordination-dependent surface atomic contraction in nanocrystals revealed by coherent diffraction. Nat Mater. (2008) 7:308–13. doi: 10.1038/nmat2132
59. Wang Y, Nie YG, Wang LL, Sun CQ. Atomic-layer- and crystal-orientation-resolved 3d(5/2) binding energy shift of Ru(0001) and Ru(1010) surfaces. J Phys Chem C (2010) 114:1226–30. doi: 10.1021/jp909549f
61. Matsui F, Matsushita T, Kato Y, Hashimoto M, Inaji K, Guo FZ, et al. Atomic-layer resolved magnetic and electronic structure analysis of Ni thin film on a Cu(001) surface by diffraction spectroscopy. Phys Rev Lett. (2008) 100:207201. doi: 10.1103/PhysRevLett.100.207201
65. An B, Fukuyama S, Yokogawa K, Yoshimura M. Surface superstructure of carbon nanotubes on highly oriented pyrolytic graphite annealed at elevated temperatures. Jpn J Appl Phys. (1998) 37(Pt 1):3809–11.
70. Zhang YH, Zhou KG, Xie KF, Zeng J, Zhang HL, Peng Y. Tuning the electronic structure and transport properties of graphene by noncovalent functionalization: effects of organic donor, acceptor and metal atoms. Nanotechnology (2010) 21:065201. doi: 10.1088/0957-4484/21/6/065201
82. Ferro Y, Allouche A. Interpretation of STM images of graphite with an atomic vacancy via density-functional calculations of electronic structure. Phys Rev B (2007) 75:155438. doi: 10.1103/PhysRevB.75.155438
90. Zhang X, Nie Y, Zheng W, Kuo JL, Sun CQ. Discriminative generation and hydrogen modulation of the Dirac-Fermi polarons at graphene edges and atomic vacancies. Carbon (2011) 49:3615–21. doi: 10.1016/j.carbon.2011.04.064
98. Liu Y, Bo M, Yang X, Zhang P, Sun CQ, Huang Y. Size modulation electronic and optical properties of phosphorene nanoribbons: DFT-BOLS approximation. Phys Chem Chem Phys. (2017) 19:5304–9. doi: 10.1039/C6CP08011A
Keywords: Hubbard model, tight-binding, 2D material, nanoribbon, BOLS, edge effect, electronic structure
Citation: Zhang X, Wang T, Chen W, Wang S and Peng D (2017) Edge-Corrected Mean-Field Hubbard Model: Principle and Applications in 2D Materials. Front. Phys. 5:13. doi: 10.3389/fphy.2017.00013
Received: 05 March 2017; Accepted: 25 April 2017;
Published: 16 May 2017.
Edited by:Jiabao Yi, University of New South Wales, Australia
Reviewed by:Zhimin Ao, Guangdong University of Technology, China
Zhaoqiang Bai, Western Digital Technologies, United States
Copyright © 2017 Zhang, Wang, Chen, Wang and Peng. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xi Zhang, email@example.com