Front. Phys. Frontiers in Physics Front. Phys. 2296-424X Frontiers Media S.A. 10.3389/fphy.2017.00013 Physics Review Edge-Corrected Mean-Field Hubbard Model: Principle and Applications in 2D Materials Zhang Xi 1 * Wang Tianlei 1 Chen Wencong 1 Wang Sanmei 2 Peng Da 1 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

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

*Correspondence: Xi Zhang zh0005xi@szu.edu.cn

This article was submitted to Physical Chemistry and Chemical Physics, a section of the journal Frontiers in Physics

 16 05 2017 2017 5 13 05 03 2017 25 04 2017 Copyright © 2017 Zhang, Wang, Chen, Wang and Peng. 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.

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.

 Hubbard model tight-binding 2D material nanoribbon BOLS edge effect electronic structure National Natural Science Foundation of China10.13039/501100001809 Natural Science Foundation of Guangdong Province10.13039/501100003453 Introduction

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 [1] and, bare or Hydrogen ended [2]; 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 [11]. Antimonene is stable at high temperature as high as 1,000 K [12] and becomes semiconducting when it is a one atomic layer [13].

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 [17].

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 [27], 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 [33]. Edge effects were also modeled for a specific system, like graphene, from the geometrical perturbation [34], curvature-strain [35, 36], and enhancement of hopping integral at edge [37] 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.

 Principle

The starting point for any discussion of the tight-binding method for electronic and atomic structure calculations must be Slater and Koster [27]. 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 [29]. TB parameters and codes have been developed by various groups, such as Harrison [40], Papaconstantopoulos with NRL-TB code [28], Bowler with DensEL code [41], and Seifert with DFTB code [30, 42].

Table 1 summarized the TB methods of Slater and Koster [27], Harrison [40], Papaconstantopoulos with NRL-TB code [28], Bowler with DensEL code [41], 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.

Comparison of current TB methods with BOLS-HM.

Basis set System Hamiltonian matrix parameter Electron-electron interaction Self-consistent charge
Slater Löwdin orthogonal orbitals Periodic system “Disposable constants” firstly proposed No No
Harrison LCAO; Augmented plane wave for metals Periodic system Fit from experiments No No
NRL-TB code Augmented plane wave Periodic system Fit from DFT Embedded atom method No
DensEL code LCAO Surface structure Fit from DFT No No
DFTB code LCAO Periodic and non-periodic system Fit from DFT Electron repulsion and exchange-correlation Yes
BOLS-HM LCAO Nanomaterial under-coordination system Fit from BOLS correlation Hubbard repulsion term Yes

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.

Conventional TB

TB is firstly based on two assumptions: Born-Oppenheimer approximation [38, 39] and Single electron approximation [27]. Without considering electron-electron interaction, the simplest TB Hamiltonian of one electron in a unit cell is written as:

H0TB(r) = −ℏ22m∇2+V(r) = Hatom(r−Ri) + Vcry(r)

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 [27]:

“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 [43] 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: H0TBck=εkSck, with overlap integral matrix Sivjv′=∫ϕiv*·ϕjv′dr and Hamiltonian matrix element:

Hijvv′TB={ Δivifi=j,v=v′tijvv′if i and j are neighbors0otherwise

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 [44] provided a scheme for creating an orthogonal basis set, Löwdin functions, which are defined as

uv,i = ∑j,v′Sivjv′−1/2ϕjv′

where S is the overlap matrix.

While Bloch functions are periodic and non-localized, Wannier functions [45] construct orthogonal “atomic” wave functions which can take non-periodic terms into account. Wannier function can be expressed as Wannier [45]:

av,i(r−Ri) = N−1/2∑kexp(−ikRi)ψv,k(r)

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 [45].

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 [46] and used in LMTO-TB calculation [47, 48] and Gaussian-type orbital (GTO) proposed by Boys [49]. 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 H0TB 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) [50]:

