ELECTRONIC PROPERTIES OF COMPLEX NANOSTRUCTURES By Zhen Zhu A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Physics-Doctor of Philosophy 2015 ABSTRACT ELECTRONIC PROPERTIES OF COMPLEX NANOSTRUCTURES By Zhen Zhu Nanostructured materials have brought an unprecedented opportunity for advancement in many fields of human endeavor and in applications. Nanostructures are a new research field which may revolutionize people's everyday life. In the Thesis, I have used theoretical methods including density functional theory (DFT), molecular dynamic simulations (MD) and tight-binding methods to explore the structural, mechanical and electronic properties of various nanomaterials. In all this, I also paid attention to potential applications of these findings. First, I will briefly introduce the scientific background of this Thesis, including the motivation for the study of a boron enriched aluminum surface, novel carbon foam structures and my research interest in 2D electronics. Then I will review the computational techniques I used in the study, mostly DFT methods. In Chapter 3, I introduce an effective way to enhance surface hardness of aluminum by boron nanoparticle implantation. Using boron dimers to represent the nanoparticles, the process of boron implantation is modeled in a molecular dynamics simulation of bombarding the aluminum surface by energetic B2 molecules. Possible metastable structures of boron-coated aluminum surface are identified. Within these structures, I find that boron atoms prefer to stay in the subsurface region of aluminum. By modeling the Rockwell indentation process, boron enriched aluminum surface is found to be harder than the pristine aluminum surface by at least 15%. In Chapter 4, I discuss novel carbon structures, including 3D carbon foam and related 2D slab structures. Carbon foam contains both sp2 and sp3 hybridized carbon atoms. It forms a 3D honeycomb lattice with a comparable stability to fullerenes, suggesting possible existence of such carbon foam structures. Although the bulk 3D foam structure is semiconducting, an sp2 terminated carbon surface could maintain a conducting channel even when passivated by hydrogen. To promote the experimental realization of this novel foam structure, I also propose a growth model. I postulate that preferred growth should occur near the grain boundary of a carbon saturated polycrystal of transition metal. These findings are supported by a calculation of carbon diffusion in the solid. 2D semiconductors of group V elements are discussed in Chapters 5, 6, 7, and 8, including different phosphorus and arsenic structural phases. Structural and electronic properties of bulk and few-layer black phosphorus, so-called phosphorene, are studied in Chapter 5. In Chapter 6, I propose a new 2D structural phase of phosphorus, with the name blue phosphorus related to its wide predicted fundamental band gap. Then I move down in the periodic table and investigate the properties of grey arsenic in Chapter 7. Finally, I propose a tiling model to identify and categorize these structural phases in Chapter 8. 7R P\ ZLIH < Ψ|Ψ > . (2.4) In systems with many electrons, especially in solids with the number of electrons of the order of 1023 , it is impossible to solve the many-body electronic Schr¨odinger equation exactly in nowadays computers. Many applicable approaches with different approximations have been proposed to simplify the complex many-body problem, including the Hartree and Hartree-Fock (HF) self-consistent methods. The Hartree method describes the motion of an electron in the mean field of all the other electrons. While useful to get basic insights into the system, it is not precise enough. Its lack of taking into account the antisymmetric form of the all-electron wave function, a necessary prerequisite to describe the quantum nature of the electrons as Fermi particles, the electronic wave function |Ψ > is represented by a large Slater determinant in the Hartree-Fock method. The shortcoming of this approach is its failure to adequately describe electron correlation. In Quantum Chemistry, this is achieved by considering a set of Slater determinants in the Configuration Interaction (CI) scheme. This approach has been limited to systems of very few atoms and is not useful for infinite 8 solids. A very different approach that does not make use of the all-electron wave function is used in this Thesis. This approach, which considers only the total charge density, is called the Density Functional Theory (DFT) [25]. As I show in the following, DFT is superior to HF by taking into account not only exchange, but also correlation of electrons and their effect on the total energy. This is achieved exactly in model systems and approximately in solids of interest in this Thesis. DFT is suitable to describe solids with translation symmetry and has become the technique of choice to study the fundamental properties of regular solids and nanostructures. Electron density ρ(r), rather than the electronic wave function Ψ, has been introduced as a fundamental quantity in the Hohenberg-Kohn theorems [25], which were first introduced by Hohenberg and Kohn in 1964. After these conceptional foundations have been provided, DFT has been made applicable to solids of interest through the Kohn-Sham approach [26], which maps a many-body problem onto an independent quasi-electron problem by expressing the complex many-body interaction in terms of an exchange-correlation functional of the electron density. Since then, DFT has been very successful to study electronic and structural properties of atoms, molecules, and solids. I will start with briefly reviewing the HohenbergKohn theorems and then introduce the Kohn-Sham approach by deriving the Kohn-Sham equation. Finally, I will briefly introduce the computer codes including SIESTA [27] and VASP [28, 29, 30, 31] that I used in this Thesis to study the structural, mechanical and electronic properties of selected complex nanostructures. 9 2.1 Hohenberg-Kohn theorem DFT is a remarkable theory that, in the electronic ground state, replaces the complicated formalism using a many-body wave function and the related Schr¨odinger equation by a much simpler approach based on the electron density ρ(r). The associated total energy functional then depends only on ρ(r). The idea to use the electron density ρ(r) as the basic variable originates from the Thomas-Fermi model and has been legitimized by Hohenberg and Kohn in the Hohenberg-Kohn theorem. The Hohenberg-Kohn theorem, the heart of DFT, states that in ground state, the external potential v(r), which affects the motion of all electrons, could be determined by the electron density ρ(r), up to a trivial additive constant. The ground-state electron density ρ(r) is sufficient to determine the total energy of the system in the ground state. Thus, in the ground state, it turns out that the total electronic energy in the system can be subdivided into the kinetic energy of electrons T [ρ], the interaction between electrons Vint [ρ], and the interaction with the external potential v(r) as E[ρ] = T [ρ] + Vint [ρ] + v(r)ρ(r)dr. (2.5) Once the expression for the kinetic energy T [ρ] and interaction energy Vint are known, the ground state electron density ρ and total energy E[ρ] could be calculated based on variational principle that minimizes the total energy while keeping the total number of electrons N constant. Vint contains the repulsive Hartree energy and the remaining energy terms related to exchange and correlation interactions, which are the most challenging to determine. A viable method to approximate Vint will be provided in the following Section. 10 2.2 Kohn-Sham equations The Hohenberg-Kohn theorem justifies using the ground state electron density ρ(r) as the variable to determine the ground state properties of the system. However, as exact expressions for the kinetic energy T [ρ] and Vint [ρ] are unknown, it is unclear how to obtain E[ρ]. In the Thomas-Fermi model, a direct approach has been developed to construct approximate forms for T [ρ] and Vint [ρ] as explicit functions of the electron density ρ. These approximations greatly simplify the expressions to evaluate the total energy E[ρ], but are too crude to provide reliable results. In particular, molecules would not bind if described by the Thomas-Fermi model. Instead of using a many-body wave function or considering an over-simplified non-interacting system, Kohn and Sham introduced an ingenious method to project the real Physics of a many-body interacting system onto an independent quasi-electrons system, with the two systems maintaining the same electron density ρ in the ground state. This is the basis of the so-called Kohn-Sham (KS) approach. In this approach, the real many-body wave functions are represented in terms of ρ(r) by the wave functions of non-interacting quasi-electrons. The major part of the kinetic energy could be estimated using the quasi-electrons wave functions |Ψi >, N Ts [ρ] = 1 < ψi | − ∇2 |ψi > . 2 i (2.6) Also the electron density is related to the quasi-electrons wave functions, N ρ(r) = i s 11 |ψi (r, s)|2 , (2.7) where r refers to the spatial and s to the spin coordinates. The Coulomb energy of the quasi-electrons, the Hartree energy, is determined as VHartree = 1 2 drdr ρ(r)ρ(r ) . |r − r | (2.8) Then, ignoring the trivial nucleus-nucleus interaction, the ground state energy functional of the system could be rewritten as EKS [ρ] = Ts [ρ] + VHartree [ρ] + Exc [ρ] + drvext ρ(r). (2.9) In the above equation, the non-trivial many-body interaction is put into the Exc term. Exc [ρ] describes the energy caused by the exchange interaction and by electron correlation in the electron-electron interaction energy Vint [ρ] besides the Hartree energy, and also the difference in the kinetic energy between the interacting electrons and non-interacting quasi-electrons. Thus, we get Exc [ρ] = T [ρ] + Vint [ρ] − Ts [ρ] − VHartree [ρ]. (2.10) The Kohn-Sham independent quasi-particle problem in the ground state could be solved using the variational principle that minimizes the energy functional EKS [ρ] with respect to the electron density ρ(r) while keeping the total charge constant. The quasi-electron wave functions < ψi |ψj > are orthonormal. Then, variational minimization could be performed using Lagrange multipliers, giving δ[EKS [ρ] − i (< ψi |ψi > −1)] = 0. 12 (2.11) This leads to the set of equations 1 [− ∇2 + vef f (r)]ψi (r) = i ψi (r). 2 (2.12) Here, vef f is effective potential depending on ρ and has the form vef f (r) = vext (r) + δVHartree δExc + = vext (r) + vHartree (r) + vxc (r). δρ δρ (2.13) Equations (2.7) and (2.12) are called Kohn-Sham equations. The eigenvalues of the KohnSham equations are just Lagrange multipliers related to the quasi-electrons, rather than the true eigenvalues of the interacting electrons in the system [26]. The merits of the Kohn-Sham equations are that they separate the complex many-body interaction energy into the independent quasi-particle kinetic energy Ts , the Hartree energy VHartree , the interaction with an external potential Vext , and exchange-correlation energy Exc . The first three terms contribute a large portion to the total energy and could be determined directly from the quasi-electrons wave functions. The exchange-correlation functional Exc [ρ] contains many-body interactions, including the electron exchange energy and the electron electron correlation energy. In most implementation, Exc [ρ] could be approximated by a function of ρ(r) and its gradients. Once the formula of Exc [ρ] is determined, the Kohn-Sham equations could be solved in a self-consistent way. 2.3 Exchange-correlation functionals The Kohn-Sham equations are solvable when suitable forms of exchange-correlation functional Exc [ρ] and the corresponding exchange-correlation potential vxc [ρ] are determined. 13 Approximations at different levels have been made to construct suitable forms of exchangecorrelation functionals, including the Local Density Approximation (LDA) [26, 32] and the Generalized Gradient Approximation (GGA) [33]. In LDA, the exchange-correlation energy is a local function of the local electron density ρ. Then, the exchange-correlation energy could be estimated using LDA [ρ] = Exc drρ(r) xc (ρ). (2.14) Here, xc (ρ) is the exchange and correlation energy density of a homogenous electron gas with the density ρ. Consequently, LDA is exact −by construction− for the uniform electron gas. It is approximately correct when the charge density does not change much. Within the LDA, it is further possible to separate the contributions from exchange and correlation as xc (ρ) = x (ρ) + c (ρ). (2.15) In a uniform electron gas, the exchange part obeys the expression xc (ρ) 1 3 3 1 = − ( ) 3 ρ(r) 3 . 4 π (2.16) The remaining part of xc is the correlation energy c , which has been determined accurately using quantum Monte Carlo calculations of a uniform electron gas with different density by Ceperley and Alder [32]. Application of LDA to solids and nanostructures is based on the assumption that the exchange-correlation energy of the nonuniform system could be obtained by separating the system into infinitely many small portions with constant electron density ρ, and by treating these as a uniform electron gas of the same density. In this Thesis, most 14 of the studies related to boron and carbon nanostructures are conducted based on LDA. Although LDA provides a good approximation to the exchange-correlation energy, ignoring the variation of the electron density will lead to errors especially to systems with significant variations in the electron density. To bridge the gap between the uniform electron gas and real systems with varying charge density, it is quite natural to consider Exc also as a functional of the electron density gradient. This is the essential idea of so-called generalized gradient approximation (GGA). Under the approximation of GGA, the exchange-correlation energy has the form GGA [ρ] = Exc drf (ρ(r), ∇ρ(r)). (2.17) It is useful to mention that the exchange-correlation energy is not just a conventional gradient expansion in a power-series, but rather a general function of the electron density and its gradient. There are many different GGAs depending on the choice of the function f (ρ, ∇ρ). Among these, the currently most popular is the PBE [33] functional, proposed by Perdew, Burke and Ernzerhof in 1996. In this Thesis, GGA-PBE is used in studies of 2D phosphorene and related systems. When compared to LDA, GGA-PBE predicts structural and electronic properties that agree well with experimental data. LDA and GGA could give satisfactory results for structural and energetic properties of solids and nanostructures. However, none of them could describe adequately the long-range dispersion (or van der Waals (vdW)) interaction in a general system very well. In layered systems, such as graphite and layered phosphorus, which are studied in this Thesis, there is a significant contribution of vdW interactions to the inter-layer bonding. In these systems, the inter-layer bonding is typically overestimated in LDA and underestimated in GGA-PBE. Various corrections and methods to introduce vdW interaction in DFT have been developed 15 recently, such as the DFT-D2 [34] and optB86b-vdW [35, 36] functionals. To represent the correct Physics of the inter-layer interactions in such so-called vdW-materials, I have also made use of such functionals. In the Kohn-Sham equations (2.12), i has been introduced as a Lagrange multiplier with no particular significance. i with the dimension of energy is often confused with the energy eigenvalue in a one-electron Schr¨odinger equation. It turns out that while this is strictly approximate, the energy spectrum of i often resembles the electronic density of states. When this is done, both LDA and GGA underestimate the fundamental electronic band gap of semiconductors. Methods like GW and hybrid functionals have been developed to better reproduce the band gap in the solids. The GW method determines the eigenvalues of the self-energy operator, which is expanded as the product of the single particle Green’s function G and the dynamically screened Coulomb interaction W . While predicting the electronic band gap value in good agreement with experiments, the calculation is rather time consuming. Hybrid functionals methods are an alternative, that combines partial HartreeFock exact exchange with partial exchange-correlation energy of DFT. The mixture of HF, which typically overestimates the band gap, with DFT, which typically underestimates the band gap, usually predicts a reasonable electronic band gap value for the system. In this Thesis, to predict the electronic band gap value of certain nanostructures, I have made use of hybrid functional methods, in particular the HSE06 approximation introduced by Heyd, Scuseria, and Ernzerhof [37]. This approach uses a Coulomb potential screened by an error function in the short range. It then expresses the exchange correlation energy as Exc = aExHF,SR (w) + (1 − a)ExP BE,SR (w) + ExP BE,LR (w) + EcP BE . 16 (2.18) P BE,SR The exchange energy of PBE is separated into a short-range term Ex P BE,LR range term, Ex and a long- . The short-range term of PBE is mixed with HF exact exchange energy. a is the mixing parameter between PBE and HF and w controls the short range of the interaction. In HSE06, a = 0.25 and w = 0.2, but these quantities are often treated as variable parameters. 2.4 Implementation of the DFT method As long as the expression of Exc is selected, the Kohn-Sham equations can be solved in a self-consistent manner. The initial electron density ρ0 (r) is guessed as an input, typically obtained from a superposition of atomic charge densities. A given density ρ(r) is used to determine the effective potential vef f . Then the Kohn-Sham equations are solved to get a set of eigenstates, which will generate a new electron density ρ(r). This new electron density is used to generate a new potential vef f . Then, the Kohn-Sham equations are solved iteratively until the self-consistency is reached. The computational procedure described above is the general concept to solve the KohnSham equations. To deal with real systems, including periodic solids, further approximations and technical treatments are needed. In a periodic solid or a nanostructure, a set of basis functions is selected to expand the eigenfunctions of the system and construct the Kohn-Sham matrix. Then the process of solving the Kohn-Sham equations is to diagonalize the KohnSham matrix. The orbital basis could be plane wave functions or localized atomic orbitals. Moreover, as the core electrons does not affect the properties of solid much, their effects on the valence electrons could be represented by a pseudopotential, which reduces the calculation efforts substantially. There are many different formulas to construct pseudopotentials, 17 including ultrasoft pseudopotential [38] and norm-conserving pseudopotential [39, 40]. In this Thesis, the theoretical calculations are conducted using the SIESTA [27] and VASP [28, 29, 30, 31] codes. SIESTA uses a localized numerical basis and norm-conserving pseudopotential, so it is relatively fast and capable of producing pretty precise results for the targeted systems. Most of the calculations in this Thesis, which address the structural and electronic properties of nanostructures using LDA or PBE, are conducted by SIESTA. VASP uses a projector augmented wave method (PAW) with an PAW basis and pseudopotentials [31, 41]. Besides providing LDA and PBE functionals like SIESTA, VASP provides advanced implementations of hybrid functionals (HSE06) and van der Waals corrections. Calculations in this Thesis related to the latter issues are mainly done by VASP. 18 Chapter 3 Enhancing mechanical hardness of aluminum surfaces by nanoboron implantation The following discussion is closely related to a publication by Zhen Zhu, Dae-Gyeon Kwon, Young-Kyun Kwon and David Tom´anek, Enhancing mechanical toughness of aluminum surfaces by nanoboron implantation: An ab initio study, Chem. Phys. Lett. 620, 25-28 (2015) [1]. This study is a collaboration with Prof. Young-Kyun Kwon’s group at Kyung Hee University. 3.1 Introduction Aluminum owes its widespread use especially in the aerospace industry to its low gravimetric density, good machinability, and high strength-to-weight ratio[8, 9]. Since this metal is rather soft, its low wear resistance poses a serious problem especially in harsh environments. The common way to enhance surface hardness of aluminum is by hard-coat anodizing, an electrolytic passivation process, which forms a brittle surface oxide layer that may chip off. Alternative ways to enhance hardness while maintaining ductility mostly involve bulk modifications, such as changing from Al to the Al/TiB2 composite [10]. Hard diamond-like 19 carbon coatings do not adhere well to the soft Al surface without adhesion promoters at the interface [11]. An intriguing alternative to carbon involves boron that is nearly as hard as diamond and that may bond more strongly to aluminum, since Al belongs to the boron group in the periodic system. Boron is known for its large flexibility in bonding, forming unusual structures in the bulk [12] and in nanoparticles [13]. The ability of boron to harden surfaces of metals including Al has been suggested by the observed increase in their abrasion resistance following bombardment by boron nanoparticles [14]. Here I study the possibility of enhancing the surface hardness of Al by implantation of boron nanoparticles using ab initio density functional calculations. My results include the equilibrium structure, stability, elastic properties and formation dynamics of a boronenriched Al surface after exposure to energetic boron nanoparticles. Combining molecular dynamics simulations with structure optimization studies, I identify structural arrangements that optimize the formation of strong covalent B-Al bonds. Nano-indentation simulations based on constrained optimization suggest that presence of boron aggregates enhances significantly the mechanical hardness and wear-resistance of aluminum surfaces. 3.2 Computational methods My calculations are based on the ab initio density functional theory (DFT) as implemented in the SIESTA code [27]. I used the Ceperley-Alder [32] exchange-correlation functional as parameterized by Perdew and Zunger [42] and norm-conserving Troullier-Martins pseudopotentials [39]. My basis consisted of pseudo-atomic orbitals (PAOs) generated by the split-valence scheme for a double-ζ polarized basis set. All calculations were performed using periodic boundary conditions. The (111) surface of aluminum was represented by a 20 five-layer slab, shown in Fig. 3.1(a), with each layer containing four Al atoms per unit cell and the slabs separated by a vacuum region of 10 ˚ A. I sampled the small quasi-2D Brillouin zone associated with the large supercells by a dense 8×8×1 k-point mesh [43]. The energy shift due to the spatial confinement of the PAO basis functions [44, 45] was limited to less than 10 meV. The charge density has been determined self-consistently on a real space mesh with a very high cutoff energy of 180 Ry, sufficient for total energy convergence to within 1 meV/atom. I studied bombardment of Al(111) by boron nanoparticles using microcanonical molecular dynamics (MD) calculations. I integrated the equations of motion with the Verlet algorithm, using Δt = 0.5 fs as time step to cover very long simulation periods of up to 2 ps. The geometry of selected structures was optimized using the conjugate gradient method. A structure was considered optimized when none of the residual forces exceeded 0.01 eV/˚ A. 3.3 Results and discussion The key objective of my study is to identify unusually stable, covalently bonded B-Al nanostructures in the surface region of Al following implantation of boron nanoparticles, which would improve wear resistance of Al. Correct description of bonding changes requires treatment by a self-consistent electronic structure calculation. This approach is computationally very demanding and, when combined with molecular dynamics simulations, poses a serious limit on the number of atoms per unit cell. Also, I will represent boron nanoparticles by B2 molecules when addressing implantation. My focus will be to investigate favorable bonding changes in specific structural arrangements, which has been neglected in previous studies. While obtaining this microscopic information, I will not be able to provide an adequate de- 21 scription of compressive strain changes at the surface, the classic mechanism for improving hardness, toughness, and wear at a metal surface. Before studying implantation of boron nanoparticles, I optimized the structure of the Al(111) surface. I found that slabs with at least five Al layers are required to reproducibly represent the surface relaxation of the topmost layers. To adjust the interface to the optimum bulk structure, the bottom two layers have been constrained in the bulk geometry. In the optimum geometry, shown in Fig. 3.1(a), the largest change in the inter-layer separation di,i+1 occurs for the topmost layer i = 1 and decreases with increasing depth below the surface as Δd12 /d12 = +2.1%, Δd23 /d23 = +1.2%, and Δd34 /d34 = +1.2%. Whereas the topmost layer contracts at most metal surfaces [46], I find a surface expansion in slabs with more than 3 layers in agreement with experimental observation [47], indicating that the close-packed Al(111) surface does not benefit from Smoluchowski smoothing [48]. Next, I need to locate stable atomic arrangements in a boron enriched Al(111) surface that would locally enhance surface hardness. This task is very hard, since there is no empirical guidance for lack of a well-behaved potential energy surface representing boron interacting with aluminum. Global optimization of the positions of 20 Al and additional B atoms in each unit cell, when viewed as a search in configurational space of 60 or more dimensions, appears nearly impossible using DFT due to the high computational requirements. To represent the experimental conditions, which partly mimic simulated annealing, I performed molecular dynamics simulations of energetic boron nanoparticles bombarding the Al(111) surface at initially T = 0 K. For the sake of simplicity, I limited my MD study to a boron dimer with kinetic energy 0.5 eV/atom that impacted on the surface vertically. The results indicate that the B2 molecule transfers its excess kinetic energy to locally melt the aluminum surface. Even 22 (a) (b) B Side View (c) d12 d23 d34 [111] Al Top View [112] [110] - [112] Ĭ=0 Ĭ=0.5 Ĭ=2 Figure 3.1 For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this Thesis. Optimized structure of pristine and B enriched Al(111) surfaces in side view (top panels) and top view (bottom panels). (a) Pristine Al(111) surface prior to bombardment by boron dimers, shown above the surface. The most stable structure of the Al surface enriched with (b) 2 B atoms and (c) 8 B atoms per unit cell containing 20 Al atoms. The unit cells used in the study are highlighted by color and shading. Reproduced from [1], c 2015 Elsevier. though the B2 atoms may separate in the hot surrounding Al metal, they eventually remain connected to each other in subsurface sites, as seen in Fig. 3.1(b). Whereas the total energy is conserved during the molecular dynamics simulation, the continuous transformation between more stable and less stable structures is reflected well in the time-dependence of the potential energy V , which is presented in Fig. 3.2(a). The time evolution of the system may be viewed as a trajectory in configurational space, and may be compared to the trajectory during a simulated annealing optimization, with potential energy minima representing particularly stable geometries. Since potential energy minima correspond to kinetic energy maxima, high velocity may prevent atoms from probing closely the most stable arrangements. Assuming that minima in the potential energy during my MD simulation are close to optimum atomic arrangements, I froze the corresponding structures 23 (a) 0 Ĭ=0.5 ǻEtot (eV) ǻV (eV) -2 (b) 0.0 -4 -6 -8 B CD -10 A 0 0.5 E F 1.0 H G 1.5 I 2.0 Time (ps) -0.2 -0.4 -0.6 -0.8 -1.0 -1.2 -1.4 -1.6 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 B exposure Ĭ Figure 3.2 Stability of boron enriched Al(111) surfaces. (a) Potential energy evolution during the microcanonical MD simulation of boron dimers impacting on the Al(111) surface. Solid circles below selected potential minima represent optimized structures at T = 0 K. (b) Diagram depicting the relative stability change ΔEtot of the B enriched Al(111) surface as a function of the B exposure Θ. The green dashed line is a guide to the eye. Changes in the potential energy ΔV and stability are given per unit cell. Reproduced from [1], c 2015 Elsevier. and optimized them using the conjugate gradient technique. I believe that these structures, labeled A − I and shown by the solid circles in Fig. 3.2(a), represent some of the most stable geometries of B enriched Al(111). I found that for the the most stable among these optimized structures, boron atoms generally occupy interstitial sites below the topmost layer and are connected to each other. This optimum bonding geometry reflects the bond strength hierarchy, with the B-B bonds being the strongest and the B-Al bonds being in-between the B-B and the weakest Al-Al bonds. Thus, the structural priority is to maintain a contiguous optimum B structure that also maximizes bonding to Al, possibly at the expense of disrupting Al-Al bonds. I expect that this can be achieved best by placing boron aggregates below the topmost Al layer. Indeed, I find that structure C in Fig. 3.2(a), which represents the most stable atomic arrangement in the B2 Al20 unit cell and is shown in Fig. 3.1(b), agrees with this stability criterion. 24 With the above working hypothesis for optimum atomic arrangements at hand, I next discuss ways to enhance the stability of the Al(111) surface by changing the concentration of boron in the surface region. The corresponding results are presented in Fig. 3.2(b), where I define the boron exposure Θ by the ratio of the total number of B atoms N and the number of surface Al atoms Ns (Al)=4 in the unit cell of the slab. I next define the stability change associated with the presence of B atoms by the total energy difference per unit cell ΔEtot = [Etot (BN Al20 ) − Etot (Al20 ) − N Etot (Bref )]/N . I chose to use Etot (Bref ) = Etot (B12 )/12 as the boron reference energy, since B12 icosahedra are not only very stable as nanoparticles, but also occur as building blocks of the bulk structure of elemental boron [12], α−B and β−B. Then, the sign of ΔEtot should reflect the energetic preference of boron atoms to either form stable isolated nanoparticles or rather to form nanostructures embedded in the metal matrix. My results indicate that structures with Θ = 0.5 are generally less stable than systems with a higher B concentration. Among the systems investigated in my study, I found the Θ = 2.0 structure depicted in Fig. 3.1(c) to be the most stable. The B8 Al20 surface compound benefits from the stability of the contiguous planar boron honeycomb lattice sandwiched between the topmost and the second Al layer, which locally resembles the structure of the stable AlB2 compound. As seen in Fig. 3.1(c), my calculations indicate a stability optimum at Θ≈2.0, corresponding to B8 Al20 . Considering stability changes per boron atom, as reflected in the definition of ΔEtot above, the B8 Al20 structure is more stable by 0.6 eV/B atom than the boron poorer B2 Al20 and by 0.3 eV/B atom than the boron richer B12 Al20 system. As suggested in the beginning, the reason for very strong bonding between B and Al is that both elements are in the same group of the periodic system. In this case, I expect only a moderate charge transfer and covalent bonds between neighboring B and Al atoms. Results 25 DOS (state D es/eV/unit cell) DOS (stattes/eV/unitt cell) D 8 (a) 6 Al 4 2 0 -15 -10 -5 EF 5 10 15 8 (b) Total Al B 6 4 2 0 -15 -10 -5 EF 5 10 15 Energy (eV) Energy (eV) (c) ǻn (el/Å3) +0.0172 +0.0114 +0 0114 +0.0057 0.0000 -0.0057 -0.0114 0 0114 -0.0172 Figure 3.3 Effect of B-Al bonding on the electronic structure. The electronic density of states (DOS) of (a) a pristine Al(111) slab and (b) the slab structure of Fig. 3.1(b) with Θ = 0.5. The partial DOS of Al is shown by the dashed and that of B is shown by the dotted line in (b). (c) Charge density difference in the system of (b), presented as a contour plot in a plane normal to the surface, indicated by the dashed line in Fig. 3.1(b). Reproduced from [1], c 2015 Elsevier. related to the electronic structure of B enriched Al are presented in Fig. 3.3. The electronic density of states (DOS) of the topmost 2 layers of the Al(111) surface, shown in Fig. 3.3(a), is close to that of a free-electron metal. The corresponding DOS of the topmost 2 layers of the B-enriched Al(111) surface, with the structure depicted in Fig. 3.1(b), is presented in Fig. 3.3(b). Comparing these two densities of states, including the projection onto Al and B sites in Fig. 3.3(b), indicates only a moderate perturbation of the Al subsystem by bonding to boron. Mulliken population results with a single-ζ basis indicate net transfer of ≈0.7 electrons from Al to B in the structure of Fig. 3.1(b). To get a more accurate idea about the charge redistribution in the system caused by the presence of B in Al(111), I 26 plotted in Fig. 3.3(c) the electron density difference Δn(r) = nBx Al1−x (r) − nAl (r) − nB (r). Taking nAl (r) and nB (r) as the electron densities of the isolated, frozen subsystems, Δn(r) is nonzero only in regions, where the total electron density deviates from a mere superposition of the subsystem charge densities. My calculated Δn(r) in Fig. 3.3(c) indicates a moderate electron transfer from Al to B, in agreement with the Mulliken analysis, with many Al neighbors of B contributing as electron donors. Close inspection of Fig. 3.3(c) reveals that presence of Al does not change the B-B bond and vice versa. A large electron accumulation in the region of Al-B bonds suggests that Al-B bonds are covalent and should be strong as anticipated. (a) F (c) d (Å) d (b) Soft Ĭ=0 Ĭ = 0.5 Ĭ = 2.0 F d Hard F (eV/Å) Figure 3.4 Effect of nano-boron on the surface hardness of Al. Schematic model of a Rockwell nano-indenter in (a) is compared to its atomic-scale counterpart in (b). (c) Displacement d as a function of the force F pressing vertically on a surface atom in Bx Al(111). Results are presented for structures with Θ = 0, Θ = 0.5 and Θ = 2.0, depicted in Fig. 3.1(a-c). Reproduced from [1], c 2015 Elsevier. There is no easy way to provide a realistic description of surface hardness, as probed by nano-indentation, in atomistic calculations, due to the complexity of processes including plastic flow associated with dislocation motion along slip planes in the sample. Instead of explicitly considering large-scale displacement of atoms during indentation, I note that plastic 27 flow is initiated at defects that require particularly low activation barriers for nucleation. It has been shown that such low nucleation barriers are commonly found at defects that form within the anharmonic regime of soft harmonic modes [49]. This correlation provides a direct link between softness in the elastic response and plastic activity. In this case, I should be able to get valuable insight into surface hardness by probing the elastic response only. In the following, I study this important quantity by replicating the way hardness is probed experimentally using nano-indentation and focussing on the initial elastic response. The schematic setup of the commonly used Rockwell nano-indentation test is shown in Fig. 3.4(a). Rockwell hardness is determined by the indentation depth d of a hard, conical nano-indenter that is pressed by force F towards the surface. The deeper the indentation at a given load, the softer the material. As an atomic-scale analogy, I studied structural rearrangements introduced by a vertical displacement of a surface atom by d with respect to the bottom of the slab, as illustrated schematically in Fig. 3.4(b). For different values of d, I optimized the position of all atoms in the unit cell except the constrained atom, considered a ‘nano-indenter’, and the bottom of the slab. I then interpreted the vertical component of the Hellmann-Feynman force acting on the nano-indenter as the load F associated with the indentation d. My results for the atomic nano-indentation process are presented in Fig. 3.4(c) for the three structures shown in Fig. 3.1, representing different degrees of boron enrichment. These results indicate a linear relationship between the load and the indentation depth for d below 0.5 ˚ A. Defining local surface hardness by the F/d ratio, I find the values 1.6 eV/˚ A2 or 25.6 J/m2 for pure Al(111) (Θ = 0), 1.8 eV/˚ A2 or 28.8 J/m2 for a low nano-boron exposure Θ = 0.5 and 2.4 eV/˚ A2 or 38.4 J/m2 for the higher nano-boron exposure Θ = 2.0. Considering the fact that the net amount of nano-boron is rather small, the increase of the surface 28 hardness by 13% for Θ = 0.5 and by 50% for Θ = 2.0 is formidable. I must note that my five-layer slab is a very limited representation of a realistic polycrystalline surface, where energetic boron nanoparticles may penetrate deeper, modifying the grain structure and reinforcing intergranular interfaces by covalent bonds. I wish my study to inspire corresponding experiments. 3.4 Summary and Conclusions In conclusion, I used ab initio density functional calculations to study the equilibrium structure, stability, elastic properties and formation dynamics of a boron-enriched Al(111) surface. I used molecular dynamics simulations to model the implantation of energetic boron nanoparticles in Al and identified structural arrangements that optimize the formation of strong covalent B-Al bonds for different concentrations of boron in the surface region. Nanoindentation simulations based on constrained optimization suggest that presence of contiguous boron nanostructures in the subsurface region may increase the mechanical hardness of aluminum surfaces by up to 50% at relatively low boron exposures. 29 Chapter 4 Carbon foam: structure, growth and electronic properties The following discussion is closely related to publications by Zhen Zhu, and David Tom´anek, Formation and stability of cellular carbon foam structures: An Ab initio study, Phys. Rev. Lett. 109, 135501 (2012) [3] and also Zhen Zhu, Zacharias Fthenakis, Jie Guan, and David Tom´anek, Topologically protected conduction state at carbon foam surfaces: An em Ab initio study, Phys. Rev. Lett. 112, 026803 (2014) [2]. 4.1 Introduction The last few decades have witnessed an unprecedented interest in carbon nanostructures, the most prominent of them being fullerenes [15], nanotubes [17] and graphene [18]. Previously postulated hybrid carbon nanofoam structures [50, 51, 52, 53] with a mixed sp2 /sp3 bonding character have received much less attention for lack of direct experimental observation. The growing body of information about the formation of carbon nanostructures including graphene [54], nanotubes [55, 56] and fibers [57] on transition metal surfaces with a particular morphology suggests ways that should favor the formation of particular carbon allotropes. I propose that previously unseen nanostructures including carbon foam may form under specific conditions on a metal substrate. 30 Inspired by previously postulated carbon foams [51, 50, 52, 53], I explore ways to grow such structures on a carbon saturated metal substrate. I use ab initio density functional calculations to investigate the equilibrium structure, structural and thermal stability and elastic properties of the growing system. The foam structures I study, which have a mixed sp2 /sp3 bonding character and resemble a bundle of carbon nanotubes fused to a contiguous 3D honeycomb structure, are rather stable even as slabs of finite thickness. The foam structure may be compressed more easily by reducing the symmetry of the honeycombs. It may accommodate the same type of defects as graphene at little energy cost, and its surface may be stabilized by terminating caps. I postulate that the foam could form under nonequilibrium conditions near grain boundaries of a carbon-saturated metal surface and should remain stable until T > 3, 500 K. Moreover, as electronic properties of foams [58, 51, 50, 59] have received much less attention than their structural stability in spite of the obvious possibility to fine-tune the fundamental band gap value in-between zero in sp2 -bonded graphene and 5.5 eV in sp3 bonded diamond by modifying the foam morphology. Here I also show my results of ab initio electronic structure and quantum conductance calculations indicating the emergence of conduction at the surface of semiconducting carbon foams. Occurrence of new conduction states in these systems is intimately linked to the topology of the surface and not limited to foams of elemental carbon. My interpretation based on rehybridization theory indicates that conduction in the foam derives from first- and second-neighbor interactions between p orbitals lying in the surface plane, which are related to p⊥ orbitals of graphene. The topologically protected conducting state occurs on bare and hydrogen-terminated foam surfaces and is thus unrelated to dangling bonds. My results for carbon foam indicate that the conductance behavior may be further significantly modified by surface patterning. 31 (a) bulk foam sp3 sp2 sp3sp2 sp2 sp3 sp3 sp2 sp3 terminated slab (b) (c) sp3 sp2 2 2 sp3 sp sp terminated slab sp3 sp 2 Figure 4.1 Geometry of carbon foam and thin foam slabs. (a) Structure of bulk foam and that of an individual foam cell in side and top view, allowing to distinguish sp2 and sp3 sites. The dashed line shows the long cell axis. The perspective view in the right panel depicts a larger bulk segment. Structure of (b) an sp3 -terminated and (c) an sp2 -terminated thin foam slab. The tilted view used in left panels depicts the structure and primitive unit cells. The right panels in (b) and (c) are side views of the structure that better illustrate the type of termination and illustrate partial hydrogen coverage. Reproduced from [2], c 2014 American Physical Society. 4.2 Computational Methods My numerical results for the equilibrium structure, stability and electronic properties of carbon foam slabs are based on density functional theory (DFT) as implemented in the SIESTA code [27] and VASP [28, 29, 30]. The different foam surfaces are represented by a periodic array of slabs, separated by a 15 ˚ A thick vacuum region. I used the CeperleyAlder [32] exchange-correlation functional as parameterized by Perdew and Zunger [42], 32 norm-conserving Troullier-Martins pseudopotentials [39], and a double-ζ basis including polarization orbitals. The reciprocal space was sampled by a fine grid [43] of 16×16×1 k-points in the Brillouin zone of the primitive surface unit cell. I used a mesh cutoff energy of 180 Ry to determine the self-consistent charge density, which provided us with a precision in total energy of ≤2 meV/atom. Equilibrium structures and energies based on SIESTA have been checked against values based on VASP code and the hybrid functional calculations [37] are done in VASP. Transport properties of the slabs were investigated using the nonequilibrium Green’s function (NEGF) approach as implemented in the TRAN-SIESTA code [60]. For structures optimized by DFT, I used a single-ζ basis with polarization orbitals, the same 180 Ry mesh cutoff energy, and a 8×8×1 k-point grid [43]. 4.3 Structure and stability The bulk carbon foam, depicted in Fig. 4.1(a), is a cellular structure resembling vaguely a fused triangular array of (6,0) zigzag nanotubes. In contrast to a nanotube array, the walls of the foam cells consist of 60% sp2 bonded atoms shared by two neighboring cells and 40% sp3 bonded atoms shared by three adjacent cells. Density-functional based tight-binding (DFTB) results indicate that the bulk structure is a semiconductor [50] with a band gap of 2.55 eV. Cleavage normal to the long cell axis may generate two different surfaces. I distinguish the sp3 surface terminated by C atoms, which were fourfold coordinated in the bulk, from the sp2 surface terminated by atoms that were threefold coordinated. For computational reasons, I will represent the surface by slabs of finite thickness with two identical surfaces 33 (a) (b) sp3 ˂Q HO‡ +0.087 sp2 +0.052 +0.017 0.000 -0.017 -0.052 -0.087 sp3 sp2 (c) (d) ˮ൹ ˮൻ HO‡ +0.12 +0.06 0.00 -0.06 -0.12 2.52 Å Figure 4.2 Structure and electronic properties of cellular carbon nanofoam. (a) Left panels depict individual cells of the foam in top and side view. Right panel shows the contiguous foam in top view, with individual cells terminated by different types of caps. (b) Electron density difference Δn(r) in a plane normal to the surface, indicated by the dotted line in (a). (c) Side view of the structure of a stable minimum-thickness foam slab. (d) Spin density distribution ρ↑ − ρ↓ in the structure shown in (c), represented in the same plane as in (b). A3 . Reproduced from [3], c 2012 American The isosurfaces are shown for ρ↑ −ρ↓ = ±0.05 el./˚ Physical Society. that are either bare or terminated by hydrogen. The optimum structure of the thinnest foam slab with sp3 termination is shown in Fig. 4.1(b) and that of the thinnest sp2 -terminated slab in Fig. 4.1(c). Fig. 4.2(a) shows a carbon foam structure with possible capping and its building block. In top view, it closely resembles the graphene honeycomb lattice with two important distinctions. I find the optimum lattice constant in the honeycomb plane of the foam to be a = 4.81 ˚ A, which is about twice the graphene value a = 2.46 ˚ A. More important, 1D carbon-carbon bonds in the 2D graphene structure correspond to 2D walls in the infinite 3D foam structure. The foam cells, shown in the left panels of Fig. 4.2(a), closely resemble (6,0) 34 carbon nanotubes The gravimetric density of the optimized foam structure, ρ = 2.4 g/cm3 , lies in-between the experimental values [61] for graphite, ρ = 2.27 g/cm3 , and diamond, ρ = 3.54 g/cm3 . I find the 3D carbon foam structure to be less stable than graphene by ΔEcoh ≈0.42 eV/atom, which is comparable to the C60 fullerene. In Fig. 4.2(b) I display the electron density difference, defined by Δn(r) = ntot (r) − natom (r) as the difference between the total electron density ntot (r) and the superposition of atomic charge densities natom (r). Charge accumulation in the bond region indicates strong covalent bonding especially between neighboring sp2 atoms. My DFT results for the electronic structure indicate that the bottom of the conduction band lies below the top of the valence band, suggesting that the infinite foam should be metallic. In reality, this finding is a well-known artifact of DFT that I correct using the hybrid functional method (HSE06) [37], which indicates semiconducting behavior of the bulk structure. The details about electronic structure will be discussed in the following sections. Besides the bulk structure, I also considered and optimized foam slabs of different thickness. I must take into account the fact that the surface terminated with sp2 -type atoms, which are shared by two honeycombs, is inequivalent to a surface with sp3 -type atoms, which are shared by three honeycombs. The thinnest stable free-standing slab, dubbed the ‘single-decker’ structure and shown in Fig. 4.2(c), has both surfaces of the sp3 -type. It has some commonalities with graphitic nanostructures that show magnetic ordering at zigzag edges [62, 63, 64, 65, 66]. Similar to the narrowest zigzag graphene nanoribbon, my system displays a flat band near EF that gives rise to spin polarization with antiferromagnetic coupling across the slab, as seen in Fig. 4.2(d). The dominating role of the surface reduces the stability of the ‘single decker structure’ by ΔEcoh = 0.95 eV/atom with respect to the bulk carbon foam. I note an even-odd alternation in the energy as a function of slab thickness in 35 terms of the number of hexagon rows, since the slab surfaces may be either identical or different. In any case, the role of the surface decreases with increasing slab thickness, and reaches a much smaller value ΔEcoh ≈0.46 eV/atom in the ‘triple-decker’, shown in Fig. 4.2(b), than in the ‘single-decker’ structure. The energy penalty due to unsaturated surfaces may be significantly reduced if the slab is attached to a substrate, or if the cells are covalently terminated by caps, similar to the dome termination of carbon nanotubes. I considered either a hexagon or two adjacent pentagons as candidate caps to terminate the honeycombs, as seen in the right panel of Fig. 4.2(a). Both caps have 6 twofold coordinated C atoms at the edge that may form covalent bonds with the surface atoms. Assuming that all honeycombs on one side are capped and using A = 20.04 ˚ A2 for the area of each honeycomb, I estimated the surface energy reduction associated with cap termination to be ΔEs = −1.03 eV/˚ A2 for hexagonal caps and ΔEs = −0.25 eV/˚ A2 for the less-stable two-pentagon caps. I need to note that this stabilization energy contains the termination energy of both the surface and the individual unsaturated caps, and that these energy terms can not easily be separated. Since epitaxy is an issue when considering the possibility of foam growth on a metal substrate, I investigated the lateral compressibility of the foam structure. My definition is analogous to the elastic response of a uniform isotropic 3D structure with volume V to hydrostatic pressure P = F/A, given by the force F per area A, which is represented by the bulk modulus B = −V (∂P/∂V )T . The elastic deformation of the area A within a 2D slab structure subject to in-plane hydrostatic pressure P2D = F/l, given by the force per length l, can be represented by an analogous 2D bulk modulus, defined by B2D = −A(∂P2D /∂A)T . Of course, I expect B2D to be nearly proportional to the slab thickness. I find this value to be quite useful, since it allows to determine the critical slab thickness for epitaxial growth 36 (a) (c) (b) 6 6 6 6 5 7 5 7 6 5 6 5 8 5 6 5 6 Figure 4.3 Defects in the foam. (a) Folding of the perfect foam, induced by applying hydrostatic pressure or by electron doping. Foam structures containing (b) 5775 and (c) 558 defects, familiar from defective graphene. Reproduced from [3], c 2012 American Physical Society. on a particular incommensurate substrate. Applying hydrostatic pressure in the plane of the layer, I find that the honeycomb structure may be compressed more easily by breaking the honeycomb symmetry than by uniformly compressing the honeycombs. The structure of the deformed foam, depicted in Fig. 4.3(a), indicates the preferential way the foam may fold. For this elastic response, I find B2D = 112.9 N/m in the ‘single-decker’ and B2D = 163.9 N/m in the ‘triple-decker’ structure. For the sake of comparison, when considering a very thick slab of thickness h, I used the bulk calculation to obtain B3D ≈B2D /h = 178 GPa. I find this value to be much smaller than that of the ideal structure with suppressed folding, which had been studied previously [50] with results similar to my value B2D /h = 299.4 GPa. Even though the possibility of folding reduces the bulk modulus, finite compressibility should still play a significant role during foam growth on lattice-mismatched or defective substrates. Interestingly, I find that foam folding occurs spontaneously when the system is doped by electrons. The structure presented in Fig. 4.3(a) can be obtained by either applying isotropic pressure in 2D or, at zero pressure, by doping with 0.2 electrons per C atom. In the latter 37 case, I find that folding induced by doping reduced the foam energy by 0.19 eV/atom for the bulk structure. I also find that the proposed foam structure may accommodate a similar type of defects as graphene with the main difference that bond rotations in graphene correspond to wall rotations in the foam. In graphene monolayers, lines of 5775 or Stone-Thrower-Wales [67, 68] and of 558 defects have been observed to accumulate near grain boundaries [69, 70, 71] and step edges [72]. Their presence reduced stress in strained free-standing layers and the lattice mismatch energy in adsorbed layers, which in this way maintained their epitaxy over large areas. The analogous 5775 or 558 defect structures in the foam are depicted in Fig. 4.3(b) and 4.3(c). Since the foam structure is rather flexible, the energy penalty associated with these types of defects is relatively small, amounting to ΔE = 0.19 eV/atom for the 5775 structure of Fig. 4.3(b) and ΔE = 0.20 eV/atom for the 558 structure of Fig. 4.3(c) with respect to the perfect infinite honeycomb lattice. With a bulk modulus B≈250 GPa, the defective 5775 and 558 foam structures are slightly more compressible than the perfect foam with suppressed cell folding. Similar to supported graphene, these types of defects should reduce the lattice mismatch energy on a particular substrate caused by different lattice constants or, on a polycrystal, across grains with different orientation. To find out whether the carbon foam may or may not decompose to a more stable allotrope under growth conditions, I studied its thermal stability by performing molecular dynamics (MD) simulations in the temperature range 500 K< T <5, 000 K. To avoid artifacts caused by small unit cells, I used supercells containing 160 carbon atoms. For these large unit cells, I used the Tersoff bond-order potential [73] in molecular dynamics simulations covering time periods of 10 ps using 0.5 fs time steps. My results indicate that the infinite foam should be stable up to a high melting temperature near 3, 700 K. Even though free38 standing slabs with finite thickness may be thermally less stable, termination by caps or attachment to a substrate should increase their thermal stability. 4.4 Growth mechanism Inspired by the observed growth of graphene [54] and carbon nanotubes [55] on cobalt saturated with carbon, I studied possible growth pathways of the foam on this substrate. To get insight into the foam-substrate interaction including optimum lattice registry, I represented the Co(0001) surface by a four-layer slab with the two bottom layers constrained in the bulk geometry. Besides the perfect Co(hcp) lattice, I also considered fcc layer stacking when discussing grain boundaries. I considered different foam terminations at the interface in order to find the optimum interface geometry. I found that the sp2 -type terminated foam surface attaches more strongly to Co(0001) than the sp3 -type terminated surface. The largest reduction of the foam surface energy by ΔEs = −0.75 eV/˚ A2 occurs, when surface C atoms occupy the hollow sites. I should note that this stabilization energy reflects the reduction of both the metal and the foam surface energy. Since a realistic representation of the growth mechanism by molecular dynamics (MD) simulations is currently not possible due to time limitations, I discuss in the following likely processes that should contribute to foam growth and judge their importance according to potential energy surfaces. To favor foam growth, I need to find a suitable substrate geometry and identify growth conditions that promote the formation of foam rather than other competing nanostructures [54, 55, 74]. Assuming that the feedstock are carbon atoms dissolved in the substrate, I consider grain boundaries and steps as preferential nucleation sites of the foam. Three competing processes contribute to the nucleation and growth of carbon 39 Energy (eV) (b) (a) hcp Top view _ [1100] b fcc b hcp o (c) fc o fc o bc fc bc fc bc _ [1120] A Side view B A [0001] _ [1120] B Figure 4.4 Surface and bulk diffusion of C atoms on a carbon saturated Co(0001) surface. Surface diffusion in (a) is compared to bulk diffusion in (b) and diffusion along a grain boundary in (c). The top panels represent energy changes per atom along the optimum diffusion path, which is indicated by the dashed line in top and side views, presented in the bottom two panels. Reproduced from [3], c 2012 American Physical Society. nanostructures on the surface: surface diffusion of carbon, bulk diffusion of carbon inside individual grains, and bulk diffusion along grain boundaries that had not been considered previously. My results for these three processes are presented in Fig. 4.4. Since surface diffusion of C atoms, depicted in Fig. 4.4(a), does not require displacement of substrate metal atoms, it occurs with a low activation barrier of only 0.41 eV and should be the fastest process of all. The optimum path involves diffusion between the more stable hollow sites, with the hcp 40 sites being energetically favored by 0.28 eV over the fcc sites, across less stable bridge sites labeled b. In bulk cobalt, carbon atoms prefer energetically the octahedral interstitial sites over the tetrahedral sites. The optimum bulk path, presented in Fig. 4.4(b), involves diffusion normal to the surface between octahedral (o) sites across barriers at the triangular face centers (fc) of the octahedra. I emphasized one triangular face of an octahedron by the white dotted line in the middle panel of Fig. 4.4(b). In this view, the barrier fc site in the center of the triangle separates the favored o sites directly below and above. Since the Co atoms are closely packed in the hcp structure, passing through the center of the triangular face requires displacing atoms, which requires a high activation energy of 3.19 eV. This value is to be considered an upper limit, since presence of defects including vacancies should reduce the activation barrier for bulk diffusion significantly [75]. In contrast to a single crystal, the atomic packing at grain boundaries is less compact. Consequently, interstitial carbon atoms may find an energetically less costly diffusion path along the grain boundary than in the perfect lattice. A possible grain boundary structure that ends in a step edge is shown in the middle and bottom panel of Fig. 4.4(c). The atomic packing in this grain boundary resembles that of a simple cubic lattice, with interstitial carbon favoring energetically the body center bc sites in the cube center. The optimum diffusion path requires passing through a square face center fc at the interface of neighboring cubes. As seen in the top panel of Fig. 4.4(c), the activation barrier for the diffusion along this grain boundary is ≈1.3 eV, less than half the single crystal value. Considering growth conditions similar to those in Ref. [54], diffusion to the surface along this grain boundary should be ≈4×1010 times faster than in the perfect crystal at T = 900 K according to Arrhenius law. 41 (a) (b) Top view [0001] . Side view [0001] Figure 4.5 Possible formation mechanism of the cellular carbon nanofoam, represented by structural snap shots in top and side view. Different grains are distinguished by color and shading. Initial formation of a graphene nanoribbon along a step edge in (a) is followed by lateral growth of honeycomb cells in (b). Reproduced from [3], c 2012 American Physical Society. With the information at hand about the diffusion rates of the carbon feedstock, I proceed to discuss a possible growth scenario. The Co structure in Fig. 4.5 schematically depicts three grains, distinguished by color and shading. It is plausible to assume that the terrace height at both sides of the grain boundary may not be the same, yielding a step structure, which is best visible in side view. Under growth conditions [54] near 600◦ C, the fastest rate of carbon diffusion to the surface is along the grain boundary towards the step edge, where carbon may aggregate to a narrow graphene nanoribbon. Since according to my studies a zigzag edge binds more strongly to Co than an armchair edge, I consider a zigzag graphene nanoribbon attached to the step edge, as seen in Fig. 4.5(a). To best conform to the substrate, the nanoribbon acquires a washboard structure, depicted in the top panel in Fig. 4.5(a). The more reactive nanoribbon atoms, which protrude towards the terrace, are more likely to form bonds with carbon atoms diffusing along the terrace, thus initiating the formation of foam cells. In the meantime, atoms or flakes diffusing along the upper terrace become the feedstock for the termination of the foam layer by caps, as seen in Fig. 4.5(b). I hope that 42 this information may encourage follow-up experimental studies aiming at synthesizing the carbon foam and related carbon allotropes. 4.5 Electronic Structure After identifying the interesting structural properties of carbon foam, I continue to investigate its electronic and transport properties. As DFT usually underestimate the fundamental band gap, electronic band structure results obtained by DFT must be interpreted very carefully. Also in carbon foam, the DFT band gaps become too small and even turn negative in the bulk system, placing the bottom of the conduction band below the top of the valence band. While this is not a point of central interest in my study, I should at least note that a better description of the band gaps may be obtained at substantial computational cost by performing calculations using the HSE06 hybrid functional [76, 37] or self-energy calculations using the GW approach. My HSE band structure results, presented in Fig. 4.6(f), show a fundamental A–H band gap of 0.5 eV. The size of the band gap is consistent with results of my GW calculations using the Vienna Ab initio Simulation Package (VASP) [28, 29, 30]. Comparison of my DFT and HSE results in Figs. 4.6(e) and 4.6(f) shows that the bulk carbon foam, according to HSE, is a semiconductor, and that the band gap reduction in DFT is an artifact of that approach. Even considering such artifact, electronic structure of the valence and the conduction band region is believed to closely represent experimental results. My DFT calculations of the band structure of the hydrogen-covered thin carbon slab with two sp3 surfaces, shown in Fig. 4.1(b), are presented in Fig. 4.7(a). These results suggest this system to be semiconducting, same as its bulk counterpart [50]. On the other hand, the Hcovered, sp2 -terminated slab, depicted in Fig. 4.1(c), is clearly metallic according to the band 43 (a) (b) (c) (d) (e) (f) 1-honeycomb thick slab DFT 3-honeycomb thick slab DFT 5-honeycomb thick slab DFT 7-honeycomb thick slab DFT Bulk foam DFT Bulk foam HSE Figure 4.6 Electronic band structure of sp3 -terminated foam slabs and of bulk carbon foam. DFT-based band structure is presented for (a) a 1-honeycomb thick, (b) a 3-honeycomb thick, (c) a 5-honeycomb thick, (d) a 7-honeycomb thick slab, and (e) for the bulk structure. (f) Band structure of bulk carbon foam based on the HSE hybrid functional. Reproduced from [2], c 2014 American Physical Society. structure results in Fig. 4.7(b). This result is surprising, since conduction in semiconducting carbon structures including diamond has so far only been observed in presence of unsaturated dangling bonds [77]. My band structure results in Fig. 4.7(b) also show a Dirac-like cone similar to graphene at the K-point in the Brillouin zone. Moderate n-doping should be able to align it with the Fermi level, providing carriers with the same desirable properties as graphene. To learn more about the character of the new conduction states, I display in Fig. 4.7(c) the charge density associated with states close to EF , shown by the shaded region in Fig. 4.7(b). These states with p character, which are oriented within the surface plane of the foam, are located only on the sp2 sublattice. This is very different from graphene, where conduction is caused by nearest-neighbor hopping between p⊥ orbitals oriented normal to the surface, which are equally occupied at all lattice sites. To judge the suitability of carbon foam for electronic applications, I calculated quantum conductance of hydrogen-covered thin sp2 - and sp3 -terminated carbon foam slabs, shown in 44 Energy (eV) (a) sp3 terminated slab sp2 terminated slab (c)) EF Near-top view (b) sp2 terminated slab (d) G/G0 Energy (eV) Left lead EF Side view 6 5 4 3 2 1 0 -4 -3 Scattering region Right lead sp2 termination sp3 termination -2 -1 0 1 2 Energy (eV) 3 4 Figure 4.7 Electronic structure of thin foam slabs. DFT-based band structure of a thin, hydrogen-covered foam slab with (a) sp3 and (b) sp2 termination on both sides. (c) Charge distribution in the sp2 -terminated slab corresponding to states in the energy range EF − 2 eV≤E≤EF + 2 eV, indicated by shading in (b). (d) Conductance G of sp3 - and sp2 terminated slabs along the armchair direction, in units of the conduction quantum G0 . Reproduced from [2], c 2014 American Physical Society. Fig. 4.1(b-c), and present the results in Fig. 4.7(d). These results reflect my band structure results, namely a large conductance at small bias values in the sp2 -terminated carbon foam slab and a conduction gap of 0.6 eV in the sp3 -terminated slab. Additional results, also for thicker slabs, suggest that conduction is linked to the sp2 -termination of the carbon foam surface and nearly isotropic. To confirm the generality of this finding and obtain insight into its origin, I used the linear combination of atomic orbitals (LCAO) technique [78, 79, 80] to study the electronic 45 (b) (c) (g) 1..+ 2. neighbors (h) M 1. neighbors only ī pŒ K Vppʌ(1) EF EF (i) (e) EF pŒ Vpp(ʌ+ı)(2) sp2 terminated slab (d) Vppʌ(1) (j) (f) bulk (k) sp3 termination Ep-Vppʌ(1) sp2 termination 1..+ 2. neighbors 1. neighbors only Ep EF Ep+Vppʌ(1) sp3 site sp2 site Vppʌ(1) Vpp(ʌ+ı)(2) normal to surface (a) bulk Figure 4.8 Interpretation of electronic structure results for the sp2 -terminated slab (a-c) and bulk foam (d-f) using the LCAO technique. Band structure results considering only nearest-neighbor interaction (a,d), results of calculations that also include second-neighbor interaction between p states (b,e), and the corresponding densities of states (c,f). (g) Highsymmetry points in the slab and bulk Brillouin zones. (h) Rehybridized orbitals used in the calculation. p orbitals, shown in darker (green) shade, are responsible for conduction. (i) Tilted top view of the second neighbor-interaction between p states that dominate band dispersion near EF . Connectivity diagram for the bulk (j) and thin foam slabs (k) that helps explain the semiconducting behavior of bulk and sp3 -terminated foam, and the origin of the conducting state at the sp2 -terminated surface. Reproduced from [2], c 2014 American Physical Society. structure of carbon foam slabs. My results for the sp2 terminated slab of interest are presented in Figs. 4.8(a-c). The band associated with conduction, obtained by considering only nearest-neighbor interactions and presented in Fig. 4.8(a), was found to be quite different from its DFT-based counterpart in Fig. 4.7(b). To find out if this difference is caused by omitting second and third neighbor interactions, I introduced these interactions in my initial Hamiltonian. By selectively modifying individual hopping parameters, I have found that (i) the most significant changes in the electronic structure, including band broadening near EF , 46 are caused by second-neighbor Vppσ (2) and Vppπ (2) interactions; (ii) Vssσ (2) and Vspσ (2) interactions between second neighbors do not affect the electronic structure near EF and can be neglected; (iii) third-neighbor interactions contribute very little to the electronic structure and can also be safely neglected. Results obtained using the Hamiltonian augmented by second-neighbor interactions between p states, presented in Fig. 4.8(b), agree much better with the DFT results of Fig. 4.7(b), in particular regarding the width of the occupied part of the conduction band. The effect of the second-neighbor interaction is even more pronounced in the electronic structure of the bulk foam, shown in Figs. 4.8(d-f). As seen from the comparison between Figs. 4.8(d) and 4.8(e), neglecting the second-neighbor interaction in Fig. 4.8(d) caused a drastic narrowing of the bulk valence band and a significant increase in the fundamental band gap. LCAO calculations for the bulk with both first and second neighbor interactions, presented in Fig. 4.8(e), indicate an indirect A − H fundamental band gap of Eg = 1.12 eV, in qualitative agreement with non-self-consistent DFTB results. Having shown that the electronic structure near the Fermi level is well reproduced by the LCAO Hamiltonian, which also considers second-neighbor interactions between p states, I proceed to identify the reason for the fundamental difference between the sp2 - terminated metallic and sp3 - terminated semiconducting slab. I use the rehybridization theory [81, 82] to identify proper hybrid orbitals, since the atomic arrangements in the foam are generally not purely linear, hexagonal or tetrahedral. Hybrid orbitals |hi with i = 1, 2, 3, 4, which are associated with an atom, are linear combinations of the |s and a |p atomic orbital at this site. The direction of the |p orbital is not that of the Cartesian axes, but rather taken as that of nearest-neighbor bonds to up to three neighbors. The relative weights of the atomic orbitals are constructed [83] by enforcing the orthogonality condition hi |hj = δij . This 47 orthogonality condition is also used to construct any remaining hybrid orbitals. Selected hybrids at fourfold coordinated sp3 sites and at threefold coordinated sp2 sites are shown in Fig. 4.8(h). Analysis of the foam eigenstates at EF indicates dominance of p hybrids, highlighted by the darker color Fig. 4.8(h), at sp2 sites only. While parallel to the slab surface, these p orbitals are not aligned along a spatially fixed direction, but rather are locally normal to the graphitic strips lining the foam cells. These p hybrids on the sp2 sublattice play the role of π orbitals in the H¨ uckel Hamiltonian that describes most of the interesting physics in this system. I will show below that considering only first- and second-neighbor interaction in the H¨ uckel Hamiltonian is sufficient to explain, why bulk foam and sp3 -terminated slabs are semiconducting, whereas sp2 -terminated slabs become metallic. The spatial distribution of the p -dominated states at EF at one of the slab surfaces, shown in Fig. 4.8(i), agrees well with the DFT results presented in Fig. 4.7(c). In a H¨ uckel Hamiltonian with up to second-neighbor interactions, the p hybrids may interact in two ways only. The first type of interaction is normal to the surface and involves pairs of adjacent sp2 sites, which interact by the first-neighbor Vppπ (1) interaction. This is shown schematically by the red solid lines in the right panel of Fig. 4.8(h) and in Figs. 4.8(j-k), which represent the connectivity diagram of the foam. The second type of interaction is in-plane and much weaker, involving only the second-neighbor Vppπ (2) and Vppσ (2) interactions between p states. This is shown schematically by the green dashed lines in the right panel of Fig. 4.8(i) and in Figs. 4.8(j-k). In contrast to graphene, where the relatively strong Vppπ (1) interaction forms a network that connects all atoms in the layer, this interaction connects only pairs of adjacent sp2 sites in the carbon foam, as illustrated in Figs. 4.8(j-k). I conclude that bulk foam states near EF resemble those of a set of carbon dimers, 48 where the Vppπ (1) interaction splits the Ep energy eigenvalue into Ep ±Vppπ (1), indicated by the blue dash-dotted lines in Figs. 4.8(d-f). This corresponds to opening up a large fundamental band gap with Eg = 2|Vppπ (1)|, which turns carbon foam to a semiconductor. The δ-function-like localized states at Ep ±Vppπ (1), originating from decoupled dimers, are clearly visible in the density of states of the system with nearest-neighbor interactions only, shown by the red dashed line in Fig. 4.8(f). Much weaker second-neighbor Vppπ (2) and Vppσ (2) interactions couple these dimers and broaden the localized states into wide valence and conduction bands, shown by the solid black lines in the density of states in Fig. 4.8(f). Now it is rather straight-forward to explain the fundamental difference between the bulk, sp3 - and sp2 -terminated surfaces in terms of conductivity. Deciding whether a structure is metallic or semiconducting boils down to the simple question, whether all of the sp2 sites are connected as first-neighbors to another sp2 site. If so, then all energy eigenvalues Ep of the p states will split by Vppπ (1) into the eigenvalue pair Ep ±Vppπ (1), which opens up a gap. This is the case for the bulk and an sp3 -terminated surface of the foam. If, on the other hand, there are at least some sp2 sites present with no first-neighbor bonds to other sp2 sites, then the energy eigenvalue Ep of these p states will not split. Presence of partly filled p states at Ep = EF indicates that such a system should be conducting. The second-neighbor interaction between p states provides a weak coupling between sp2 site pairs or lone sp2 sites at the surface, as shown schematically in Figs. 4.8(j) and 4.8(k). In the bulk or in the sp3 -terminated foam, this second-neighbor interaction broadens the pair of sharp eigenvalues into a valence and conduction band that are still separated by a band gap, keeping their semiconducting character. At the sp2 -terminated surface, the weak second-neighbor interaction between p states at sp2 sites with no sp2 nearest neighbors will broaden the state at Ep to a metallic band, as seen also in the density of states presented in 49 Fig. 4.8(c). This behavior is not restricted to carbon foam, but should also occur at surfaces of isomorphic foams of Si, Si-carbides or BN systems. This reasoning also explains, why conductance depends only on the presence of sp2 sites with a particular topological arrangement. I may thus conclude that the conducting state at sp2 -terminated surfaces is topologically protected and is independent of the fact that the surface has been represented by a finite-thickness slab. This finding is furthermore confirmed by my electronic structure and conductance results for a thicker slab. With increasing thickness, the electronic structure of finite slabs approaches that of the bulk material. My conclusions regarding occurrence of states at EF are valid not only for perfect sp2 terminated surfaces, but also for all defective structures that contain sp2 sites with no sp2 nearest neighbors. Isolated monatomic vacancies in the bulk will produce corresponding midgap defect states. Lithographic patterning of the sp3 -terminated semiconducting surface by selective removal of lines of surface atoms will yield lines with sp2 -termination that will behave as conductive nanowires. 4.6 In-plane conduction anisotropy at hydrogen-covered foam surfaces Results of my quantum transport calculations along different directions on hydrogen-covered sp3 - and sp2 -terminated 1-honeycomb thick foam slabs are shown in Fig. 4.9. To judge the level of anisotropy, the infinite slab can be thought to be subdivided into strips of constant width, aligned with the transport direction. The fact that conductance results along the armchair and the zigzag direction are nearly the same strongly suggests that the conductance along the foam surface is nearly isotropic. 50 (a) (b) 1-honeycomb-thick slab sp2 armchair zigzag G/G0 sp2 Energy (eV) (c) (d) sp3 1-honeycomb-thick slab armchair zigzag G/G0 sp3 Energy (eV) Figure 4.9 Tilted view of (a) an sp2 - and (c) an sp3 -terminated foam slab, indicating the armchair and the zigzag direction along the surface. Conduction along the armchair and the zigzag direction in (b) an sp2 - and (d) an sp3 -terminated 1-honeycomb thick foam slab. Reproduced from [2], c 2014 American Physical Society. 4.7 Conclusions In conclusion, I studied the formation and structural as well as thermal stability of cellular foam-like carbon nanostructures by performing ab initio density functional calculations. I found that these systems with a mixed sp2 /sp3 bonding character may be compressed by reducing the symmetry of the honeycomb cells. The foam may accommodate the same type of defects as graphene, and its surface may be be stabilized by terminating caps. I postulate that the foam may form under non-equilibrium conditions near grain boundaries 51 on a carbon-saturated Co surface and should be thermally stable up to ≈3, 700 K. Regarding to electronic properties of this novel carbon foam structure, I have identified theoretically an unusual conduction mechanism occurring at the sp2 -terminated surface of a semiconducting carbon foam. To obtain microscopic insight into the origin of this mechanism, I augmented ab initio electronic structure and quantum conductance calculations by rehybridization theory calculations. I found that the occurrence of new conduction states in this system is intimately linked to the topology of the surface and not limited to foams of elemental carbon. My interpretation based on rehybridization theory indicates that conduction in the foam derives from first- and second-neighbor interactions between p orbitals lying in the surface plane, which are related to p⊥ orbitals of graphene. The topologically protected conducting state occurs on bare and hydrogen-terminated foam surfaces and is thus unrelated to dangling bonds. My results for carbon foam indicate that the conductance behavior may be further significantly modified by surface patterning, allowing to create conductive paths at the surface of the semiconducting foam matrix. 52 Chapter 5 Semiconducting bulk black phosphorus and 2D phosphorene The following discussion represents the theoretical part of the publication by Han Liu, Adam T. Neal, Zhen Zhu, Zhe Luo, Xianfan Xu, David Tom´anek and Peide Ye, Phosphorene: an unexplored 2D semiconductor with a high hole mobility, ACS Nano 8, 4033-41 (2014) [4]. This study is a collaborative project between our Theory group and Prof. Peide Ye’s Experimental group at Purdue University that synthesized and characterized the system. 5.1 Introduction Phosphorene, 2D counterpart of layered black phosphorus, which could be exfoliated mechanically, similar to other layered materials like graphene [18] and MoS2 [23]. Unlike metallic graphite and graphene, bulk black phosphorus is a direct band gap semiconductor with a band gap around 0.3 eV [84, 85], and similar semiconducting character is also maintained in thin slabs of phosphorene. Electronic devices of phosphorene shows both high on/off ratio and a high carrier mobility. These two advantages shed light on phosphorene as an emerging material for 2D electronics application, surpassing well established 2D materials like graphene and MoS2 . 53 (b) Energy (eV) (a) a1 EF a2 Bulk Brillouin zone U Z U X T Γ Y S Figure 5.1 Structural and electronic properties of layered black phosphorus. (a) Calculated band structure of bulk black phosphorus. (b) Atomic arrangement in one monolayer and the Brillouin zone of the bulk system. Reproduced from [4], c 2014 American Chemical Society. 5.2 Computational methods I determined the equilibrium structure, stability and electronic properties of bulk and multilayer structures of black phosphorus using ab initio density functional theory (DFT) calculations as implemented in the SIESTA code [27]. Multilayer structures are represented by a periodic array of slabs, separated by a 15 ˚ A thick vacuum region. I used the Perdew-BurkeErnzerhof [33] exchange-correlation functional, norm-conserving Troullier-Martins pseudopotentials [39], and a double−ζ basis including polarization orbitals. The reciprocal space was sampled by a fine grid [43] of 8×8×1 k-points in the Brillouin zone of the primitive unit cell. I used a mesh cutoff energy of 180 Ry to determine the self-consistent charge density, which provided us with a precision in total energy of ≤2 meV/atom. 5.3 Results and discussion My theoretical results of bulk black phosphorus are presented in Fig. 5.1. The calculated band structure in Figs. 5.1(a) indicates that bulk is a semiconductor. The monolayer dis54 (b) Energy (eV) (a) Top view EF Side view S S Figure 5.2 (a) Electronic structure of monolayer phosphorene. (b) Density of electrons ρvb near the top of the valence band of a P monolayer. States in the energy range EF − 1 eV < E < EF , indicated by shading in (a), are shown at the isosurface value ρ = 3×10−3 e/˚ A3 . Reproduced from [4], c 2014 American Chemical Society. playing a significantly larger fundamental band gap, shown in Fig. 5.2, about 0.9 eV and I expect the experimental value of the band gap is even larger. The bulk unit cell contains two AB stacked layers and the structure of one layer is shown in Fig. 5.1(b). The calculated bulk lattice parameters a1 = 3.36 ˚ A, a2 = 4.53 ˚ A, and a3 = 11.17 ˚ A suggest that the in-layer structure is very close to that of the monolayer, which is characterized by orthogonal lattice parameters a1 = 3.35 ˚ A and a2 = 4.62 ˚ A. The optimized geometry agrees well with the observed layered structure. Even though the in-layer morphology is close to that of graphene, the layers are puckered. The relatively large value of a3 is caused by layer puckering and a relatively large inter-layer separation, close to the graphite value. Electronic band structure results obtained by DFT must be interpreted very carefully. Even though DFT generally underestimates the fundamental band gap, the electronic structure of the valence and the conduction band region is believed to closely represent experimental results. Therefore, I expect the charge density associated with states near the top of the valence band, shown in Fig. 5.2(b), to be represented accurately. These states and their 55 (b) Eg (eV) Eg (eV) (a) ˰[ ˰\ Indirect Direct gap gap 1Oൺ෱ % Number of layers Nl % % % 6WUDLQ˰ Figure 5.3 Band gap engineering of few-layer phosphorene through (a) the number of layers and (b) applied strain along the x- and y-direction of a monolayer. Change from a direct to an indirect band gap is indicated by the dashed line in (b). Reproduced from [4], c 2014 American Chemical Society. hybrids with electronic states of the contact electrodes will play a crucial role in the carrier injection and quantum transport. My calculated bulk band gap value Eg = 0.04 eV agrees with published results [86, 87], but is significantly smaller than the experimental value [84] of 0.31 eV. More interesting than its precise value is its dependence on the number of layers Nl in a multi-layer slab, shown in Fig. 5.3(a). My finding that Eg may vary between 0.9 eV in a monolayer and 0.04 eV in the bulk promises the possibility of tuning the electronic properties of this system. Equally interesting is the sensitive dependence of the gap on in-layer stress along different directions, shown in Fig. 5.3(b). Of particular importance is my finding that even a moderate ≈3% in-plane compression, possibly caused by epitaxy on a substrate, turns a direct-gap to an indirect-gap semiconductor with a significantly smaller gap. 56 5.4 Conclusions To conclude, I have revealed the puckered structure of black phosphorus and its 2D counterpart phosphorene. I also found that few-layer phosphorene is a direct band gap semiconductor with a much larger band gap than bulk black phosphorus. The band gap value depends sensitively on the number of layers, which turn the number of layers into an effective method to engineer the electronic properties of few-layer phosphorene. Moreover, appling moderate in-layer tensile strain could modify the band gap value and induce a direct-to-indirect band gap transition. 57 Chapter 6 Semiconducting layered blue phosphorus The following discussion is closely related to a publication by Zhen Zhu, and David Tom´anek, Semiconducting layered blue phosphorus: A computational study, Phys. Rev. Lett. 112, 176802 (2014) [5]. 6.1 Introduction Elemental phosphorus is stable in a large number of structures, including the common white, red, violet and black allotropes [88, 89], with the color defined by the fundamental band gap. Most stable among them is black phosphorus, which – besides graphitic carbon – is the only layered structure of an elemental solid I know of. Individual layers of black phosphorus, shown in Fig. 6.1(a), resemble the honeycomb structure of graphene in terms of connectivity, but are non-planar. I noted that specific dislocations may convert black phosphorus, characterized by armchair ridges in the side view of the layers, to a well-defined structure with zigzag puckering, shown in Fig. 6.1(b). Assuming that the modified structure is stable, it appears worthwhile to study the equilibrium atomic arrangement in the bulk and the possibility of exfoliating individual layers. Whereas the observed fundamental band gap of 0.3-0.4 eV in black phosphorus [84, 85] is rather narrow, it is intriguing to see if 58 (a) Black phosphorus Top view (b) Blue phosphorus Top view a2 a2 a1 a1 Side view Side view (d) (c) Blue osphorus phosphorus A d=5.63 Å Side view B A Black phosphorus Figure 6.1 The layered structure of (a) black and (b) blue phosphorus in top and side view. Atoms at the top and bottom of the non-planar layers are distinguished by color and shading and the Wigner-Seitz cells are shown by the shaded region. (c) Schematic of the conversion of black to blue phosphorus by dislocations, highlighted by the shaded regions and arrows. (d) Equilibrium structure of AB stacked blue phosphorus in side view. Reproduced from [5], c 2014 American Physical Society. the modified phosphorus structure is a semiconductor with a wider gap. If also the carrier mobility were high, few-layer phosphorus in the new phase would become a worthy contender in the emerging field of post-graphene 2D electronics. Here I use ab initio calculations to investigate this previously unknown phase of phosphorus that shares its layered structure with the most stable black phosphorus allotrope. I find this structural phase, which I call ‘blue phosphorus’, to be nearly as stable as black phosphorus. Whereas the in-plane hexagonal structure and bulk layer stacking of blue phosphorus 59 are closely related to graphite, the main advantage of blue phosphorus is its wide fundamental band gap in excess of 2 eV. Due to the weak inter-layer interaction, blue phosphorus should exfoliate easily to form quasi-2D structures for potential electronic applications. I study a likely transformation pathway from black to blue phosphorus and discuss possible ways to synthesize the postulated structure. 6.2 Computational techniques My computational approach to gain insight into the equilibrium structure, stability and electronic properties of blue phosphorus is based on ab initio density functional theory (DFT) as implemented in the SIESTA [27] and VASP [28] codes. I used periodic boundary conditions throughout the study, with multilayer structures represented by a periodic array of slabs separated by a 15 ˚ A thick vacuum region. I used the Perdew-Burke-Ernzerhof [33] exchange-correlation functional, norm-conserving Troullier-Martins pseudopotentials [39], and a double-ζ basis including polarization orbitals. The reciprocal space was sampled by a fine grid [43] of 8×8×1 k-points in the Brillouin zone of the primitive unit cell. I used a mesh cutoff energy of 180 Ry to determine the self-consistent charge density, which provided us with a precision in total energy of < 2 meV/atom. All geometries have been optimized by SIESTA using the conjugate gradient method [90], until none of the residual Hellmann-Feynman forces exceeded 10−2 eV/˚ A. My SIESTA results for the optimized geometry, interlayer interactions and electronic structure were found to be in general agreement with VASP calculations. To verify the stability of the system at elevated temperatures, I performed a canonical molecular dynamics calculation of the monolayer and free-standing flakes using the SIESTA code. I used the Verlet integration algorithm to cover time periods 60 up to 1 ps with the time step of 1 fs. 6.3 Structural results The optimized reference structure of black phosphorus, shown in Fig. 6.1(a), is presented next to the proposed blue phosphorus structure in Fig. 6.1(b). The top view of both structures illustrates their similarity with the honeycomb lattice of graphite, which contains two atoms per layer per unit cell. Both phosphorus allotropes differ from graphite in the non-planar structure of their layers. In top view, the isotropic structure of blue phosphorus in Fig. 6.1(b) differs significantly from the anisotropic structure of black phosphorus in Fig. 6.1(a). As seen in side view in Figs. 6.1(a) and 6.1(b), the puckered zigzag structure in the cross-section of blue phosphorus differs from the distinct armchair ridges that cause the anisotropy of black phosphorus. The puckering in the blue phosphorus monolayer is similar to the postulated structure of single-wall phosphorus nanotubes [91]. The structural relationship between blue and black phosphorus is illustrated schematically in Fig. 6.1(c). A dislocation is introduced in a monolayer of black phosphorus by flipping specific P atoms from a ‘down’ to an ‘up’ position without changing the local bond angles, as described below and indicated by the arrows in Fig. 6.1(c). Subjecting every fourth row of P atoms to this transformation converts a monolayer of black to blue phosphorus. The location of dislocation lines in the monolayer structures is emphasized by the shaded regions in Fig. 6.1(c). Atoms in the layers of blue phosphorus are covalently bonded at the equilibrium distance of 2.27 ˚ A, resulting in a large binding energy of 3.28 eV/atom for the monolayer, comparable with experimental value [61] for bulk red phosphorus. This value differs from that of black 61 phosphorus by less than 2 meV/atom, suggesting that blue and black phosphorus are equally stable A weak inter-layer interaction of 6 meV/atom holds the layered structure together at the inter-layer distance d = 5.63 ˚ A, as seen in Fig. 6.1(d), which indicates the possibility of easy exfoliation. The AB hexagonal stacking and the ABC rhombohedral stacking of layers differ energetically by less than 1 meV/atom. The optimized hexagonal unit cell of an isolated blue phosphorus monolayer, shown in Fig. 6.1(c), is spanned by lattice vectors a1 and a2 , with a = |a1 | = |a2 | = 3.33 ˚ A. The influence of the inter-layer interaction on the in-layer structure is small, causing only a negligible change from a = 3.324 ˚ A in the bulk to a = 3.326 ˚ A in the isolated monolayer. With its higher symmetry, the smaller hexagonal Wigner-Seitz cell of blue phosphorus, which contains two atoms and is shown in Fig. 6.1(b), differs from the rectangular Wigner-Seitz cell of the anisotropic black phosphorus monolayer with 4 atoms/unit cell, as shown in Fig. 6.1(a). Still, monolayers of blue and black phosphorus may form an ‘ideal’ in-layer connection by the dislocation illustrated in Fig. 6.1(c). 6.4 Vibrational structure results One way to compare the stability and structural rigidity of the different phosphorus allotropes is by studying the vibration spectrum. My results for the vibration spectra of blue and black phosphorus monolayers are presented in Fig. 6.2. I find the vibration spectra to be rather similar, reflecting a very similar bonding character. Acoustic and optical modes are well separated in both black and blue phosphorus. The harder longitudinal optical modes reflect a higher in-plane rigidity of the blue phosphorus allotrope in comparison to the accordion-like black phosphorus structure. The calculated flexural rigidity value D = 62 (a) Black phosphorus Y Wave number (cm-1) Œ ˁ Y S X (b) Blue phosphorus S X ˁ M Œ ˁ M K K ˁ Figure 6.2 Vibrational band structure ω(k) of a monolayer of (a) black and (b) blue phosphorus. The slope of the dashed lines along the longitudinal acoustic branches near Γ corresponds to the speed of sound and the in-plane stiffness. Reproduced from [5], c 2014 American Physical Society. 0.84 eV in blue phosphorus is lower than the D = 1.51 eV value in black phosphorus. These values provide a quantitative explanation for the dispersion of the flexural acoustic modes ZA near the Γ point [92] in the two allotropes. The high rigidity of a free-standing blue phosphorus monolayer was also confirmed by my molecular dynamics calculations at nonzero temperatures. I next compare the slopes of the longitudinal acoustic branches near Γ, which correspond to the speed of sound and reveal the in-plane stiffness. As seen in Fig. 6.2(a), the speed of sound along the Γ − Y direction in black phosphorus, vsΓ−Y = 7.8 km/s, is significantly higher than the vsΓ−X = 3.8 km/s value along the Γ − X direction, reflecting an anisotropy in the elastic constants. The lower rigidity along the Γ − X direction, corresponding to the a1 direction in Fig. 6.1(a), reflects the fact that compression along a1 requires primarily bond bending, which comes at a lower energy cost than bond stretching. In strong contrast to those findings, my results in Fig. 6.2(b) indicate that the in-plane elastic response of blue 63 (b) Energy (eV) (a) (c) Top view M Γ K Side view Γ Γ 0 (d) 2 Bulk 8765 4 3 2 1 1 (e)2.5 AB stacking ABC stacking Egap(eV) 2 Egap(eV) 1.8 2 Density of states Number of layers N 1.6 1.5 1.4 1 1.2 1 0 0.25 0.5 0.75 1 0.5 -5% 0% 5% Strain˰ 1/(Number of layers N) Figure 6.3 (a) Electronic band structure and (b) density of states of a blue phosphorus monolayer. (c) Electron density ρvb in the 0.1 eV wide energy range near the top top of the valence band in blue phosphorus, indicated by the green shaded region in (a) and (b). ρvb is represented at the isosurface value ρ = 7×10−3 e/˚ A3 and superposed with a ball-and-stick model of the structure. Dependence of the fundamental band gap Eg on (d) the number of layers and (e) the in-plane uniform stretching. Reproduced from [5], c 2014 American Physical Society. phosphorus is nearly isotropic, with nearly the same value vs = 7.7 km/s for the speed of sound along the Γ − M and the Γ − K direction. The predicted finite in-layer compressibility of these non-planar structures is advantageous when accommodating lattice mismatch during Chemical Vapor Deposition (CVD) growth on a substrate. 64 6.5 Electronic structure results My DFT results for the electronic structure of blue phosphorus are presented in Fig. 6.3. The calculated band structure and density of states, presented in Figs. 6.3(a) and 6.3(b), indicate that a blue phosphorus monolayer is a semiconductor with an indirect band gap Egap ≈2 eV. Comparison with more proper self-energy calculations based on the GW approach, performed by VASP, indicate that the DFT band gap is underestimated by ≈1.0±0.2 eV in mono- and multi-layers of blue P as a common shortcoming of DFT. Still, the electronic structure of the valence and the conduction band region in DFT is believed to closely represent experimental results. Therefore, I expect the charge density associated with states near the top of the valence band, shown in Fig. 6.3(c), to be represented accurately. These states correspond to the energy range highlighted by the green shading in Figs. 6.3(a) and 6.3(b), which extends from mid-gap to 0.1 eV below the top of the valence band. These states cause the inter-layer band dispersion in few-layer systems, and their hybrids with electronic states of the contact electrodes will play a crucial role in the carrier injection and quantum transport. More interesting than the precise value of the fundamental band gap is its dependence on the number of layers N in a multi-layer slab, shown in Fig. 6.3(d). This is a consequence of the inter-layer dispersion near the Fermi level in the bulk material. Independent of the type of stacking, I find Egap to be inversely proportional to N between one monolayer and the bulk structure. I conclude that modifying the slab thickness may change the value of Egap by up to a factor of 2, which is very important for electronic applications. The N -dependent band gap value of ≤3 eV lies just above the photon energy of visible blue light. I derive the name ‘blue phosphorus’ from this absorption edge, which plays a key role in the optical appearance. 65 As seen in Fig. 6.3(e), the band gap Egap depends also sensitively on the applied in-layer strain σ. In view of the non-planarity of the structure, the discussed strain range −10% < σ < +10% can be achieved without a large energy penalty, as discussed in connection with the vibration spectra. The possibility to change the band gap value by < 50% in strained epitaxial geometries on different substrates is one more indication that blue phosphorus may find intriguing applications in nanoelectronics. 6.6 Conversion from black to blue phosphorus I next studied the possibility of converting a monolayer of black to blue phosphorus by introducing an array of dislocations, depicted schematically in Fig. 6.1(c). My discussion of the conversion process including energy estimates is presented in Fig. 6.4. The starting black phosphorus, the final blue phosphorus, and an intermediate structure are depicted in side view in Fig. 6.4(a). To estimate changes in the atomic arrangement during the conversion process, I changed and constrained the out-of-plane displacement of specific atoms in the unit cell to follow the ‘black-to-blue’ trajectory and relaxed all other degrees of freedom. This provided us with a sequence of structures ‘1-8’ in-between black and blue phosphorus, which loosely define the reaction coordinate in Fig. 6.4(b). The change in the unit cell from black to blue phosphorus has been imposed between steps ‘1’ and ‘2’ by deforming it from the initial size and shape in black phosphorus to the final rectangular supercell in blue phosphorus, followed by atomic relaxation. My results for the relative energy with respect to the black phosphorus structure ‘1’ again illustrate my finding that the blue phosphorus structure ‘8’ is equally stable as ‘1’. The activation barrier for the conversion process is likely overestimated due to the constraints imposed on the intermediate structures and may further be lowered 66 (a) (b) 0.6 Side view Energy (eV) 0.5 0.4 0.3 0.2 0.47 eV/atom 0.1 0 -0.1 1 2 3 4 5 6 7 8 Reaction coordinate Figure 6.4 (a) Series of structural snapshots depicting the transition from black to blue phosphorus. (b) Total energy change during the transformation from black to blue phosphorus. The reaction coordinate range ‘1-2’ describes initial in-layer stretching along the horizontal direction, identified as a1 in Fig. 6.1(a), followed by out-of-plane displacement of atoms at a fixed value of a1 in the range ‘2-8’. Reproduced from [5], c 2014 American Physical Society. by stretching black phosphorus layers. Thus the true value should be below the already low value of 0.47 eV/atom, indicating the relative ease of mechanical conversion from black to blue phosphorus. As suggested earlier, the weak inter-layer interaction should allow mechanical exfoliation of blue phosphorus in analogy to the black allotrope [4]. As a matter of fact, depositing mechanically a monolayer of black phosphorus onto a stepped substrate may cause formation of dislocation lines and thus formation of narrow domains of blue phosphorus in the monolayer structure along the step edges of the substrate. In analogy to graphene [93, 94] and silicene [95], also layered phosphorus structures may be grown by CVD on specific substrates. As suggested in this study, blue and black phosphorus should be equally stable. The preferential phase should thus be determined by the lattice constant and the symmetry of the substrate in order to maximize the adsorption energy. Consequently, I expect blue phosphorus to form preferentially on substrates with hexagonal symmetry and a matching 67 lattice constant, such as MoS2 , or the (0001) surfaces of Zr and Sc, whereas black phosphorus should start growing on substrates with a rectangular lattice. Due to the low energy of forming a dislocation that connects blue and black phosphorus, both allotropes could coexist on particular substrates, including stepped surfaces, to optimize the adlayer-substrate interaction. One of the main reasons for the interest in 2D semiconductors including graphene for electronic applications is the observed high mobility of carriers. Related quasi-2D systems, including MoS2 , bring the benefit of a nonzero band gap, but display lower intrinsic carrier mobility due to enhanced electron-phonon coupling, primarily caused by the presence of heavy elements such as Mo [96]. It appears likely that the blue phosphorus structure, similar to black phosphorus [4], may exhibit a higher carrier mobility than MoS2 . The combination of a significant band gap and high carrier mobility would turn blue phosphorus into an excellent contender for a new generation of nano-electronic devices. Infinite lattice (a) t = 0 ps (b) T = 300 K t = 0.15 ps (c) T = 300 K t = 0.3 ps (d) T = 1500 K t = 0.5 ps (e) T = 1500 K t = 1.0 ps Top view Side view Finite flake (f) t = 0 ps (g) T = 300 K (h)T = 300 K t = 0.5 ps t = 1.0 ps (i) T = 1500 K t = 0.2 ps (j) T = 1500 K t = 0.7 ps Top view Side view Figure 6.5 Snap shots of canonical molecular dynamics simulations depicting structural changes in (a-e) a contiguous monolayer and (f-j) a free-standing, finite flake of blue phosphorus at different temperatures. The unit cell of the lattice and the flake contain 64 P atoms. The initial geometries are shown in (a) and (f). Reproduced from [5], c 2014 American Physical Society. 68 6.7 Stability at elevated temperatures To study the structural stability of a blue phosphorus monolayer, I also performed canonical molecular dynamics simulations of a periodic layer structure and a free-standing finite flake and present my results in Fig. 6.5. At room temperature, the infinite layer undergoes only minimal changes. The flake exhibits reconstruction at the edge, but otherwise keeps the structure in the middle of the layer intact. The higher temperature of 1500 K is well above the melting point of phosphorus (TM = 863 K for red phosphorus) [61]. My results indicate significant structural changes, which are somewhat suppressed in the infinite lattice structure. There is no indication of a preferential structural change to the atomic arrangement found in black phosphorus. 6.8 Conclusions In conclusion, I have conducted ab initio calculations to investigate a previously unknown layered phase of phosphorus, which I call ‘blue phosphorus’. I find blue phosphorus to be nearly as stable as black phosphorus, the most stable phosphorus allotrope. While sharing the atomic connectivity with the honeycomb lattice of graphene, layers of blue phosphorus are non-planar. Whereas the bulk layer stacking of blue phosphorus is closely related to graphite, the main advantage of this allotrope is its wide fundamental band gap in excess of 2 eV. Due to the weak inter-layer interaction, blue phosphorus should exfoliate easily to form quasi-2D structures for potential electronic applications. I have investigated a likely transformation pathway from black to blue phosphorus and show that the postulated structure may form spontaneously by CVD on a lattice-matched substrate or may be result by stretching black phosphorus. Monolayers of blue phosphorus may form a structurally ideal connection to 69 monolayers of black phosphorus with a different electronic structure. 70 Chapter 7 Metal-semiconductor transition in 2D few-layer grey arsenic The following discussion is closely related to a publication by Zhen Zhu, Jie Guan and David Tom´anek, Unusual electronic structure of few-layer grey arsenic: A computational study, Phys. Rev. B 91, 161404(R) (2015) [6]. 7.1 Introduction There is growing interest in two-dimensional semiconductors with a significant fundamental band gap and a high carrier mobility. Whereas obtaining a reproducible and robust band gap has turned into an unsurmountable obstacle for graphene [22, 21], the presence of heavy transition metal atoms in layered dichalcogenide compounds limits their carrier mobility [97]. Few-layer structures of layered phosphorus allotropes, such as black phosphorus, are rapidly attracting attention due to their combination of high mobility and significant band gaps [4, 5, 98, 99]. I find it conceivable that other isoelectronic systems, such as arsenic, may display similar structural and electronic properties as few-layer phosphorene while being chemically much less reactive [100]. In this respect, the most abundant grey arsenic allotrope is the structural counterpart of the layered A7 or blue phosphorus [5]. Arsenic is commonly known for its toxicity, which is highest for the yellow As allotrope and should not be of concern for 71 few-layer nanostructures. Whereas crystalline grey arsenic displays rhombohedral stacking of layers and is semimetallic [101, 102, 103, 104, 105, 106, 107], loss of crystallinity opens a fundamental band gap in the amorphous structure [108, 102]. Even though few-layer grey arsenic has not been studied yet, analogies with blue phosphorene make arsenic monolayers and bilayers plausible candidates for 2D semiconductors. Here I study the equilibrium geometry and electronic structure of thin films of layered grey arsenic using ab initio density functional theory. In contrast to bulk grey As that is semimetallic, thin films display a significant band gap that depends sensitively on the number of layers, in-layer strain, layer stacking and inter-layer spacing. I find that metallic character can be introduced by increasing the number of layers beyond two or by subjecting semiconducting monolayers and bilayers to moderate tensile strain. The strain-induced semiconductor-metal transition is triggered by changes in the band ordering near the top of the valence band that causes an abrupt change from σ to π character of the frontier states. 7.2 Computational techniques My computational approach to gain insight into the equilibrium structure, stability and electronic properties of arsenic structures is based on ab initio density functional theory as implemented in the SIESTA [27] and VASP [28, 29, 30, 31] codes. I use periodic boundary conditions throughout the study, with multilayer structures represented by a periodic array of slabs separated by a vacuum region > 15 ˚ A. Unless specified otherwise, I use the PerdewBurke-Ernzerhof (PBE) [33] exchange-correlation functional for most calculations. Selected results are compared to the Local Density Approximation (LDA) [32, 42] and other functionals including the optB86b-vdW functional [35, 36] that provide a better description of van der 72 (a) (b) Monolayer grey As Bulk grey As (ABC stacking) A d a2 B a1 C A Side view Top view Figure 7.1 (a) Side view of the rhombohedrally (ABC) stacked layered structure of bulk grey arsenic. (b) Top view of the buckled honeycomb structure of a grey arsenic monolayer. Atoms at the top and bottom of the non-planar layers are distinguished by color and shading and the Wigner-Seitz cell is shown by the shaded region. Reproduced from [6], c 2015 American Physical Society. Waals interactions and the HSE06 [76, 37] hybrid functional. In my SIESTA calculations I use norm-conserving Troullier-Martins pseudopotentials [39], and a double-ζ basis including polarization orbitals. The reciprocal space is sampled by a fine grid [43] of 16×16×1 k-points in the Brillouin zone of the primitive unit cell for 2D structures and 16×16×3 k-points for the bulk. I use a mesh cutoff energy of 180 Ry to determine the self-consistent charge density, which provides us with a precision in total energy of ≤2 meV/atom. All geometries have been optimized using the conjugate gradient method [90], until none of the residual Hellmann-Feynman forces exceeded 10−2 eV/˚ A. Equilibrium structures and energies based on SIESTA have been checked against values based on the VASP code. 7.3 Structural results In contrast to the AB-stacked isoelectronic black phosphorus, bulk grey arsenic prefers the rhombohedral (or ABC) layer stacking, with the optimized structure shown in Fig. 7.1. The monolayer of grey As, depicted in top view in Fig. 7.1(b), resembles the honeycomb lattice 73 of graphene with two atoms per unit cell. Unlike planar graphene, however, the unit cell is buckled, similar to blue (or A7) phosphorus [5], and very different from layered black phosphorus [4]. Interatomic interactions within a monolayer are covalent, resulting in a nearest-neighbor distance of 2.53 ˚ A. Observed and calculated structural and cohesive properties in grey arsenic are summarized in Table 7.1. The calculated inter-layer separation in the ABC-stacked bulk system is d = 3.58 ˚ A, similar to the observed value [109]. The interlayer interaction energy of ≈0.02 eV/atom, based on PBE, is slightly higher than in blue phosphorus [5]. While this value is likely underestimated, the optB86b value of 0.17 eV/atom and the LDA value of 0.16 eV/atom likely overestimate the interlayer interaction. The low interlayer interaction energy, similar to graphite and black phosphorus, suggests that few-layer As may be obtained by mechanical exfoliation from the bulk structure. The small difference in the length of the in-plane lattice vectors a = |a1 | = |a2 | = 3.64 ˚ A in the isolated monolayer and a = 3.85 ˚ A in the bulk structure can be traced back to a small difference in buckling of the layers that is introduced by the weak interlayer interaction. Changes in buckling are characterized by the pyramidalization angle [81]. In agreement with the experiment, I find AA-stacked grey arsenic to be less stable than the ABC-stacked structure, even though the energy difference lies within 10 meV/atom. The optimum inter-layer separation in the less favorable AA stacking increases to d = 5.15 ˚ A. The optimum geometry is usually rather independent of the DFT functional in most covalent solids, but this is often not the case in layered materials with a significant van der Waals interlayer interaction. In Figs. 7.2(a) and 7.2(b) I show the dependence of the inter-layer spacing d and the in-layer lattice constants a = a1 = a2 on the number of layers N for different types of stacking and different DFT functionals. With the bilayer being an 74 exception, the trends in d(N ) and a(N ) are independent of the DFT functional. (PBE) AA (LDA) AA (optB86b-vdW) 0.1 (b) AA ABC ABC ABC 3.9 a1=a2 (Å) d (Å) 6 4 AA ABC ABC ABC (PBE) AA (LDA) AA (optB86b-vdW) 3.8 5 4 Expt 3.7 0.07 (c) N=2 0.05 0.06 AB (PBE) AB (LDA) ˂(H9DWRP (a) ˂(H9DWRP 7 0 -0.05 3.6 3 1 2 3 4 5 6 Number of layers N ∞ 3.5 1 2 3 4 5 6 Number of layers N ∞ -0.1 3 4 d (Å) 5 6 N=1 (PBE) 0.04 0.03 0.02 0.01 Expt (d) 0.05 0 tensile compressive -5% -4% -3% -2% -1% 0 1% 2% 3% 4% 5% Strain Figure 7.2 (a) The inter-layer spacing d and (b) the in-layer lattice constant a1 = a2 of grey arsenic slabs as a function of the number of layers N . Presented results, obtained using different DFT functionals, consider both the energetically favorable ABC and the less stable AA stacking. (c) Energy dependence of an AB stacked As bilayer on the interlayer spacing d as determined by PBE and LDA. Equilibrium geometries are indicated by arrows and ΔE = 0 refers to two isolated grey arsenic monolayers. (d) Energy dependence of a monolayer on in-layer strain, based on PBE. The dashed line represents harmonic behavior. Reproduced from [6], c 2015 American Physical Society. I find that the Local Density Approximation (LDA) [32] as well as the van der Waalscorrected optB86b-vdW DFT functional [35, 36] underestimate the interlayer spacing with respect to the observed bulk value, whereas the Perdew-Burke-Ernzerhof (PBE) [33] functional overestimates this value. Deviations of LDA and PBE results from the observed inter-layer distance are expected, since none of these functionals is designed to take into account van der Waals interactions. Much less expected is the error in the inter-layer spacing value predicted by the optB86b-vdW functional, which is supposed to treat van der Waals interactions more accurately. In my study of the isoelectronic black phosphorus [98], I also found that optB86b-vdW overestimates the inter-layer interaction energy when compared to more precise Quantum Monte Carlo (QMC) calculations. In view of this uncertainty, I present results based on PBE as the default functional in this Chapter. My results in Fig. 7.2(a) indicate a very weak dependence of d on the number of layers. The inter-layer separations are rather uniform within a given slab, but somewhat larger near the surface, especially in ABC stacked slabs. As expected intuitively, the inter-layer 75 spacing in the energetically less favorable AA stacking is generally larger than in the observed rhombohedral or ABC stacking. While the inter-layer interaction may be small, I find that it still changes the in-layer structure, primarily by changing the buckling of the layers. I characterized the degree of non-planarity or buckling by the pyramidalization angle θP , defined in Fig. 7.3(a). As seen in Fig. 7.3(b), increasing the number of layers from a monolayer to the bulk system decreases the pyramidalization angle, but does not change the interatomic bond length noticeably. It is this change of the pyramidalization angle that is chiefly responsible for the dependence of the lattice constants a1 = a2 on the number of layers N , shown in Fig. 7.2(b). The observed decrease of the pyramidalization angle with increasing N explains the decrease in the interlayer distance d with increasing N , seen in Fig. 7.2(a), and a corresponding increase of the inter-layer interaction. 35° (a) (b) (c) N=1 34° 33° θP θσπ 32° 31° 30° ∞ 2 Number of layers N 1 0% 1% 2% 3% 4% 5% Strain Figure 7.3 (a) Definition of the pyramidalization angle θp in a buckled monolayer of arsenic with sp3 hybridization. The σ orbitals extend along the interatomic bonds, and the direction of the π orbital is indicated by the dashed line. Dependence of the pyramidalization angle on (b) the number of layers in a thin film and (c) in-layer tensile strain in a monolayer, as obtained using the PBE, LDA and optB86b-vdW functionals. The dashed lines in (c) serve as guides to the eye. Reproduced from [6], c 2015 American Physical Society. 76 As seen in Fig. 7.2(b), the in-layer lattice constant shows a much stronger dependence on the number of layers in the favorable ABC stacking than in the less favorable AA stacking. The changes are significant, amounting to a 5% reduction in the in-layer lattice constant in the monolayer with respect to the bulk structure. As common in DFT calculations using nonlocal exchange-correlation functionals including PBE, the interlayer interaction is underestimated and the interlayer distance larger than in LDA. This contrast is particularly large in layered grey arsenic structures, as evidenced in the bilayer results in Figs. 7.2(a) and 7.2(c). As mentioned above, the PBE value (in contrast to the LDA value) of the interlayer distance in the AB stacked bilayer does not follow the trend of d(N ) for N > 2. To make sure that this is not an artifact of the optimization, I present in Fig. 7.2(c) the energy of the bilayer as a function of the interlayer separation d. Both LDA and PBE indicate the presence of a well-defined single optimum structure. In as similar way, I find that the deformation energy ΔE of a grey arsenic monolayer subject to in-layer strain, presented in Fig. 7.2(d), shows no deviation from the expected near-parabolic behavior, suggesting that the buckled structure represents a single optimum geometry. As seen in Fig. 7.3(c), also the pyramidalization angle in a stretched monolayer decreases uniformly with tensile strain, with no indications of a structural bistability. 7.4 Electronic structure results In agreement with observations [102], my DFT results indicate that bulk grey arsenic is semimetallic. My corresponding DFT results for the electronic structure of a monolayer of grey arsenic are presented in Fig. 7.4. In stark contrast to the bulk, the monolayer structure 77 (b) Monolayer grey As (c) Top view -1 ˁ0.ˁ DOS σ Side view 1 EF Top p view -1 -2 σ ˁ0.ˁ Side view Energy (eV) EF Energy (eV) Energy (eV) 1 -2 -5% strain 2 2 +5% strain (d) 2 2.0 1 1.5 tensile 1.0 EF Top view -1 -2 compressive Eg (eV) (a) π Side view ˁ0.ˁ Monolayer 0.5 0 -10% -5% Bilayer 0 5% 6WUDLQ 10% Figure 7.4 Electronic structure of (a) a relaxed monolayer of grey As, and the same monolayer subject to a uniform in-layer strain of (b) -5% and (c) +5%. The energy range between the Fermi level EF and 0.2 eV below the top of the valence band is green shaded in the band structure in the left panels. Electron density ρvb associated with these states, superposed with a ball-and-stick model of the structure, is shown in the right panels. (d) Dependence of the fundamental band gap Eg on the in-layer strain in a monolayer and an AB-stacked bilayer of grey As. Reproduced from [6], c 2015 American Physical Society. is semiconducting with an indirect fundamental band gap Eg ≈1.71 eV. Comparison with more precise HSE06 [76, 37] hybrid functional calculations, indicates that the PBE value of the band gap is likely underestimated by > 0.4 eV in few-layer grey arsenic as a common shortcoming of DFT. Still, the electronic structure of the valence and the conduction band region in DFT is believed to closely represent experimental results. Therefore, I expect the charge density associated with frontier states near the top of the valence band, shown in Fig. 7.4(a), to be represented accurately. These states correspond to the energy range highlighted by the green shading in the band structure plots, which extends from mid-gap to 0.2 eV below the top of the valence band. The E(k) plot in Fig. 7.4(a) shows that states near the top of the valence band at Γ display a strong dispersion. This is a signature of a very low hole mass, caused by frontier states forming a network of σ-orbitals connecting neighboring atoms, as seen in the right panel of Fig. 7.4(a). This is quite different from black phosphorene, where the frontier valence band states are dominated by atomic out-of-plane p-orbitals with little overlap, which reduces the band dispersion and thus increases the hole mass. I find that the effective mass near Γ in 78 -1 d ˁ0.ˁ DOS π Side view 1.5 1 Top view EF -1 -2 ˁ0.ˁ 2.0 d σ Side view Top view EF ABC ABC AA AA (PBE) (LDA) (PBE) (LDA) 1.0 0.5 -1 -2 Eg (eV) Energy (eV) Energy (eV) EF (d) 2 1 Top view +5% strain (c) 2 1 -2 -5% strain (b) Bilayer grey As 2 Energy (eV) (a) π ˁ0.ˁ d 0 1 2 3 4 5 6 Side view ∞ 1XPEHURIOD\HUV1 Figure 7.5 Electronic structure of (a) a relaxed AB-stacked bilayer of grey As, and the same bilayer subject to a uniform in-layer strain of (b) -5% and (c) +5%. The energy range between the Fermi level EF and 0.2 eV below the top of the valence band is green shaded in the band structure in the left panels. Electron density ρvb associated with these states, superposed with a ball-and-stick model of the structure, is shown in the right panels. The inter-layer spacing d in the side views is reduced for clarity. (d) Dependence of the fundamental band gap Eg on the number of layers N and the DFT functional in few-layer As with ABC and AA stacking. The lines in (d) are guides to the eye. Reproduced from [6], c 2015 American Physical Society. few-layer arsenic is not only lower, but – in contrast to black phosphorene [110, 4] – also isotropic. Since higher hole mobility values have been reported in bulk grey arsenic than bulk black phosphorus [102], I believe that also arsenic monolayers and bilayers should display a higher hole mobility than monolayers and bilayers of black phosphorus. As seen in Fig. 7.4(b), a uniform 5% in-layer compression reduces the band gap, but keeps it indirect and does not change drastically the character of the frontier states. According to Fig. 7.4(d), compression in excess of 10% would close the band gap, turning the monolayer metallic. Interestingly, also a uniform in-layer stretch reduces the band gap in the monolayer significantly. As seen in Fig. 7.4(c), stretching the monolayer by 5% moves the bottom of the conduction band near Γ. The less dispersive valence band with an energy eigenvalue −1.3 eV at Γ in the relaxed monolayer, shown in Fig. 7.4(a), moves up and becomes the top valence band in the stretched layer in Fig. 7.4(c). This band gradually flattens near Γ upon stretching the monolayer, and the system becomes a direct-gap semiconductor at tensile strain values ≥+6%. 79 More important is the change in the character of the frontier states caused by the motion of this band relative to the other valence bands near Γ. In contrast to the relaxed and compressed layer depicted in Figs. 7.4(a) and 7.4(b), the frontier states in the stretched layer are lone-pair p-orbitals of tetrahedral As that are normal to the layer. The lower overlap between these atomic p-orbitals is also the cause of the low dispersion of the corresponding π-like band near Γ. As I will see later, the change in character of the frontier states from σto π-like has a profound effect on the electronic structure near EF in multi-layer systems. The profound effect of the change in character of the frontier states can be seen most easily when comparing the fundamental band gap Eg in a monolayer and a bilayer in Fig. 7.4(d). Eg is essentially the same in both systems during compression, since the frontier states have σ-character and show little overlap between layers. The character change of the frontier states to π-like upon stretching increases significantly the inter-layer overlap between these states, causing a drastic reduction in Eg of the bilayer as compared to the monolayer. Quite significant in this respect is my finding that a bilayer should turn metallic for tensile strains 5%. This critical tensile strain value is likely underestimated in my PBE study and is projected to increase to < 7% according to my HSE06 calculations. I note that similar tensile strain values have been achieved experimentally on suspended graphene membranes that are more resilient to stretching due to their planar geometry [111, 112, 113]. A closer look at the behavior of frontier states in a bilayer subject to different levels of strain is offered in Fig. 7.5. As seen in Fig. 7.5(a), the top valence band in the relaxed bilayer has a similar low dispersion in the electronic structure near Γ as a stretched monolayer in Fig. 7.4(c). Also the right panels of the two sub-figures confirm the similar π character of the frontier states in a stretched monolayer and a relaxed bilayer. Since the bottom of the conduction band in the relaxed bilayer is not at Γ, the relaxed As bilayer is an indirect-gap 80 semiconductor with a narrower gap than the relaxed monolayer. Upon compression, the band ordering near the top of the valence band changes in the bilayer. As seen in Fig. 7.5(b), the frontier states in a system subject to a 5% uniform compression change their character to σ-like, similar to a compressed monolayer. The reduced overlap between frontier states on neighboring layers increases the inter-layer separation d. A compressed arsenic bilayer remains an indirect-gap semiconductor. As seen in the right panel of Fig. 7.5(c), a uniform +5% tensile strain in the bilayer changes the character of the frontier states to π-like, similar to my findings for the stretched monolayer reported in Fig. 7.4(c). The frontier states, which are primarily distributed in the inter-layer region, bring the two layers closer together. The increased hybridization of inter-layer frontier states at a decreased interlayer distance causes a significant band gap reduction over the monolayer subject to the same level of tensile strain. As mentioned above, the bilayer should turn semi-metallic at moderate tensile strain values. The overall dependence of the fundamental band gap on the number of layers N and the stacking sequence, depicted in Fig. 7.5(d), shows a uniform trend of band gap reduction with increasing value of N , which has been noted also for the different layered allotropes of phosphorus [4, 5]. In particular, I find ABC-stacked arsenic slabs with N > 2 to be metallic. I also observe notable changes in the optimum inter-layer distance d with changing number of layers and stacking sequence, which strongly affect the electronic structure. I find the optimum inter-layer distance in AA-stacked structures to be much larger than in ABCstacked structures, which slows down the reduction of the band gap with growing N . To check the validity of this trend, I reproduced the band gap values obtained using both PBE and LDA exchange-correlation functionals for structures optimized by PBE in Fig. 7.5(d). It is well known that the fundamental band gap is usually underestimated in DFT cal81 (a) EF (b) Strained arsenic bilayer (7%) Arsenic monolayer PBE LDA HSE06 PBE HSE06 EF ˁ0.ˁ ˁ0.ˁ Figure 7.6 Electronic band structure of (a) a grey arsenic monolayer and (b) bilayer under 7% strain, obtained by DFT calculations using the LDA, PBE and HSE06 functionals. Reproduced from [6], c 2015 American Physical Society. culations. Even though the proper way to determine the electronic band structure involves advanced methods such as GW or QMC, I have not used these approaches, as they are computationally very demanding. I rather made use of the HSE06 hybrid functional [76, 37], which is believed to reproduce the fundamental band gap correctly. I compare the electronic band structure of a relaxed grey As monolayer, obtained using LDA, PBE and HSE06, in Fig. 7.6. For the sake of fair comparison, I use the identical geometry, namely that of a PBE-optimized monolayer, in the three sets of results obtained using the VASP code [28, 29, 30, 31]. My results indicate that both LDA and PBE give a nearly identical band structure with a fundamental band gap near 1.6 eV, in close agreement with SIESTA results [27]. HSE06 opens up the band gap to about 2.0 eV by essentially rigidly shifting the conduction band up with respect to the valence band. Results based on HSE06 suggest that also a relaxed bilayer is a semiconductor. According to HSE06, relaxed grey arsenic structures with more than two layers are either semimetallic or metallic. To better estimate the strain required to induce a semiconductor-metal transition in a 82 bilayer of grey arsenic, I performed HSE06 hybrid functional calculations and compared them to PBE results. Both PBE and HSE06 results indicate that the fundamental band gap decreases upon stretching the monolayer and bilayer. Similar to the monolayer, I find the HSE06 functional to increase the band gap value by essentially rigidly shifting the conduction states up, thus delaying the closure of the fundamental gap in the strained system. As seen in Fig. 7.6(b), the fundamental band gap closes as the arsenic bilayer is stretched by 7% according to HSE06, very close to the 5% value based on PBE. This small numerical difference between vastly different functionals indicates that a metal-insulator transition in a stretched bilayer should be observable at moderate strain values. The weak inter-layer interaction in layered grey arsenic should allow for a mechanical exfoliation of monolayers and bilayers in analogy to graphene and phosphorene. More appealing for large-scale production is the reported synthesis of thin films of grey arsenic by Molecular Beam Epitaxy (MBE) two decades ago [114], which should be also capable of producing monolayers and bilayers. Chemical Vapor Deposition (CVD), which had been used successfully in the past to grow graphene [93, 94] and silicene [95], may become ultimately the most common approach to grow few-layer grey arsenic on specific substrates. Substrates such as Ag(111), or even Zr(0001) and Hf(0001) should be advantageous to minimize the lattice mismatch during MBE or CVD growth. From the viewpoint of electronic applications, an ideal 2D semiconductor should combine a sizeable fundamental band gap with a high carrier mobility and chemical stability. Equally important is identifying a suitable way to make electrically transparent contacts. Graphene, with the exception of its vanishing band gap, satisfies the latter three criteria ideally. Transition metal dichalcogenides, including MoS2 , bring the benefit of a nonzero band gap, but display lower intrinsic carrier mobility due to enhanced electron-phonon coupling, primarily 83 caused by the presence of heavy elements such as Mo, and suffer from high tunneling barriers at contacts to the chalcogen atoms. Few-layer systems of black phosphorus and arsenic do show a sizeable band gap and the ability to form transparent contacts to metal leads. The higher carrier mobility in bulk grey As and black P over bulk MoS2 [102] has also been observed in ultrathin black phosphorus films [4], and the same is expected for monolayers and bilayers of arsenic as well. Of the two group V elements, the heavier arsenic appears more resilient to oxidation [100]. If indeed grey arsenic monolayers and bilayers turn out to be chemically stable, arsenic may become an excellent contender for a new generation of 2D nano-electronic devices. 7.5 Conclusions In conclusion, I have used ab initio density functional theory to study the equilibrium geometry and electronic structure of thin films of layered grey arsenic. I found the PBE, LDA and optB86b-vdW DFT functionals as well as the hybrid HSE06 functional to predict consistent trends for the interlayer distance and interaction as well as the fundamental band gap and character of the frontier states in monolayers, bilayers and few-layer systems of grey arsenic. In contrast to bulk grey As that is semimetallic, thin films display a significant band gap that depends sensitively on the number of layers, in-layer strain, layer stacking and inter-layer spacing. I find that metallic character can be introduced by increasing the number of layers beyond two or by subjecting semiconducting monolayers and bilayers to moderate tensile strain. The strain-induced semiconductor-metal transition is triggered by changes in the band ordering near the top of the valence band that causes an abrupt change from σ to π character of the frontier states. Due to the weak inter-layer interaction, grey arsenic should 84 exfoliate easily to form few-layer structures. Alternative ways to synthesize few-layer arsenic include MBE and CVD. 85 Table 7.1 Observed and calculated properties of layered grey arsenic. a = |a1 | = |a2 | is the in-plane lattice constant and d is the interlayer separation, as defined in Fig. 7.1. Ecoh is the cohesive energy and Eil is the interlayer interaction energy per atom. Structure Bulk(ABC) Bulk(ABC) (expt.) (theory) a ˚ a ( A) 3.76 3.85b – 3.85c – 3.82d d (˚ A) 3.52a 3.58b – 3.46c – 3.47d Ecoh 2.96e 2.86b (eV/atom) – 3.60c Eil – 0.02b (eV/atom) – 0.16c – 0.17d a Bulk(AA) (theory) 3.65b 3.64c 3.62d 5.15b 4.20c 4.31d 2.85b 3.53c 0.01b 0.10c 0.13d Monolayer (theory) 3.64b 3.61c 3.58d – – – 2.84b 3.45c – – – Experimental data of Ref. [109]. Results based on the DFT-PBE functional [33]. c Results based on the LDA [32]. d Results based on the optB86b van der Waals functional [35, 36]. e Experimental data of Ref. [61]. b 86 Chapter 8 Insight into structural phase coexistence of group V elements by a tiling model The following discussion is closely related to a publication by Jie Guan, Zhen Zhu, and David Tom´anek, Tiling phosphorene, ACS Nano 12, 12763-12768 (2014) [7]. 8.1 Introduction Phosphorene, a monolayer of black phosphorus, is emerging as a viable contender in the field of two-dimensional (2D) electronic materials [85, 84, 4]. In comparison to the widely discussed semi-metallic graphene, phosphorene displays a significant band gap while still maintaining a high carrier mobility [99, 4, 115, 116]. The flexible structure of semiconducting phosphorene [110, 117] is advantageous in applications including gas sensing [118], thermoelectrics [119], and Li-ion batteries [119]. Unlike flat sp2 -bonded graphene monolayers, the structure of sp3 -bonded phosphorene is buckled. There is a large number of sp3 -bonded layered phosphorene structures, including blue-P, γ-P, and δ-P [5, 98], which are nearly as stable as the related black phosphorene structure but exhibit very different electronic properties. I believe that the above list of stable phosphorene structures is still incomplete, giving 87 ˠ-P (N=1) Blue-P (N=0) (a) Black-P (N=2) (b) (c) (e) (f) (h) h) (i) i) Side view (d) Top view (g) g) Tiling model top bottom Figure 8.1 Atomic structure of different phosphorene allotropes in (a-c) side view, (d-f) top view, and (g-i) in a tiling model representation. Blue-P (a,d,g), γ-P (b,e,h), and black-P (c,f,i) can be distinguished by the structural index N . Primitive unit cells are emphasized by shading and delimited by black dashed lines in (d-i). Atoms at the top and bottom of the layer, as well as the corresponding tiles, are distinguished by color. Reproduced from [7], c 2014 American Chemical Society. rise to an unprecedented richness in terms of polymorphs and their electronic structure. Here I introduce a scheme to categorize the structure of different layered phosphorene allotropes by mapping the non-planar 3D structure of three-fold coordinated P atoms onto a two-color 2D triangular tiling pattern. In the buckled structure of a phosphorene monolayer, I assign atoms in “top” positions to dark tiles and atoms in “bottom” positions to light tiles. Optimum sp3 bonding is maintained throughout the structure when each triangular tile is surrounded by the same number N of like-colored tiles, with 0≤N ≤2. My ab initio density functional calculations indicate that both the relative stability and electronic properties depend primarily on the structural index N . Common characteristics of allotropes 88 with identical N suggest the usefulness of the structural index for categorization. The proposed mapping approach may also be applied to phosphorene structures with non-hexagonal rings, counterparts of planar haeckelite [120, 121], to point and line defects [122], and to 2D quasicrystals with no translational symmetry, which I predict to be nearly as stable as the hexagonal network. 8.2 Computational methods My computational approach to gain insight into the equilibrium structure, stability and electronic properties of various phosphorene structures is based on ab initio density functional theory (DFT) as implemented in the SIESTA [27]. I used periodic boundary conditions throughout the study. I used the Perdew-Burke-Ernzerhof (PBE) [33] exchange-correlation functional, norm-conserving Troullier-Martins pseudopotentials [39], and a double-ζ basis including polarization orbitals. Selected PBE results were compared to results based on the Local Density Approximation (LDA) [32, 42]. The reciprocal space was sampled by a fine grid [43] of 8×8×1 k-points in the Brillouin zone of the primitive unit cell. I used a mesh cutoff energy of 180 Ry to determine the self-consistent charge density, which provided us with a precision in total energy of ≤ 2 meV/atom. All geometries have been optimized by SIESTA using the conjugate gradient method [90], until none of the residual HellmannFeynman forces exceeded 10−2 eV/˚ A. 8.3 Tiling pattern of phosphorene The non-planar atomic structure of selected sp3 -bonded phosphorene allotropes is depicted in side and top view in Fig. 8.1(a-f). I find it convenient to map the 3D structure of 89 a phosphorene monolayer with threefold coordinated atoms onto a 2D tiling pattern by assigning a triangular tile to each atom, as shown in Fig. 8.1(g-i). There is a one-to-one correspondence between structures and tiling patterns, so that different structures can be distinguished by different tiling patterns. Dark-colored tiles are associated with atoms at the top and light-colored tiles with atoms at the bottom of the layer. Since each atom has 3 neighbors, each triangular tile is surrounded by 3 neighboring tiles, N of which have the same color. It is obvious that 0≤N ≤2 provides the atom associated with the central tile with a tetrahedral neighbor coordination associated with the favorable sp3 bonding. In my tiling model, N = 3 would represent the planar structure of an energetically unfavorable sp2 bonded lattice that, according to my findings, would spontaneously convert to a non-planar sp3 -bonded allotrope. As I will show in the following, different allotropes with N = 0, N = 1 and N = 2 share similar characteristics. Therefore, the structural index N is useful for primary categorization of the allotropes. In each structure depicted in Fig. 8.1, N maintains an identical value throughout the lattice, keeping the favorable sp3 bonding at all sites. I believe that this is the underlying reason for my finding that these structures are nearly equally stable [5, 98]. In the first category characterized by N = 0, all neighbors of a given atom have the same, but different height within the layer, as seen in Figs. 8.1(a) and 8.1(d). This translates into a tiling pattern, where all adjacent tiles have a different color, as seen in Fig. 8.1(g). There is only one structural realization within the N = 0 category, namely the blue-P allotrope. In the second category characterized by N = 1, each atom has one like neighbor at the same height and two unlike neighbors at a different height within the layer, as seen in Figs. 8.1(b) and 8.1(e). Besides the γ-P structure in Figs. 8.1(b) and 8.1(e), there is a θ-P allotrope, depicted Fig. 8.2(a), with the same structural index N = 1. The tiling patterns 90 (a) (b) ˥-P (N=1) ˡ-P (N=2) Atomic structure (c) Tiling model (d) a a b b Figure 8.2 Atomic structure of N = 1 and N = 2 phosphorene allotropes in top view (a-b) and the corresponding tiling model representation (c-d). The N = 1 θ-P allotrope in (a,c) and the N = 2 allotrope in (b,d) are structurally different from the allotropes with the same N in Fig. 8.1. Primitive unit cells are emphasized by shading and delimited by black dashed lines. Atoms at the top and bottom of the layer, as well as the corresponding tiles, are distinguished by color. The two different orientations of bonds between like atoms, indicated by the double arrows as guides to the eye, are denoted by “a” and “b”. Reproduced from [7], c 2014 American Chemical Society. of γ-P and θ-P, shown in Figs. 8.1(h) and 8.2(c), are characterized by a diamond harlequin pattern. Each diamond, formed of two adjacent like-colored triangles, is surrounded by unlike-colored diamonds. As a guide to the eye, I indicate the orientation of the diamonds, same as the direction of the atomic bonds, by the double arrows in Fig. 8.2(c). The shape of the primitive unit cells shown in Figs. 8.1 and 8.2 is chosen to see more easily the correspondence between the atomic structure and the tiling pattern. The primitive unit cell of γ-P contains 4 atoms according to Fig. 8.1(h) and that of θ-P contains 8 atoms, as seen 91 in Fig. 8.2(c). As indicated in Fig. 8.2(c), the orientation of diamonds in a row may be distinguished by the letters “a” or “b”. Whereas the perfect γ-P structure in Fig. 8.1(h) could be characterized by the sequence “aaaa . . .” and the structure of θ-P by the sequence “abab . . .”, an infinite number of different sequences including “abaa . . .” would result in an infinite number of N = 1 phosphorene structures. The most stable and best-known phosphorene allotrope is black-P, depicted in Fig. 8.1(c) and 8.1(f). Each atom in this structure has two like neighbors at the same height and one unlike neighbor at a different height, yielding a structural index N = 2. The tiling model of this structure type, shown in Fig. 8.1(i), contains contiguous arrays of like-colored diamonds. These arrays may be either straight, as in Fig. 8.1(i) for black-P, or not straight, as in Fig. 8.2(d) for the structurally different δ-P allotrope with the atomic structure shown in Fig. 8.2(b). Describing diamond orientation by letters “a” and “b” as in the case of N = 1, I may characterize black-P in Fig. 8.1(i) by the sequence “aaaa . . .” and δ-P in Fig. 8.2(d) by the sequence “abab . . .”. As in the case of N = 1, an infinite number of different sequences including “abaa . . .” would result in an infinite number of N = 2 phosphorene structures. The structural similarity and energetic near-degeneracy of N = 2 and N = 1 structures stems from the fact that a structural change from N = 2 to N = 1 involves only a horizontal shift of every other row, indicated by the horizontal lines in Figs. 8.2(c) and 8.2(d), by one tile. It is even possible to generate structural domains with different values of N . The energy cost of domain wall boundaries may be extremely low [98] if optimum sp3 bonding is maintained at the boundaries. As mentioned above, there is only one allotrope with N = 0, but infinitely many structures with N = 1 and N = 2. Of these, I identified and optimized all lattices with up to 28 atoms per primitive unit cell and selected other structures with up to 32 atoms per unit 92 (a) (b) 2.0 N=0 N=1 N=2 0.15 0.10 1.0 0.5 0.05 0 N=0 N=1 N=2 1.5 Eg (eV) ΔE (eV/atom) 0.20 024 0 8 12 16 20 24 28 32 Number of atoms/cell 024 8 12 16 20 24 28 32 Number of atoms/cell Figure 8.3 (a) Relative stability ΔE of different phosphorene allotropes with respect to black phosphorus and (b) their fundamental band gap Eg . The distribution of the ΔE and Eg values, provided in the right panels of the subfigures, indicates presence of three distinguishable groups that may be linked to the different values of the structural index N . The dashed and dotted lines are guides to the eye. Reproduced from [7], c 2014 American Chemical Society. cell. For each lattice, I identified the relative stability ΔE with respect to the most stable black phosphorene allotrope on a per-atom basis and plotted the values in Fig. 8.3(a). The electronic band structure of systems with large unit cells is very dense and hard to interpret in comparison to that of the allotropes discussed in Figs. 8.1 and 8.2. For each of these structures, though, I identified the value Eg of the fundamental band gap and provide the results in Fig. 8.3(b). I find that neither ΔE nor Eg display a general dependence on the size of the unit cell. I found all structures to be relatively stable. The small values ΔE < 0.15 eV/atom indicate a likely coexistence of different allotropes that would form under nonequilibrium conditions. All band gap values, which are typically underestimated in DFT-PBE calculations [27, 33], occur in the range between 0.3 eV and 2.0 eV, similar to the allotropes discussed in Figs. 8.1 and 8.2. Rather surprisingly, the distribution of ΔE and Eg values, shown in the right panels of the sub-figures Figs. 8.3(a) and 8.3(b), exhibits three peaks that can be associated with the structural index N , with a rather narrow variance caused by the differences between the 93 allotrope structures. I find the energetically near-degenerate blue phosphorene (N = 0 with 2 atoms per unit cell) and black phosphorene (N = 2 with 4 atoms per unit cell) structures to be the most stable, followed by other N = 2 structures with more than 4 atoms per unit cell. I found N = 1 structures to be the least stable of all. Similarly, the N = 0 blue phosphorene allotrope has the largest band gap, N = 2 allotropes have the smallest band gap, and N = 1 allotropes lie in between. The higher stability of N = 2 phosphorene structures in comparison to N = 1 allotropes indicates an energetic preference for phosphorus atoms forming zigzag chains at the same height rather than forming isolated dimers. Among the N = 2 structures, δ-P is the least stable, with ΔE≈0.07 eV/atom. All the other N = 2 structures fall in-between δ-P and black phosphorus in terms of stability. This finding is easy to understand, since all these structures are combination of black phosphorus and δ-P. For both N = 2 and N = 1 allotropes, I find structures with the same orientation of diamonds in the tiling pattern to be more stable. The γ-P structure, with all diamonds aligned in the same direction in the tiling pattern, is the most stable N = 1 phosphorene allotrope, but still less stable by 0.09 eV/atom than the N = 2 black phosphorene. At the other extreme of the relative stability range, θ-P with disordered diamond orientations in the tiling pattern is the least stable N = 1 allotrope, being 0.14 eV/atom less stable than black phosphorene. In analogy to what I concluded for N = 2 structures, all N = 1 phosphorene allotropes can be viewed as a combination of γ-P and θ-P, with their stability in-between the above limiting values. As mentioned above, also the distribution of Eg values, shown in the right panel of Fig. 8.3(b), indicates three distinct groups that can be associated with the structural index N . The largest band gap value of 2.0 eV in the only N = 0 structure, blue phosphorene, is 94 well separated from the band gap distribution of N = 1 and N = 2 structures that form a double-hump shape. I note that the two peaks in the band gap distribution of N = 1 and N = 2 allotropes are not as well separated as the two peaks in the stability distribution in Fig. 8.3(a), so the trends in the band gap value are not as clear as trends in the relative stability. In systems with large unit cells, band gaps of N = 1 structures are grouped around 0.8 eV, whereas band gaps of N = 2 structures are grouped around 0.5 eV. The largest spread in Eg values is in systems with very small unit cells. Among N = 1 allotropes, I find the smallest value Eg ≈0.5 eV in the structure with 4 atoms/unit cell (γ-P) and the largest value Eg ≈1.2 eV in the structure with 8 atoms/unit cell (θ-P). Band gap values of other N = 1 structures range between these two values. N = 2 structures have generally the lowest band gap values of the three groups. Among N = 2 systems, I find the largest value Eg ≈0.9 eV in the structure with 4 atoms/unit cell (black phosphorene) and Eg ≈0.3 eV in a system with 8 atoms/unit cell, the smallest gap value among several metastable structures of δ-P. Band gap values of other N = 2 structures range between these two values. As discussed earlier [5, 98], my PBE-based band gap values are generally underestimated. More precise quasiparticle calculations beyond DFT, including the GW formalism, indicate that the band gap values should be about 1 eV larger than the PBE values presented here [5, 123]. As the unit cell size of N = 1 and N = 2 structures grows infinitely large, I gradually approach amorphous phosphorene. Assuming that my findings in Fig. 8.3 are universal and not limited to the finite sizes addressed by my study, I conclude that the stability and the fundamental band gap of such amorphous structures should also be found in the range suggested by their structural index N . 95 (a) 4-8 (N=0) (b) (c) 3-12 (N=2) 5-7 (N=2) (d) 5-8 (N=2) Figure 8.4 Equilibrium structures (top) and corresponding tiling patterns (bottom) of 2D phosphorene with (a) 4-8, (b) 3-12, (c) 5-7, and (d) 5-8 membered rings. The color of the atoms and corresponding tiles distinguishes positions at the top and the bottom of the layer. Reproduced from [7], c 2014 American Chemical Society. 8.4 Defective structures The one-to-one mapping between 3D structures of periodic systems and 2D tiling patterns is not limited to a honeycomb lattice with 6-membered rings, but can equally well be applied to lattices with 3-, 4-, 5-, 7-, 8- and 12-membered rings found in planar haeckelites [120, 121]. The corresponding geometries and tiling patterns are shown in Fig. 8.4. Among these structures, 4-8 phosphorene has the highest symmetry, a relatively small unit cell with the shape of a square and a tiling pattern composed of right triangles. Besides the N = 0 structure depicted in Fig. 8.4(a), I can identify allotropes with 4-8 rings with structural indices N = 1 and N = 2. Other allotropes with 3-12, 5-7 and 5-8 rings, shown in Figs. 8.4(bd), may not exist in all the variants of the structural index N due to their lower symmetry. For example, the allotrope with 5-7 rings does not have a structure with N = 0. I find structures with non-hexagonal rings to be generally less stable than the most 96 stable black phosphorene, but the energy differences ΔE < 0.2 eV/atom are very small. Consequently, I expect that such structures should coexist with black phosphorene as either pure phases, or as local defects at domain wall boundaries, or as finite-size domains in the host layer. I find all phosphorene allotropes with non-hexagonal rings to be semiconducting, with the band gap determined primarily by the structural index N . 8.5 Quantum dot Phosphorene may also form aperiodic structures with no translational symmetry. Examples of such systems with only rotational symmetry are shown in Fig. 8.5. Fig. 8.5(a) depicts a phosphorene structure of type N = 2 with a C6v point group symmetry and the corresponding tiling pattern. In this structure, arrays of neighboring atoms form an alternating circular pattern about the center that can cover an infinite plane. The analogous N = 2 structure with C5v symmetry is depicted in Fig. 8.5(b), and analogous structures with Cnv symmetry could be imagined as well. To judge the stability of these aperiodic structures, I optimized finite-size flakes that were terminated by hydrogen atoms at the exposed edge. I found these structures to be semiconducting and as stable as the periodic structures discussed in Fig. 8.3(a), with ΔE = 0.07 eV/atom for the C6v and ΔE = 0.04 eV/atom for the C5v structure falling into range expected for N = 2. These findings indicate that my classification scheme and tiling model is useful to characterize monolayers of three-fold coordinated, sp3 -hybridized phosphorus atoms arranged in periodic or aperiodic patterns. Due to structural similarities between layered structures of group-V elements, I believe that my findings regarding relative stability, electronic structure and fundamental band gap will likely also apply to other systems including monolayers of 97 (a) (b) C6v (N=2) C5v (N=2) Atomic structure Tiling model Figure 8.5 Equilibrium structures (top panels) and corresponding tiling patterns (bottom panels) of phosphorene structures with (a) a C6v and (b) a C5v symmetry. Edges of the finitesize flakes are terminated by hydrogen atoms, colored in white. The color of the phosphorus atoms and the corresponding tiles distinguishes positions at the top or the bottom of the layer. Reproduced from [7], c 2014 American Chemical Society. arsenic, antimony and bismuth. Since the cohesive energy differences are rather small, I must consider the possibility that the stability ranking of the different allotropes at T = 0 and related properties [124] may depend on the DFT functional. I have compared PBE results for the relative stability of the different allotropes with LDA results and found the maximum difference in the relative stabilities of the different allotropes to be 0.02 eV/atom, which does not change the energy ranking of the allotropes. Since phosphorene structures will likely be synthesized at nonzero temperatures, the rel- 98 ative abundance of different allotropes will depend on their free energy at that temperature. Consequently, my total energy results for stability differences at T = 0 need to be corrected by also addressing differences in entropy at T > 0. Even though the decrease in free energy with increasing temperature should be similar in the different allotropes due to their similar vibration spectra [5, 98], minute differences in vibrational entropy may become important in view of the small differences between stabilities of the allotropes at T = 0, and could eventually change tfhe free energy ranking at high temperatures. 8.6 Conclusions In conclusion, I have introduced a scheme to categorize the structure of different layered phosphorene allotropes by mapping the non-planar 3D structure of three-fold coordinated P atoms onto a two-color 2D triangular tiling pattern, which could essentially be generalized to other group V elements which would share similar structural and bonding characters. In the buckled structure of a phosphorene monolayer, I assign atoms in “top” positions to dark tiles and atoms in “bottom” positions to light tiles. I found that optimum sp3 bonding is maintained throughout the structure when each triangular tile is surrounded by the same number N of like-colored tiles, with 0≤N ≤2. My ab initio density functional calculations indicate that both the relative stability and electronic properties depend primarily on the structural index N . Common characteristics of allotropes with identical N suggest the usefulness of the structural index for categorization. The proposed mapping approach may also be applied to phosphorene structures with non-hexagonal rings and to 2D quasicrystals with no translational symmetry, which I predict to be nearly as stable as the hexagonal network. 99 BIBLIOGRAPHY 100 BIBLIOGRAPHY [1] Z. Zhu, D.-G. Kwon, Y.-K. Kwon, and D. Tomanek, Enhancing mechanical toughness of aluminum surfaces by nano-boron implantation: An ab initio study, Chem. Phys. Lett. 620, 25 (2015). [2] Z. Zhu, Z. G. Fthenakis, J. Guan, and D. Tomanek, Topologically Protected Conduction State at Carbon Foam Surfaces: An ab initio Study, Phys. Rev. Lett. 112, 026803 (2014). [3] Z. Zhu and D. Tomanek, Formation and Stability of Cellular Carbon Foam Structures: An Ab Initio Study, Phys. Rev. Lett. 109, 135501 (2012). [4] H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tomanek, and P. D. Ye, Phosphorene: An Unexplored 2D Semiconductor with a High Hole Mobility, ACS Nano 8, 4033 (2014). [5] Z. Zhu and D. Tomanek, Semiconducting Layered Blue Phosphorus: A Computational Study, Phys. Rev. Lett. 112, 176802 (2014). [6] Z. Zhu, J. Guan, and D. Tomanek, Strain-induced metal-semiconductor transition in monolayers and bilayers of gray arsenic: A computational study, Phys. Rev. B 91, 161404 (2015). [7] J. Guan, Z. Zhu, and D. Tomanek, Tiling Phosphorene, ACS Nano 8, 12763 (2014). [8] W. S. Miller, L. Zhuang, J. Bottema, A. J. Wittebrood, P. De Smet, A. Haszler, and A. Vieregge, Recent Development in Aluminium Alloys for the Automotive Industry, Mat. Sci. Eng. A 280, 37 (2000). [9] M. Ashby, Y. Brechet, D. Cebon, and L. Salvo, Selection Strategies for Materials and Processes, Materials & Design 25, 51 (2004). [10] A. V. Smith and D. D. L. Chung, Titanium Diboride Particle-Reinforced Aluminium with High Wear Resistance, J. Mat. Sci. 31, 5961 (1996). [11] C. Srividya and S. V. Babu, Corrosion Resistance of Diamond-like Carbon-Coated Aluminum Films, Chem. Mat. 8, 2528 (1996). [12] M. Widom and M. Mihalkovic, Symmetry-Broken Crystal Structure of Elemental Boron at Low Temperature, Phys. Rev. B 77, 064113 (2008). [13] I. Boustani, Z. Zhu, and D. Tomanek, Search for the Largest Two-Dimensional Aggregates of Boron: An Ab Initio Study, Phys. Rev. B 83, 193405 (2011). [14] R. Becker, 2010, Rick Becker (private communication). [15] H. W. Kroto, J. R. Heath, S. C. O'Brien, R. F. Curl, and R. E. Smalley, C60: Buckminsterfullerene, Nature 318, 162 (1985). 101 [16] J. Zhang, Y. Feng, H. Ishiwata, Y. Miyata, R. Kitaura, J. E. P. Dahl, R. M. K. Carlson, H. Shinohara, and D. Tomanek, Synthesis and Transformation of Linear Adamantane Assemblies inside Carbon Nanotubes, ACS Nano 6, 8674 (2012). [17] S. Iijima, Helical microtubules of graphitic carbon, Nature 354, 56 (1991). [18] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsove, Electric Field Effect in Atomically Thin Carbon Films, Science 306, 666 (2004). [19] Y. Zhang, Y.-W. Tan, H. L. Stormer, and P. Kim, Experimental observation of the quantum Hall effect and Berry's phase in graphene, Nature 438, 201 (2005). [20] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81, 109 (2009). [21] D. C. Elias, R. R. Nair, T. M. G. Mohiuddin, S. V. Morozov, P. Blake, M. P. Halsall, A. C. Ferrari, D. W. Boukhvalov, M. I. Katsnelson, A. K. Geim, and K. S. Novoselov, Control of Graphene's Properties by Reversible Hydrogenation: Evidence for Graphane, Science 323, 610 (2009). [22] M. Y. Han, B. Ozyilmaz, Y. Zhang, and P. Kim, Energy Band-Gap Engineering of Graphene Nanoribbons, Phys. Rev. Lett. 98, 206805 (2007). [23] B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, and A. Kis, Single-layer MoS2 transistors, Nature Nanotech. 6, 147 (2011). [24] M. Born and R. Oppenheimer, Zur Quantentheorie der Molekeln, Annalen der Physik 389, 457 (1927). [25] P. Hohenberg and W. Kohn, Inhomogeneous Electron Gas, Phys. Rev. 136, B864 (1964). [26] W. Kohn and L. J. Sham, Self-Consistent Equations Including Exchange and Correlation Effects, Phys. Rev. 140, A1133 (1965). [27] E. Artacho, E. Anglada, O. Dieguez, J. D. Gale, A. Garcia, J. Junquera, R. M. Martin, P. Ordejon, J. M. Pruneda, D. Sanchez-Portal, and J. M. Soler, The SIESTA Method; Developments and Applicability, J. Phys. Cond. Mat. 20, 064208 (2008). [28] G. Kresse and J. Furthmuller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996). [29] G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B 47, 558 (1993). [30] G. Kresse and J. Hafner, Ab initio molecular-dynamics simulation of the liquid-metal amorphous-semiconductor transition in germanium, Phys. Rev. B 49, 14251 (1994). [31] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmentedwave method, Phys. Rev. B 59, 1758 (1999). 102 [32] D. M. Ceperley and B. J. Alder, Ground State of the Electron Gas by a Stochastic Method, Phys. Rev. Lett. 45, 566 (1980). [33] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77, 3865 (1996). [34] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys. 132, (2010). [35] J. Klimes, D. R. Bowler, and A. Michaelides, Chemical accuracy for the van der Waals density functional, J. Phys.: Cond. Matt. 22, 022201 (2010). [36] J. Klimes, D. R. Bowler, and A. Michaelides, Van der Waals density functionals applied to solids, Phys. Rev. B 83, 195131 (2011). [37] A. V. Krukau, O. A. Vydrov, A. F. Izmaylov, and G. E. Scuseria, Influence of the exchange screening parameter on the performance of screened hybrid functionals, J. Chem. Phys. 125, (2006). [38] D. Vanderbilt, Soft self-consistent pseudopotentials in a generalized eigenvalue formalism, Phys. Rev. B 41, 7892 (1990). [39] N. Troullier and J. L. Martins, Efficient pseudopotentials for plane-wave calculations, Phys. Rev. B 43, 1993 (1991). [40] G. B. Bachelet, D. R. Hamann, and M. Schluter, Pseudopotentials that work: From H to Pu, Phys. Rev. B 26, 4199 (1982). [41] P. E. Blochl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994). [42] J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Phys. Rev. B 23, 5048 (1981). [43] H. J. Monkhorst and J. D. Pack, Special points for Brillouin-zone integrations, Phys. Rev. B 13, 5188 (1976). [44] O. F. Sankey and D. J. Niklewski, Ab initio multicenter tight-binding model for molecular-dynamics simulations and other applications in covalent system, Phys. Rev. B 40, 3979 (1989). [45] E. Artacho, D. Sanchez-Portal, P. Ordejon, A. Garcia, and J. M. Soler, Linear-Scaling ab-initio Calculations for Large and Complex System, Phys. Stat. Sol. (b) 215, 809 (1999). [46] A. Zangwill, Physics at Surfaces, Cambridge University Press, Cambridge, U.K., 1988. [47] D. W. Jepsen, P. M. Marcus, and F. Jona, Low-Energy-Electron Diffraction from Several Surfaces of Aluminum, Phys. Rev. B 6, 3684 (1972), and Phys. Rev. B 8, (1973) 1786(E). 103 [48] R. Smoluchowski, Anisotropy of the Electronic Work Function of Metals, Phys. Rev. 60, 661 (1941). [49] J. Rottler, S. S. Schoenholz, and A. J. Liu, Predicting Plasticity with Soft Vibrational Modes: from Dislocations to Glasses, Phys. Rev. E 89, 042304 (2014). [50] A. Kuc and G. Seifert, Hexagon-preserving carbon foams: Properties of hypothetical carbon allotropes, Phys. Rev. B 74, 214104 (2006). [51] K. Umemoto, S. Saito, S. Berber, and D. Tomanek, Carbon foam: Spanning the phase space between graphite and diamond, Phys. Rev. B 64, 193409 (2001). [52] Z. Zhao, B. Xu, L.-M. Wang, X.-F. Zhou, J. He, Z. Liu, H.-T. Wang, and Y. Tian, Three Dimensional Carbon-Nanotube Polymers, ACS Nano 5, 7226 (2011). [53] A. N. Enyashin and A. L. Ivanovskii, Graphene allotropes, Phys. Stat. Sol. (b) 248, 1879 (2011). [54] J. A. Rodriguez-Manzo, C. Pham-Huu, and F. Banhart, Graphene Growth by a MetalCatalyzed Solid-State Transformation of Amorphous Carbon, ACS Nano 5, 1529 (2011). [55] J. A. Rodriguez-Manzo, M. Terrones, H. Terrones, H. W. Kroto, L. Sun, and F. Banhart, In situ nucleation of carbon nanotubes by the injection of carbon atoms into metal particles, Nature Nanotech. 2, 307 (2007). [56] M. Lin, J. P. Y. Tan, C. Boothroyd, K. P. Loh, E. S. Tok, and Y.-L. Foo, Dynamical Observation of Bamboo-like Carbon Nanotube Growth, Nano Lett. 7, 2234 (2007). [57] S. Helveg, C. Lopez-Cartes, J. Sehested, P. Hansen, B. Clausen, J. Rostrup-Nielsen, F. Abild-Pedersen, and J. Norskov, Atomic-scale imaging of carbon nanofibre growth, Nature 427, 426 (2004). [58] H. Wang, K. Sun, F. Tao, D. J. Stacchiola, and Y. H. Hu, 3D Honeycomb-Like Structured Graphene and Its High Efficiency as a Counter-Electrode Catalyst for DyeSensitized Solar Cells, Angew. Chem. Int. Ed. 52, 9210 (2013). [59] X. Jiang, J. Zhao, Y.-L. Li, and R. Ahuja, Tunable Assembly of sp3 Cross-Linked 3D Graphene Monoliths: A First-Principles Prediction, Adv. Funct. Mater. 23, 5846 (2013). [60] M. Brandbyge, J.-L. Mozos, P. Ordejon, J. Taylor, and K. Stokbro, Density-functional method for nonequilibrium electron transport, Phys. Rev. B 65, 165401 (2002). [61] C. Kittel, Introduction to Solid State Physics, Wiley, Hoboken, NJ, eighth edition, 2004. [62] S. Okada and A. Oshiyama, Nanometer-scale ferromagnet: Carbon nanotubes with finite length, J. Phys. Soc. Jpn. 72, 1510 (2003). [63] Y. Higuchi, K. Kusakabe, N. Suzuki, S. Tsuneyuki, J. Yamauchi, K. Akagi, and Y. Yoshimoto, Nanotube-based molecular magnets with spin-polarized edge states, J. Phys. Cond. Mat. 16, S5689 (2004). 104 [64] A. Mananes, F. Duque, A. Ayuela, M. J. Lopez, and J. A. Alonso, Half-metallic finite zigzag single-walled carbon nanotubes from first principles, Phys. Rev. B 78, 035432 (2008). [65] O. Hod and G. E. Scuseria, Half-Metallic Zigzag Carbon Nanotube Dots, ACS Nano 2, 2243 (2008). [66] Y.-H. Kim, J. Choi, K. J. Chang, and D. Tomanek, Defective fullerenes and nanotubes as molecular magnets: An ab initio study, Phys. Rev. B 68, 125420 (2003). [67] A. J. Stone and D. J. Wales, Theoretical studies of icosahedral C60 and some related species, Chem. Phys. Lett. 128, 501 (1986). [68] P. A. Thrower, The study of defects in graphite by transmission electron microscopy, in Chemistry and physics of carbon, Vol. 5, edited by P. L. Walker Jr., pages 217-320, Marcel Dekker, New York, 1969. [69] P. Y. Huang, C. S. Ruiz-Vargas, A. M. van der Zande, W. S. Whitney, M. P. Levendorf, J. W. Kevek, S. Garg, J. S. Alden, C. J. Hustedt, Y. Zhu, J. Park, P. L. McEuen, and D. A. Muller, Grains and grain boundaries in single-layer graphene atomic patchwork quilts, Nature 469, 389 (2011). [70] F. Banhart, J. Kotakoski, and A. V. Krasheninnikov, Structural Defects in Graphene, ACS Nano 5, 26 (2011). [71] J. C. Meyer, C. Kisielowski, R. Erni, M. D. Rossell, M. F. Crommie, and A. Zettl, Direct Imaging of Lattice Atoms and Topological Defects in Graphene Membranes, Nano Lett. 8, 3582 (2008). [72] J. Lahiri, Y. Lin, P. Bozkurt, I. I. Oleynik, and M. Batzill, An extended defect in graphene as a metallic wire, Nature Nano. 5, 326 (2010). [73] J. Tersoff, Empirical Interatomic Potential for Carbon, with Applications to Amorphous Carbon, Phys. Rev. Lett. 61, 2879 (1988). [74] J. Gao, J. Yip, J. Zhao, B. I. Yakobson, and F. Ding, Graphene Nucleation on Transition Metal Surface: Structure Transformation and Role of the Metal Step Edge, J. Am. Chem. Soc. 133, 5009 (2011). [75] B. Li, Q. Zhang, L. Chen, P. Cui, and X. Pan, Vacancy-mediated diffusion of carbon in cobalt and its influence on CO activation, Phys. Chem. Chem. Phys. 12, 7848 (2010). [76] J. Heyd, G. E. Scuseria, and M. Ernzerhof, Hybrid functionals based on a screened Coulomb potential, J. Chem. Phys. 118, 8207 (2003). [77] J. Sone and S. Okada, Massless electrons on hexagonal dangling bond network on hydrogen deposited diamond (111) and Si(111) surfaces, J. Phys. Soc. Jpn. 82, 064706 (2013). 105 [78] J. C. Slater and G. F. Koster, Simplified LCAO method for the periodic potential problem, Phys. Rev. 94, 1498 (1954). [79] D. Tomanek and M. A. Schluter, Growth regimes of carbon clusters, Phys. Rev. Lett. 67, 2331 (1991). [80] D. A. Papaconstantopoulos, editor, Handbook of the band structure of elemental solids, Plenum Press, New York and London, 1986. [81] R. C. Haddon, Hybridization and the orientation and alignment of π-orbitals in nonplanar conjugated organic molecules: π-orbital axis vector analysis (POAV2), J. Am. Chem. Soc. 108, 2837 (1986). [82] R. C. Haddon and L. T. Scott, π-orbital conjugation and rehybridization in bridged annulenes and deformed molecules in general: π-orbital axis vector analysis, Pure and Appl. Chem. 58, 137 (1986). [83] Z. G. Fthenakis, R. W. A. Havenith, M. Menon, and P. W. Fowler, Structural and electronic properties of the fullerene isomers of Si38: A systematic theoretical study, Phys. Rev. B 75, 155435 (2007). [84] Y. Maruyama, S. Suzuki, K. Kobayashi, and S. Tanuma, Synthesis and some properties of black phosphorus single crystals, Physica B+C 105, 99 (1981). [85] S. Narita, Y. Akahama, Y. Tsukiyama, K. Muro, S. Mori, S. Endo, M. Taniguchi, M. Seki, S. Suga, A. Mikuni, and H. Kanzaki, Electrical and optical properties of black phosphorus single crystals, Physica B+C 117&118, 422 (1983). [86] O. Prytz and E. Flage-Larsen, The influence of exact exchange corrections in van der Waals layered narrow bandgap black phosphorus, J. Phys.: Cond. Mat. 22, 015502 (2010). [87] Y. Du, C. Ouyang, S. Shi, and M. Lei, Ab initio studies on atomic and electronic structures of black phosphorus, J. Appl. Phys. 107, 093718 (2010). [88] R. Hultgren, N. S. Gingrich, and B. E. Warren, The Atomic Distribution in Red and Black Phosphorus and the Crystal Structure of Black Phosphorus, J. Chem. Phys. 3, 351 (1935). [89] H. Thurn and H. Kerbs, Crystal Structure of Violet Phosphorus, Angew. Chem. Int. Ed. Engl. 5, 1047 (1966). [90] M. R. Hestenes and E. Stiefel, Methods of Conjugate Gradients for Solving Linear Systems, J. Res. Natl. Bur. Stand. 49, 409 (1952). [91] G. Seifert and E. Hernandez, Theoretical prediction of phosphorus nanotubes, Chem. Phys. Lett. 318, 355 (2000). 106 [92] D. Sanchez-Portal, E. Artacho, J. M. Soler, A. Rubio, and P. Ordejon, Ab initio structural, elastic, and vibrational properties of carbon nanotubes, Phys. Rev. B 59, 12678 (1999). [93] K. S. Kim, Y. Zhao, H. Jang, S. Y. Lee, J. M. Kim, K. S. Kim, J.-H. Ahn, P. Kim, J.-Y. Choi, and B. H. Hong, Large-scale pattern growth of graphene films for stretchable transparent electrodes, Nature 457, 706 (2009). [94] A. Reina, X. Jia, J. Ho, D. Nezich, H. Son, V. Bulovic, M. S. Dresselhaus, and J. Kong, Large area, few-layer graphene films on arbitrary substrates by chemical vapor deposition, Nano Lett. 9, 30 (2009). [95] P. Vogt, P. De Padova, C. Quaresima, J. Avila, E. Frantzeskakis, M. C. Asensio, A. Resta, B. Ealet, and G. Le Lay, Silicene: Compelling Experimental Evidence for Graphenelike Two-Dimensional Silicon, Phys. Rev. Lett. 108, 155501 (2012). [96] M. M. Perera, M.-W. Lin, H.-J. Chuang, B. P. Chamlagain, C. Wang, X. Tan, M. M.-C. Cheng, D. Tomanek, and Z. Zhou, Improved carrier mobility in few-layer MoS2 fieldeffect transistors with ionic-liquid gating, ACS Nano 7, 4449 (2013). [97] M. S. Fuhrer and J. Hone, Measurement of mobility in dual-gated MoS2 transistors, Nature Nanotech. 8, 146 (2013). [98] J. Guan, Z. Zhu, and D. Tomanek, Phase Coexistence and Metal-Insulator Transition in Few-Layer Phosphorene: A Computational Study, Phys. Rev. Lett. 113, 046804 (2014). [99] L. Li, Y. Yu, G. J. Ye, Q. Ge, X. Ou, H. Wu, D. Feng, X. H. Chen, and Y. Zhang, Black phosphorus field-effect transistors, Nature Nanotech. 9, 373 (2014). [100] N. Burford, Y.-Y. Carpenter, E. Conrad, and C. D. L. Saunders, The Chemistry of Arsenic, Antimony and Bismuth, pages 1-17, John Wiley & Sons, Ltd, 2010. [101] D. Bullett, Density of states calculation for crystalline As and Sb, Solid State Commun. 17, 965 (1975). [102] O. Madelung, Semiconductors: data handbook, Springer, Berlin, third edition, 2004. [103] J. H. Xu, E. G. Wang, C. S. Ting, and W. P. Su, Tight-binding theory of the electronic structures for rhombohedral semimetals, Phys. Rev. B 48, 17271 (1993). [104] H. Tokailin, T. Takahashi, T. Sagawa, and K. Shindo, Electronic band structure of rhombohedral arsenic studied by highly-angle-resolved ultraviolet photoelectron spectroscopy, Phys. Rev. B 30, 1765 (1984). [105] X. Gonze, J.-P. Michenaud, and J.-P. Vigneron, First-principles study of As, Sb, and Bi electronic properties, Phys. Rev. B 41, 11827 (1990). [106] L. M. Falicov and S. Golin, Electronic Band Structure of Arsenic. I. Pseudopotential Approach, Phys. Rev. 137, A871 (1965). 107 [107] S. Golin, The Electronic Band Structure of Arsenic. II. Self-Consistent Approach, Phys. Rev. 140, A993 (1965). [108] M. Kelly and D. Bullett, Calculation of the Electronic Structure of Crystalline and Amorphous Arsenic, Solid State Commun. 18, 593 (1976). [109] D. Schiferl and C. S. Barrett, Crystal Structure of Arsenic at 4.2 K, 78 K and 299 K, J. Appl. Cryst. 2, 30 (1969). [110] R. Fei and L. Yang, Strain-Engineering the Anisotropic Electrical Conductance of FewLayer Black Phosphorus, Nano Lett. 14, 2884 (2014). [111] C. Lee, X. Wei, J. W. Kysar, and J. Hone, Measurement of the Elastic Properties and Intrinsic Strength of Monolayer Graphene, Science 321, 385 (2008). [112] I. W. Frank, D. M. Tanenbaum, A. M. van der Zande, and P. L. McEuen, Mechanical properties of suspended graphene sheets, J. Vac. Sci. Techn. B 25, 2558 (2007). [113] M. Huang, H. Yan, T. F. Heinz, and J. Hone, Probing Strain-Induced Electronic Structure Change in Graphene by Raman Spectroscopy, Nano Lett. 10, 4074 (2010). [114] R. Boca, P. Hajko, L. Benco, I. Benkovsky, and D. Faktor, Thin layers of grey arsenic: A molecular orbital study, Czechoslovak J. Phys. 43, 813 (1993). [115] S. P. Koenig, R. A. Doganov, H. Schmidt, A. H. Castro Neto, and B. Ozyilmaz, Electric Field Effect in Ultrathin Black Phosphorus, Appl. Phys. Lett. 104, 103106 (2014). [116] F. Xia, H. Wang, and Y. Jia, Rediscovering Black Phosphorus: A Unique Anisotropic 2D Material for Optoelectronics and Electronics, Nature Commun. 5, 4458 (2014). [117] A. S. Rodin, A. Carvalho, and A. H. Castro Neto, Strain-Induced Gap Modification in Black Phosphorus, Phys. Rev. Lett. 112, 176801 (2014). [118] L. Kou, T. Frauenheim, and C. Chen, Phosphorene as a Superior Gas Sensor: Selective Adsorption and Distinct I-V Response, J. Phys. Chem. Lett. 5, 2675 (2014). [119] R. Fei, A. Faghaninia, R. Soklaski, J.-A. Yan, C. Lo, and L. Yang, Enhanced Thermoelectric Efficiency via Orthogonal Electrical and Thermal Conductances in Phosphorene, Nano Lett. 14, 6393 (2014). [120] Z. G. Fthenakis, Z. Zhu, and D. Tomanek, Effect of Structural Defects on the Thermal Conductivity of Graphene: From Point to Line Defects to Haeckelites, Phys. Rev. B 89, 125421 (2014). [121] H. Terrones, M. Terrones, E. Hernandez, N. Grobert, J.-C. Charlier, and P. M. Ajayan, New Metallic Allotropes of Planar and Tubular Carbon, Phys. Rev. Lett. 84, 1716 (2000). [122] Y. Liu, F. Xu, Z. Zhang, E. S. Penev, and B. I. Yakobson, Two-Dimensional MonoElemental Semiconductor with Electronically Inactive Defects: The Case of Phosphorus, Nano Lett. 14, 6782–6786 (2014). 108 [123] V. Tran, R. Soklaski, Y. Liang, and L. Yang, Layer-Controlled Band Gap and Anisotropic Excitons in Few-Layer Black Phosphorus, Phys. Rev. B 89, 235319 (2014). [124] X. Han, H. M. Stewart, S. A. Shevlin, C. R. A. Catlow, and Z. X. Guo, Strain and Orientation Modulated Bandgaps and effective Masses of Phosphorene Nanoribbons, Nano Lett. 14, 4607 (2014). 109