E=∑v<Ψv(r)|−ℏ22m∇2+V(r)+12∫n(r′)| r−r′ |dr′|Ψv(r)>          + EXC[n(r)]+12∑i,jNZiZj|Ri−Rj|

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 E0DFTB becomes,

E0DFTB=∑v<Ψv(r)|−ℏ22m∇2+V(r)+12∫n0(r′)| r−r′ |dr′                   + VXC[n0]|Ψv(r)>+ 12∑i,jNZiZj|Ri−Rj|              = ∑v<Ψv(r)|H0KS|Ψv(r)>+ Erep

This non-self-consistent-charge (non-SCC) DFTB approach has been successfully applied to various problems in different systems and materials [33]. Moreover, second-order SCC-DFTB [50] 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 [51]. However, DFTB [50] 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 [52], 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 [52]:

H=∑i,σΔici,σ+ci,σ+∑<i,j>,σtij(ci,σ+cj,σ+cj,σ+ci,σ)           + ∑ijkl,σ<ij|1/r|kl>ci,σ+cj,−σck,σ+cl,−σWith<ij|1/r|kl>     =e2∫aσ(r−Ri)*a−σ(r−Rj)aσ(r′−Rk)*a−σ(r′−Rl)r−r′          drdr′

where spin sign σ = ±1/2 and ci,σ+andci,σ 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 H0TB. The last term is the inter-electron potential energy, which is usually a multi-center integral. Hubbard [52] 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:

H=∑i,σΔici,σ+ci,σ + ∑<i,j>,σtij(ci,σ+cj,σ+cj,σ+ci,σ)+∑i,σUni,σni,−σ

where ni,σ=ci,σ+ci,σ 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 [53]. For example, Hubbard-Mott insulators [54] 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 [55]. 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 [56] and those from a tight-binding model are called Tamm states [57]. 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.

Simplified one-dimensional potential model.

The one-dimensional surface wave function at E state can be expressed as:

Ψ(x) = Aexp[ −2m(Vg−E)xℏ ]

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]:

{ cz=di/dB=2/{ 1+exp [ (12−z)/(8z) ] }                                                           (Bond Strain)cz−m=Ei/EB                                      (Bond Strength)zzBcz−m=Eic/EBc                                  (Atomic Cohesive  Energy)

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 [14]. The second equation indicates the single bond energy increases by cz-m as length decreases. The atomic cohesive energy Ec will change as the sum over all the bonds of the atom.

(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 [14]. (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 [14].

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 [14] and listed in Table 2.

Formulation of the measurable quantities on the bonding parameters [14].

Detectable quantity (Q) q(z, m, d, Ei)
Critical temperature (Tc) ∝zEi
Young's modulus (Y) ∝Eid−3
Core level shift (Ev(z)−Ev(0)) ∝Ei
Raman optical shift(ω) ∝zd(Eiμ)1/2

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:

Hijvv′TB={ ∑jeik·(Rj−Ri)∫ϕiv∗(r−Ri) · (Hatom+Vcry)·ϕjv′(r−Rj)dr  (1−,2−D)∫ϕiv∗(r−Ri) ·(Hatom+Vcry) ·  ϕjv′(r−Rj)dr                                    (0−D)

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 Hijvv′TB becomes:

{ Δiv= <iv|Hatom+Vcry|iv> =      εv+α                         (Diagonal Elements)tijvv′={ ∑jeik·(Rj−Ri)< iv|Vcry|jv′>=f(r)·βijvv′                                    (Off - diagonal)<iv|Hatom+Vcry|jv′> = βijvv′                                             (0 - D Off - diagonal)

with,

{ α = <iv|Vcry|iv>                 (Onsite integral)βijvv′ = <iv|Vcry|jv′>       (Hopping integral)

where the diagonal element Δi consists of the eigen-energy εv =< iv|Hatom|iv > of the vth level plus the onsite integral α. βijvv′ is the hopping exchange integral of the nearest neighbor jth atom. f (r) is the periodic factor, ∑jeik·(Rj-Ri). Since the atomic orbitals of a same atom are orthogonal to each other, <iv|Hatom|iv′> =0; also, due to the localization of atomic Hamiltonian, <iv|Hatom|jv′> =0.

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 [58].

 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 βijvv′ into consideration, the modification of under-coordinated site i to the bulk B can be expressed as,

{  Cz(zi) = di(zi)diB=21+exp(12−zi8zi)(Bond Contract) Cz−m = Ei(zi)EiB=V(zi)VB(Bond Energy; Potential)  α(zi)αB = <φi(r)|V(zi)|φi(r)><φi(r)|VB|φi(r)>         ≅V(zi)VB≅Cz−m(Onsite Integral)βij(zi)βijB = <φi(r)|V(zi)|φj(r)><φi(r)|VB|φj(r)>    =Ei(zi)EiB = Cz−m(Hopping Integral)

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 Cz-m; 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 Cz-m 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 [67], 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 [33]. 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.

Correlation between the βij and the bond length d derived by BOLS–TB compared with previous TB fitting from the DFT results [33].

 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 [52] 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:

{ Cz(z)=di(zi)/diB(Bond Contraction)Ui∝di(zi)−3(Volume)Ui(zi)/Uj(zj)=[Cz(zi)/Cz(zj)]−3(Repulsion Energy)

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, [Cz(zi)/Cz(zj)]-3.

 Application of BOLS-HM in 2D materials Electronic structure of graphene nanoribbons Challenges

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 [37] 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 [69]. 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

{ Δi=ε2pz+<2pz|Vcry|2pz>=ε2pz + α                            (Diagonal Elements)tij=fij(x) ·<2pzi|Vcry|2pzj> =fij(x)·βij          (Off−diagonal Elements)

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 fij(x)=∑jexp[ik·(Xi-Xj)]. 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.

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. [71].

 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.

{ βij(zi)= <φi(r)|V−Vi|φj(r)  >∝Ei(zi)(Overlap integral)β1β= Cz(2)−2.56Cz(3)−2.56=1.49(Bare edge integral)β2β= Cz(2.5)−2.56Cz(3)−2.56=1.17(Bare edge-interior integral)β3β= CzH−2.56=1.11(H-edge integral)

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 [67], 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 [68]. Hence, the Hamiltonian as an example, of AGNR is

H=| ABCDEF Δ′β1f1∗β2f2000⋯β1f1Δ′0β2f2*0βf1 β2f2*0Δβf1βf2*0 0β2f2βf1*Δ0βf2 00βf20Δβf1* 000βf2*βf1Δ ⋮     ⋱ |

With

f1=∑Bexp(ik·dBA)=cos(k/3)+isin(k/3)f2=∑Cexp(ik·dAC)=cos(k/23)+isin(k/23)

The first Brillion zone is k∈(π/3)[-1, 1]. The direction of k is along the x axis. 3a is the lattice constant of the unit cell of AGNR.

Considering Hamiltonian matrix, state vector matrix Ck=|cv(k),cv′(k),…> at a specific k point applies the equation:

CkTHCk=Σ

where Σ is a diagonal matrix with diagonal elements of εv,εv′, etc. Ck is an orthogonal matrix with <cv(k)|cv′(k)> = 0. Through diagonalizing the Hamiltonian Matrix, by using QR (orthogonal and upper triangular matrix decomposition) iteration [72], Arnoldi Iteration [73], or Jacobi method [74], 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.

 E<sub>G</sub> 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 [19]. Dangling bonds generates localized energy state near EF (blue line), but does not determine EG [2, 75].

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. [71].

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 [77]. Hydrogenation does not determine the EG opening, but the hopping integral enhancement does.

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. [71].

 Spin and magnetism of GNR Challenges

There was a puzzle about the mechanism of the edge-discriminative generation of the (Dirac-Fermi Polarons) DFPs [78]. 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 [81]. It has been suggested that the interlayer interaction should generate the DFPs [82]. However, the mechanism only applies to graphene of two or more layers [83–85]. Clarifying the principle of the DFPs generating [86] is a burning desire.

 Spin-polarized BOLS-HM parameters of graphene

In the TB Hubbard model, the spin-polarized Hamiltonian can be expressed as:

Hσ=∑i,σΔici,σ+ci,σ+∑<i,j>,σβij(fijci,σ+cj,σ+fij∗cj,σ+ci,σ)            +∑iUin^i,σn^i,−σ

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 [87], the hopping integral tD between the dangling σ–electron orbital and the pz orbital of the nearest neighbors is about −0.5 eV [88]. 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:

gi(E) = 1A∑v,kcivk2exp[ −12(E−εvkσs)2 ]

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.

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 [89].

The structure of vacAGNR-20 (up) and the BOLS-HM derived LDOS of the defective AGNR (down). Reprinted with permission from Zhang et al. [90].

 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.

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. [90].

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 (3d) and leave the dangling σ bond which may provide edge Dirac states.

Hybrid sp2 orbitals of carbon atom are plotted at edges of (A) AGNR and (B) ZGNR.

 E<sub>G</sub> expansion of black phosphorene Challenges

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.

 BOLS-HM scheme

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 [31]:

Δ=Vcry(N)−Vcry(∞)Vcry(∞)=Eg(N)−Eg(∞)Eg(∞)     =∑i≤3γi(EiEb−1)=∑i≤3γi(Ci−m−1)γi=ViV=τCiN−1

γ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 [98]. 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.

BOLS reproduction of the DFT calculated size-dependency of (A) zigzag PNRs and (B) armchair PNRs. Reprint with permission from Liu et al. [98].

 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 [11]. Antimonene is stable at high temperature as high as 1,000 K [12] and becomes semiconducting when it is a one atomic layer [99]. Because there are a large number of dangling bonds and broken bonds in the edge of antimonene nanoribbon (SbNR) [100], 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.

DFT-derived band gap dependence of width of the zigzag SbNRs. The red lines correspond to the BOLS theoretical prediction.

 Summary

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:

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.

Construct the Bloch sums ψv,k(r) of wave functions for low-dimentional nanomaterial.

Construct the Hamiltonian matrix with strong correlation Hubbard model of the nanosystem.

Solve the spin-polarized electronic structure self-consistently to get eigenenergy and eigenvectors Ck=|cv(k),cv′(k),…> for each k point.

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.

 Author contributions

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.

 References Hod O Barone V Peralta JE Scuseria GE. Enhanced half-metallicity in edge-oxidized zigzag graphene nanoribbons. Nano Lett. (2007) 7:2295–9. 10.1021/nl070892217628112 Barone V Hod O Scuseria GE. Electronic structure and stability of semiconducting graphene nanoribbons. Nano Lett. (2006) 6:2748–54. 10.1021/nl061703317163699 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. 10.1103/PhysRevLett.93.18720215525202 Pereira VM Guinea F Lopes dos Santos JMB Peres N Castroneto AH. Disorder induced localized states in graphene. Phys Rev Lett. (2006) 96:036801. 10.1103/PhysRevLett.96.03680116486750 Palacios JJ Fernández-Rossier J Brey L. Vacancy-induced magnetism in graphene and graphene ribbons. Phys Rev B (2008) 77:195428. 10.1103/PhysRevB.77.195428 Cervenka J Katsnelson MI Flipse CFJ. Room-temperature ferromagnetism in graphite driven by two-dimensional networks of point defects. Nat Phys. (2009) 5:840–4. 10.1038/nphys1399 Enoki T Kobayashi Y Fukui KI. Electronic structures of graphene edges and nanographene. Int Rev Phys Chem. (2007) 26:609–45. 10.1080/01442350701611991 Soldano C Mahmood A Dujardin E. Production, properties and potential of graphene. Carbon (2010) 48:2127–50. 10.1016/j.carbon.2010.01.058 Lu SB Miao LL Guo ZN Qi X Zhao CJ Zhang H . Broadband nonlinear optical response in multi-layer black phosphorus: an emerging infrared and mid-infrared optical material. Opt Express (2015) 23:11183–94. 10.1364/OE.23.011183 Guo Z Zhang H Lu S Wang Z Tang S Shao J . From black phosphorus to phosphorene: basic solvent exfoliation, evolution of raman scattering, and applications to ultrafast photonics. Adv Funct Mater. (2016) 25:6996–7002. 10.1002/adfm.201502902 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. 10.1186/s11671-016-1666-427757940 Lei T Liu C Zhao JL Li JM Li YP Wang JO . Electronic structure of antimonene grown on Sb2Te3 (111) and Bi2Te3 substrates. J Appl Phys. (2016) 119:015302. 10.1063/1.4939281 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. 10.1002/anie.20141124625564773 Sun CQ. Size dependence of nanostructures: impact of bond order deficiency. Prog Solid State Chem. (2007) 35:1–159. 10.1016/j.progsolidstchem.2006.03.001 Sun CQ. Oxidation electronics: bond-band-barrier correlation and its applications. Prog Mater Sci. (2003) 48:521–685. 10.1016/S0079-6425(03)00010-0 Sun CQ. Thermo-mechanical behavior of low-dimensional systems: the local bond average approach. Prog Mater Sci. (2009) 54:179–307. 10.1016/j.pmatsci.2008.08.001 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. 10.1021/jp076531s 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. 10.1021/jp807580t Sun CQ Sun Y Ni Y Zhang X Pan J Wang XH . Coulomb repulsion at the nanometer-sized contact: a force driving superhydrophobicity, superfluidity, superlubricity, and supersolidity. J Phys Chem C (2009) 113:20009–19. 10.1021/jp907726b Sun CQ Nie Y Pan J Zhang X Ma SZ Wang Y . 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. 10.1039/c2ra00512c Sun CQ. Dominance of broken bonds and nonbonding electrons at the nanoscale. Nanoscale (2010) 2:1930–61. 10.1039/c0nr00245c20820643 Ghoniem NM Busso EP Kioussis N Huang H. Multiscale modelling of nanomechanics and micromechanics: an overview. Philos Mag. (2003) 83:3475–528. 10.1080/14786430310001607388 Dimitri DV. Multiscale modelling of nanostructures. J Phys Condens Matter (2004). 16:R1537. 10.1088/0953-8984/16/50/R01 Fu CC Torre JD Willaime F Bocquet JL Barbu A. Multiscale modelling of defect kinetics in irradiated iron. Nat Mater. (2005) 4:68–74. 10.1038/nmat1286 Hohenberg P Kohn W. Inhomogeneous electron gas. Phys Rev. (1964) 136:B864. 10.1103/PhysRev.136.B864 Kohn W Sham LJ Self-consistent equations including exchange and correlation effects. Phys Rev. (1965) 140:A1133. 10.1103/PhysRev.140.A1133 Slater JC Koster GF. Simplified LCAO method for the periodic potential problem. Phys Rev. (1954) 94:1498. 10.1103/PhysRev.94.1498 Papaconstantopoulos DA Mehl MJ. The Slater–Koster tight-binding method: a computationally efficient and accurate approach. J Phys Condens Matter (2003) 15:R413. 10.1088/0953-8984/15/10/201 Goringe CM Bowler DR Hernández E. Tight-binding modelling of materials. Rep Prog Phys. (1997) 60:1447. 10.1088/0034-4885/60/12/001 Seifert G Joswig JO. Density-functional tight binding—an approximate density-functional theory method. Wiley Interdiscip Rev Comput Mol Sci. (2012) 2:456–65. 10.1002/wcms.1094 Froyen S Harrison WA. Elementary prediction of linear combination of atomic orbitals matrix elements. Phys Rev B (1979) 20:2420–22. 10.1103/PhysRevB.20.2420 Harrison WA. Overlap interactions and bonding in ionic solids. Phys Rev B (1986) 34:2787–93. 10.1103/PhysRevB.34.27879939976 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. 10.1103/PhysRevB.51.12947 Lammert PE Crespi VH. Geometrical perturbation of graphene electronic structure. Phys Rev B (2000) 61:7308–11. 10.1103/PhysRevB.61.7308 Karasev MV. Graphene as a quantum surface with curvature-strain preserving dynamics. Russ J Math Phys. (2011) 18:64–72. 10.1134/S1061920811010079 Wakabayashi K Okada S Tomita R Fujimoto S Natsume Y. Edge states and flat bands of graphene nanoribbons with edge modification. J Phys Soc Jpn. (2010) 79:034706. 10.1143/JPSJ.79.034706 Wang ZF Li Q Zheng H Ren H Su H Shi QW . Tuning the electronic structure of graphene nanoribbons through chemical edge modification: a theoretical study. Phys Rev B Condens Matter (2007) 75:113406. 10.1103/PhysRevB.75.113406 Yu PY Cardona M. editors. Tight-binding or LCAO approach to the band structure of semiconductors. Fundamentals of Semiconductors, 3rd Edn. Berlin, Heidelberg: Springer (2005), p. 68–71. 26237840 Szabó A Ostlund NS. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory. Mineola, TX; New York, NY: Courier Dover Publications (1989). Harrison WA. Electronic Structure and the Properties of Solids - The Physics of Chemical Bond. San Francisco: CA: Freeman (1980). Bowler D DensEl: A Tight Binding Code. Available online at: http://www.ucl.ac.uk/~ucapdrb/DensEl.html Seifert G The DFTB Website. Available online at: http://www.dftb.org Bloch F. Über die quantenmechanik der elektronen in kristallgittern. Z Phys A Hadrons Nuclei (1929) 52:555–600. Lowdin PO. On the non-orthogonality problem connected with the use of atomic wave functions in the theory of molecules and crystals. J Chem Phys. (1950) 18:365–75. 10.1063/1.1747632 Wannier GH. The structure of electronic excitation levels in insulating crystals. Phys Rev. (1937) 52:191–7. 10.1103/PhysRev.52.191 Slater JC. An augmented plane wave method for the periodic potential problem. Phys Rev. (1953) 92:603–8. 10.1103/PhysRev.92.603 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. 10.1103/PhysRevB.34.52539940356 Xie ZL Dy KS Wu SY. First-principles real-space tight-binding LMTO calculation of electronic structuresfor atomic clusters. Phys Rev B (1997) 55:1748. 10.1103/PhysRevB.55.1748 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. 10.1098/rspa.1950.0036 Elstner M Porezag D Jungnickel G Elsner J Haugk M Frauenheim T . Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties. Phys Rev B (1998) 58:7260–8. 10.1103/PhysRevB.58.7260 Koskinen P Häkkinen H Seifert G Sanna S Th F Moseler M. Density-functional based tight-binding study of small gold clusters. New J Phys. (2006) 8:9. 10.1088/1367-2630/8/1/009 Hubbard J. Electron correlations in narrow energy band. Proc R Soc Lond A Math Phys Sci. (1963) 276:238–57. 10.1098/rspa.1963.0204 Hubbard J. The magnetism of iron. Phys Rev B (1979) 19:2626–36. 10.1103/PhysRevB.19.2626 Mott NF Peierls R. Discussion of the paper by de Boer and Verwey. Proc Phys Soc. (1937) 49:72. 10.1088/0959-5309/49/4S/308 Kuiper P Kruizinga G Ghijsen J Sawatzky GA Verweij H. Character of holes in Li_{x}Ni_{1-x}O and their magnetic behavior. Phys Rev Lett. (1989) 62:221–4. 10.1103/PhysRevLett.62.22110039954 Shockley W. On the surface states associated with a periodic potential. Phys Rev. (1939) 56:317–23. 10.1103/PhysRev.56.317 Tamm I. On the possible bound states of electrons on a crystal surface. Phys Z Soviet Union (1932) 1:733. 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. 10.1038/nmat213218327263 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. 10.1021/jp909549f Sun CQ. Surface and nanosolid core-level shift: Impact of atomic coordination-number imperfection. Phys Rev B (2004) 69:045105. 10.1103/PhysRevB.69.045105 Matsui F Matsushita T Kato Y Hashimoto M Inaji K Guo FZ . 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. 10.1103/PhysRevLett.100.20720118518574 Feibelman PJ. Relaxation of hcp(0001) surfaces: a chemical view. Phys Rev B (1996) 53:13740. 10.1103/PhysRevB.53.137409983126 Falvo MR Clary GJ Taylor RM Chi V Brooks FP Washburn S . Bending and buckling of carbon nanotubes under large strain. Nature (1997) 389:582–4. 10.1038/392829335495 Wong EW Sheehan PE Lieber CM. Nanobeam mechanics: elasticity, strength, and toughness of nanorods and nanotubes. Science (1997) 277:1971–5. 10.1126/science.277.5334.1971 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. Sun CQ Bai HL Tay BK Li S Jiang EY. Dimension, strength, and chemical and thermal stability of a single C-C bond in carbon nanotubes. J Phys Chem B (2003) 107:7544–6. 10.1021/jp035070h Son YW Cohen ML Louie SG. Half-metallic graphene nanoribbons. Nature (2006) 444:347–9. 10.1038/nature0518017108960 Son YW Cohen ML Louie SG. Energy gaps in graphene nanoribbons. Phys Rev Lett. (2006) 97:216803. 10.1103/PhysRevLett.97.21680317155765 Liu L Shen Z. Bandgap engineering of graphene: a density functional theory study. Appl Phys Lett (2009) 95:252104. 10.1063/1.3276068 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. 10.1088/0957-4484/21/6/06520120057033 Zhang X Kuo JL Gu M Bai P Sun CQ. Graphene nanoribbon band-gap expansion: broken-bond-induced edge strain and quantum entrapment. Nanoscale (2010) 2:2160–3. 10.1039/c0nr00273a20697611 Francis JGF. The QR transformation, a unitary analogue to the LR transformation I. Comput J. (1961) 4:265–71. 10.1093/comjnl/4.3.265 Arnoldi WE. The principle of minimized iterations in the solution of the matrix eigenvalue problem. Q Appl Math. (1951) 9:17–29. 10.1090/qam/42792 Rutishauser H. Handbook series linear algebra: the Jacobi method for real symmetric matrices. Numer Math. (1966) 9:1–10. 10.1007/BF02165223 Street RA. Hydrogenated Amorphous Silicon. New York, NY: Cambridge University Press (1991). Yu SS Wen QB Zheng WT Jiang Q. Electronic properties of graphene nanoribbons with armchair-shaped edges. Mol Simul. (2008) 34:1085–90. 10.1080/08927020801958795 Han MY Ozyilmaz B Zhang Y Kim P. Energy band-gap engineering of graphene nanoribbons. Phys Rev Lett. (2007) 98:206805. 10.1103/PhysRevLett.98.20680517677729 Fujita M Wakabayashi K Nakada K Kusakabe K. Peculiar localized state at zigzag graphite edge. J Phys Soc Jpn. (1996) 65:1920. 10.1143/JPSJ.65.1920 Okada S Oshiyama A. Magnetic ordering in hexagonally bonded sheets with first-row elements. Phys Rev Lett. (2001) 87:146803. 10.1103/PhysRevLett.87.14680311580669 Nakada K Fujita M Dresselhaus G Dresselhaus MS. Edge state in graphene ribbons: nanometer size effect and edge shape dependence. Phys Rev B (1996) 54:17954–61. 10.1103/PhysRevB.54.179549985930 Zhang YB Brar VW Girit C Zettl A Crommie MF. Origin of spatial charge inhomogeneity in graphene. Nat Phys. (2009) 5:722–6. 10.1038/nphys1365 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. 10.1103/PhysRevB.75.155438 Ma S Xia JH Srikanth V Sun X Staedler T Jiang X . Magnetism of amorphous carbon nanofibers. Appl Phys Lett. (2009) 95:263105. 10.1063/1.3272940 Xu B Yin J Xia YD Wan XG Jiang K Liu ZG. Electronic and magnetic properties of zigzag graphene nanoribbon with one edge saturated. Appl Phys Lett. (2010) 96:163102. 10.1063/1.3402762 Zanolli Z Charlier JC. Spin transport in carbon nanotubes with magnetic vacancy-defects. Phys Rev B (2010) 81:165406. 10.1103/PhysRevB.81.165406 Koskinen P Malola S Häkkinen H. Evidence for graphene edges beyond zigzag and armchair. Phys Rev B (2009) 80:073401. 10.1103/PhysRevB.80.073401 Zheng F Sasaki KI Saito R Duan W Gu BL. Edge states of zigzag boron nitride nanoribbons. J Phys Soc Jpn. (2009) 78:074713. 10.1143/JPSJ.78.074713 Yazyev OV Helm L. Defect-induced magnetism in graphene. Phys Rev B (2007) 75:125408. 10.1103/PhysRevB.75.125408 Ugeda MM Brihuega I Guinea F Gómez-Rodríguez JM. Missing atom as a source of carbon magnetism. Phys Rev Lett. (2010) 104:096804. 10.1103/PhysRevLett.104.09680420367003 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. 10.1016/j.carbon.2011.04.064 Wang X Jones AM Seyler KL Tran V Jia Y Zhao H . Highly anisotropic and robust excitons in monolayer black phosphorus. Nat Nanotechnol. (2015) 10:517–21. 10.1038/nnano.2015.7125915195 Li L Yu Y Ye GJ Ge Q Ou X Wu H . Black phosphorus field-effect transistors. Nat Nanotechnol. (2014) 9:372–7. 10.1038/nnano.2014.3524584274 Qiao J Kong X Hu ZX Yang F Ji W. Few-layer black phosphorus: emerging 2D semiconductor with high anisotropic carrier mobility and linear dichroism. Nat Commun. (2014) 5. Han X Stewart HM Shevlin SA Catlow CRA Guo ZX. Strain and orientation modulated bandgaps and effective masses of phosphorene nanoribbons. Nano Lett. (2014) 14:4607–14. 10.1021/nl501658d24992160 Ramasubramaniam A Muniz AR. Ab initio. Phys Rev B (2014) 90:085424. 10.1103/PhysRevB.90.085424 Baringhaus J Ruan M Edler F Tejeda A Sicot M Talebibrahimi A . Exceptional ballistic transport in epitaxial graphene nanoribbons. Nature (2013) 506:349–54. 10.1038/nature1295224499819 Yoon Y Guo J. Effect of edge roughness in graphene nanoribbon transistors. Appl Phys Lett. (2007) 91:666. 10.1063/1.2769764 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. 10.1039/C6CP08011A28154850 Huang YL Zhang X Ma ZS Zhou YC Zheng WT Zhou J . Hydrogen-bond relaxation dynamics: resolving mysteries of water ice. Coord Chem Rev. (2015) 285:109–65. 10.1016/j.ccr.2014.10.003 Wang Y Ding Y. Electronic structure and carrier mobilities of arsenene and antimonene nanoribbons: a first-principle study. Nanoscale Res Lett. (2015) 10:955. 10.1186/s11671-015-0955-726058516