Native Point Defects in Mono-- and Bi--layer Phosphorene
NNative Point Defects in Mono– and Bi–layer Phosphorene
Sudipta Kundu, Mit H. Naik,
1, 2, 3 and Manish Jain ∗ Centre for Condensed Matter Theory, Department of Physics,Indian Institute of Science, Bangalore 560012, India Department of Physics, University of Californiaat Berkeley, Berkeley, California 94720, USA Material Sciences Division, Lawrence BerkeleyNational Laboratory, Berkeley, California 94720, USA
Abstract
We study the stability and electronic properties of intrinsic point defects, vacancy and self–interstitial, in mono– and bi–layer phosphorene. We calculate the formation energies, quasiparticledefect states and charge transition levels (CTLs) of these defects using ab initio density functionaltheory (DFT) and GW approximation to the electron self–energy. Using the DFT + GW two pathsformalism for studying interstitial in monolayer phosphorene, we show that with the inclusionof electrostatic corrections CTLs can be calculated reliably. Our calculations show that all thenative point defects have low formation energies 0.9–1.6 eV in neutral state. Furthermore, wefind that vacancy in phosphorene behaves as an acceptor–like defect which can explain the p–type conductivity in phosphorene. On the other hand, interstitial can show both acceptor– anddonor–like behaviour. ∗ [email protected] a r X i v : . [ c ond - m a t . m t r l - s c i ] J a n . INTRODUCTION The discovery of graphene [1, 2] has evoked tremendous research in two–dimensional(2D) layered materials. Apart from introducing new physics, 2D materials show promisingapplications in opto–photonic, nano–photonic, sensor devices etc [3, 4]. However, the semi–metallic nature of graphene limits its wide application in electronic devices despite its veryhigh carrier mobility [5]. Transition metal dichalcogenides [6–8] are another importantmember of this 2D materials family. Unlike graphene, they have a finite band gap, buttheir low mobility [9] constrains their applications. A few years ago, another 2D material,phosphorene, was exfoliated [10–13] from bulk black phosphorus. Phosphorene, the singleor few layer form of black phosphorus, has been drawing much attention since then. In alayer of phosphorene, a phosphorus atom is bonded to three other phosphorus atoms via sp hybridization. This gives rise to a puckered honeycomb structure [14]. The uniquenessof phosphorene is its structural anisotropy which manifests in anisotropic optical, thermaland transport properties [15–17]. Phosphorene has a direct band gap at the center of theBrillouin zone which ranges from ∼ ∼ / Vs and a high on/off current ratio up to 10 [10]. The tunabilityof band gap and mobility depending on the number of layers make the material suitable forelectronic, opto–electronic, photo–electric, and FET device applications [20, 21].Presence of defect affects the material properties in different ways. Depending on thecharge state of the defect, it can induce free carriers or trap charge and act as scattererin the system. Electronic, optical and magnetic properties can also be altered due to thepresence of defects. Hence, an extensive study of defects and their possible charge statesis required for defect engineering and proper understanding of material properties. Thepuckered structure of phosphorene facilitates the formation of various defects [22]. Defectsin phosphorene have low formation energy [23] in comparison to defects in graphene [24].There is experimental evidence of vacancies [25] in phosphorene. While defects can escalatethe degradation process of devices [26, 27] such as vacancy can increase the oxidization ofphosphorene, they can also enhance device performance by introducing desired properties;for example, emission of photons at new frequencies at room temperature [28]. Vacancy and2elf–interstitial are two low energy native point defects in phosphorene. In general, DFTworks well for calculating ground state properties like structures of defects. However, theestimation of the correct charge transition level (CTL) requires the calculation of the energyassociated with the change in the electron number at the defect site. This is an excited stateproperty and is not expected to be estimated correctly with Kohn–Sham DFT. While Heyd–Scuseria–Ernzerhof (HSE) functional [29] is an improvement over the generalized gradientapproximation [30] and works well for bulk materials, it fails to capture the anisotropicscreening in two dimensional materials properly [31, 32]. The combined formalism of manybody perturbation theory within GW approximation [33, 34] and DFT has emerged as areliable method for calculation of CTL [35, 36]. Furthermore, previous DFT calculationshave either been performed with a small supercell size [37] or have not taken into account thefull anisotropic (not only out-of-plane but also in-plane) dielectric constant of phosphorenewhile correcting for spurious electrostatic interaction [38–40].We study structural and electronic properties of vacancy and self–interstitial in mono–and bi–layer phosphorene. We calculate the formation energies and CTLs of these defectsusing the DFT + GW formalism. Since we employ periodic boundary conditions in ourcalculations, the charged defect calculations suffer from spurious electrostatic interactionswith their periodic neighbours. This problem is more prominent in a 2D material due toreduced screening along the out-of-plane direction. We correct for this spurious interactionsystematically by properly modeling the anisotropic dielectric medium. Further, we validateour calculation of CTLs by evaluating CTL of interstitial in monolayer phosphorene followingtwo different paths. With the electrostatic correction to defect levels, the two paths agreewithin ±
100 meV. We find that defects in phosphorene have low formation energies, in therange 0.9–1.6 eV, which is consistent with previous calculations [22]. Our study also showsthat vacancy can induce p-doping supporting the experimental finding of intrinsic p-typebehaviour of phosphorene [10, 11, 21]. In contrast, interstitial in mono– and bi–layer canact as both p-type acceptor and n-type donor. Furthermore, we find that the unoccupieddefect levels of neutral defect calculated with DFT and GW line up with respect to thevacuum level while the occupied defect levels shift downward in GW compared to DFTcalculated levels. This manuscript is organized as follows. Section 2 describes our methodsof calculation which include the details about DFT, GW and the corrections for chargeddefect calculations. Calculation and results of defects in mono– and bi–layer phosphorene3re presented in section 3 and 4 respectively.
II. COMPUTATIONAL DETAILS AND METHODOLOGY:
All the DFT calculations in this study are performed using the Quantum Espresso package[41]. We use a norm–conserving pseudopotential [42] and the exchange–correlation potentialis approximated by the generalized gradient approximation proposed by Perdew, Burke andErnzerhof (GGA–PBE) [30]. The van der Waals interaction between layers is described by“Grimme–D2” method [43]. The wave functions are expanded in a plane wave basis, withplane waves up to energy of 60 Rydberg included in the basis. For unit cell calculationswe adopt a 14 × × × × × × − eV and 0 . / ˚A respectively.The GW [33, 34] calculations are performed using the BerkeleyGW package [45]. Thenumber of unoccupied bands are 100 and 125 for unit cell in monolayer and bilayer re-spectively. This choice ensures the convergence of the band gap to be better than 0.1 eV.The dielectric matrix is expanded in plane wave with energy up to 12 Ry and extended tofinite frequencies using generalized plasmon pole (GPP) [46] model. The Coulomb inter-action along the out of-plane-direction is truncated to compute the dielectric matrix andself–energy [47]. The static remainder technique is used to accelerate the convergence of thecalculation with the number of empty bands [48]. The Brillouin zone sampling used for GWcalculations of unit cell is 21 × ×
1. We find the quasiparticle band gap for monolayer tobe 2.07 eV which is in good agreement with previous theoretical [49–51] and experimental(scanning tunneling spectroscopy [52] and photoluminescence excitation spectroscopy [53])4tudies. The band gap for bilayer is 1.29 eV which is also consistent with previous calcu-lations [50, 51]. For the defect calculations with 7 × × × fq [ (cid:126) R q ](E F ) = E q [ (cid:126) R q ] − E pristine + N P µ P + q(E F + (cid:15) VBM ) + E corr (q) (1)where E pristine is the total energy per supercell of pristine mono– or bi–layer phosphorene,E q [ (cid:126) R q ] is the total energy of a supercell containing defect in charge state q at its relaxed coor-dinates (cid:126) R q , N P is the number of phosphorus atoms added to (or removed from) the supercellto create the defect, µ P is the chemical potential of a phosphorus atom which is calculatedfrom bulk black phosphorus, (cid:15) VBM is the energy of the valence band maximum (VBM) ofpristine cell, E F is the Fermi level with respect to (cid:15) VBM and E corr (q) is the electrostaticcorrection term.As discussed previously, simulations of charged defects suffer from an erroneous electro-static interaction between the defect cell and its images arising due to periodic boundaryconditions. As a consequence, formation energies and eigenvalues of defect levels in chargeddefects are estimated incorrectly. The correction term E corr (q) in Eqn. (1) removes thespurious electrostatic interactions and includes the potential alignment in formation energy.E corr (q) is calculated using the Freysoldt– Neugebauer– Van de Walle (FNV) scheme [56] asimplemented in CoFFEE code [57]. To calculate the correction term, mono– and bi–layerphosphorene are modeled as a dielectric slab of width 5.22 ˚A and 10.44 ˚A respectively. Thedielectric constants of the slab [58] are calculated using density functional perturbation the-ory of the Quantum Espresso package [41]. The in-plane dielectric constants for monolayerphosphorene are 12.5 and 18.0 along the zigzag direction and armchair direction, respec-tively. The out-of-plane dielectric constant is 1.9. These quantities for bilayer are 13.5,27.4 and 1.9 respectively. The difference in the in-plane dielectric constants is due to thestructural anisotropy of phosphorene. The charge is modeled as Gaussian distribution whoseparameters are taken from DFT calculations. Using this model, the electrostatic energy iscalculated for different cell sizes with a uniform scaling parameter, α , and extrapolated toinfinite limit using 5th order polynomial to obtain the value for an isolated defect [58]. The5ifference between the isolated value and the value for a particular cell size gives the requiredcorrection in formation energy for the corresponding cell. To verify our model, we calculatethe electrostatic energies of negatively charged interstitial defect in bilayer phosphorene witha model charge for several supercell sizes. Our starting cell is the 7 × α = 7)for which the in-plane lattice parameters are almost equal. We uniformly scale the cell sizealong all the three directions and calculate the electrostatic energy (blue curve in Fig. 1(a)).Further, we perform a different set of scaling calculations starting from a supercell suchthat the vacuum scaling is different. The scaling factor in these cells along the out-of-planedirection is 1.5 times of those in the in-plane directions (red curve in Fig. 1(a)). We observethat the electrostatic correction does not change monotonically with the inverse of super-cell length which necessitates the several calculations of large supercell to reach the infinitelength limit. As a result, a simple 1/L extrapolation from a small supercell size can leadto incorrect isolated value [59, 60]. Fig. 1(b) shows the formation energies of negativelycharged interstitial in bilayer phosphorene without and with correction for supercell sizes7 × , × ×
10. The corrections have been obtained following the blue curveof Fig. 1(a). It can be clearly seen that the corrected formation energies are same for allthe three supercells. For comparison we have further included the formation energies ofneutral defect in bilayer phosphorene for the same supercell sizes as mentioned above. Theformation energies of neutral defect also have no supercell size dependence.The defect CTL ε (q / q (cid:48) ) is defined as the Fermi–level position for which the formationenergies of charge states q and q (cid:48) are equal: ε (q / q (cid:48) ) = E fq [ (cid:126) R q ]( E F = 0) − E fq (cid:48) [ (cid:126) R q (cid:48) ]( E F = 0)q (cid:48) − q (2)The CTL can be computed from formation energies within DFT. However, due to band gapunderestimation within DFT and the fact that CTL involves change in electron number,significant error arises in the computed CTL. Within the combined DFT and GW approach,the CTL [61] is written as: ε (q / q (cid:48) ) = (E fq [ (cid:126) R q ] − E fq (cid:48) [ (cid:126) R q ]) + (E fq (cid:48) [ (cid:126) R q ] − E fq (cid:48) [ (cid:126) R q (cid:48) ])q (cid:48) − q (3)by adding and subtracting term E fq (cid:48) [ (cid:126) R q ] in the numerator. If the charge states q and q’differ by ±
1, the term (E fq [ (cid:126) R q ] − E fq (cid:48) [ (cid:126) R q ]) can be identified as the quasiparticle energy (E QP )which is calculated using GW. This term accounts for an electron removal or addition to the6 .00 0.04 0.08 0.12 M o d e l E n e r g y ( e V ) (a) F o r m a t i o n E n e r g y ( e V ) UncorrectedCorrectedNeutral (b)
FIG. 1. (a) shows the variation of electrostatic energy of bilayer phosphorene with differentsupercell sizes. α is the scaling factor. Blue and red curves represent the variation with scaling α × α × α and α × α × . α respectively. In (b) the formation energies of negatively chargedinterstitial in bilayer phosphorene without and with correction for three different cell sizes areshown by blue circles and red diamonds and extrapolated to infinite limit. The formation energyof neutral interstitial are shown by green squares. system. The other term (E fq (cid:48) [ (cid:126) R q ] − E fq (cid:48) [ (cid:126) R q (cid:48) ]) captures the relaxation energy (E relax ) of thestructure due to the change in electron number and is calculated using DFT. Depending onthe addition or removal of the electron (Eqn. 3), CTL can be calculated following differentpaths (Fig. 2). The parabolas in Fig. 2 represent the formation energies as a functionof the generalized coordinates of the atoms in the cell. Along one path, we start with thedefect in charge state q at its equilibrium coordinates , (cid:126) R q . The E QP is represented by thegreen vertical arrow, and E relax by the green curved arrow. The CTL can also be computedstarting with the defect in charge state q’, as shown by the red arrows. Since CTL is athermodynamic quantity, the CTL computed using the two paths should be the same. Notethat if the GW calculation is performed on charged defects the defect eigenvalues have to becorrected. We calculate the CTLs for interstitial in monolayer phosphorene following twopaths [35] starting from neutral defect (q=0) and charged defect supercells (q’=+1 or –1).With the corrected quasiparticle values, the CTLs from different paths agree to within ±
100 meV. The details of the calculations are given in subsequent section.7 / = E f [ R ] - E f [ R ] ⃗ R q ⃗ R q ′ E fq [ ⃗ R q ] E fq ′ [ ⃗ R q ′ ] E fq E fq ′ ε ( q / q ′ ) E f ⃗ R FIG. 2. The figure shows calculation of CTLs following different paths. The straight green and redarrows represent the quasiparticle energies calculated at the equilibrium structures of the defectwith charge state q and q’ respectively. The curved arrows account for relaxation energies.
III. DEFECTS IN MONOLAYER PHOSPHORENEA. Interstitial
The most stable structure of interstitial in monolayer phosphorene is when the interstitialatom forms two symmetric bonds with two phosphorus atoms in one of the sub-layer andresides above the layer (Fig. 3(a) and 3(b)). Formation of two bonds leaves the interstitialatom with a lone electron which can induce defect states. Spin–polarized calculations revealtwo defect states in the band gap. The defect states are localized along one diagonal of thesupercell centering the interstitial atom as is shown in the plot of wave function (Fig. 3(a)and 3(b)). One of the defect states is occupied in neutral state making interstitial to bestable in all three charge states (+1,0,–1). Fig 3(c) shows the density of states (DOS) forall the charge states. Due to the dispersion along the diagonal direction, the defect bandsare broad. Upon accepting an electron, interstitial atom breaks the symmetry and thebonds become asymmetric. There is no significant relaxation when the defect is positivelycharged. In the negatively charged interstitial, the defect state appears flatter. This is dueto the fact that the defect wave function is more localized. While the defect states arespin split in the neutral state, they become degenerate in both positively and negativelycharged states as the electrons get paired in these states. The interstitial defect has lowformation energy of 1.63 eV in neutral state. To obtain the formation energy in charged8 a) (b) D O S ( s t a t e s / e V / s u p e r c e ll ) I M0 I M+1
Energy (eV) I M1 (c) Fermi Energy (eV) F o r m a t i o n E n e r g y ( e V ) I M+1 I M0 I M1 G W , V B M G W , C B M D F T , V B M D F T , C B M (d) FIG. 3. (a) Top view and (b) side view of interstitial defect in monolayer phosphorene withdefect wave function. The isosurface is plotted at the charge density value of 3.4 × − e/˚A . Theinsterstitial atom is shown in yellow color and the neighbor atoms are in red. (c) Spin polarizedDOS of the interstitial in 3 charged states. VBM is set to zero and is shown in black dottedline. Fermi level is marked in maroon dashed line. Red and blue are for two spin states. (d)CTLs calculated within DFT and DFT+GW formalism. The dashed and solid lines represent theformation energies calculated with DFT and DFT+GW respectively and red, black and blue areused for +1, 0 and –1 charge states respectively. Formation energy of neutral defect remains same.The VBM and CBM are also marked in the figure. We have set the vacuum to zero in the plot. states, we evaluate the correction term ( E corr ) using CoFFEE as detailed in section 2. Thecorrection using model charge calculation is 0.295 eV. Taking the potential alignment intoaccout, E corr is calculated to be 0.30 eV and 0.26 eV for positively and negatively chargedstates respectively. The variation of formation energies of interstitial in charged states, asthe Fermi energy is tuned, is shown in Fig. 3(d). Within DFT (red and blue dashed line),the defect changes its state from positive to neutral ( ε (+1 / ε (0 / − ABLE I. CTLs of interstitial along two different paths are reported. GW calculation on neutraldefect is denoted by path 1 and path 2 represents calculations starting from charged defects. Allthe energies are reported in eV. E QP is quasiparticle energy and E Relax is the energy associatedwith relaxing the defect. The correction for the eigenvalue is given by (cid:15) corr .Path 1 Path 2E QP E relax (cid:15) corr CTL E QP E relax (cid:15) corr CTL Avg. CTL Difference+1/0 0.71 0.25 0.00 0.96 1.78 -0.19 -0.59 1.00 0.98 0.040/-1 1.75 -0.27 0.00 1.48 0.46 0.31 0.59 1.36 1.42 0.12 energies of positively and negatively charged interstitial respectively within DFT + GW.The formation energy of neutral defect is same in DFT and DFT + GW. The DFT + GWCTLs ε (+1 /
0) and ε (0 / −
1) are at 0.96 eV and 1.48 eV above VBM respectively. In thisDFT + GW calculation, the E QP ’s are calculated from the GW calculation of the ionizationpotential and electron affinity of neutral interstitial. CTLs can be also calculated startingfrom charged defects (Eqn. 3) [35] as described before. The CTLs from two paths and theirconstituent energies are reported in Table 1. To obtain ε (+1 / ε (+1 /
0) is calculated to be 1.00 eV. Similarly, starting from negatively charged interstitialand calculating ionization energy we obtain ε (0 / −
1) to be 1.36 eV. The difference betweenthe CTLs obtained from the neutral and charged cell calculations lie within the error of GWcalculation (0.1 eV). As a result, the CTLs due to interstitial in monolayer phosphorene aredeep in the gap. The presence of both ε (+1 /
0) and ε (0 / −
1) CTLs implies that interstitialin monolayer phosphorene can show both donor– and acceptor–type behaviors.
B. Vacancy:
Upon removing one phosphorus atom from phosphorene, there are two possible mono–vacancy structures: MV-(55 |
66) (Fig. 4(a) and 4(b)) and MV-(5 |
9) (Fig. 4(c) and 4(d)) [22].MV-(55 |
66) has a symmetric structure. In this structure, the atom closest to the vacancyis connected with 4 P atoms instead of the usual coordination of 3. In contrast, the closestatom to the vacancy in MV-(5 |
9) moves and the system rearranges such that it bonds with 310 a) (b) (c) (d)
Reaction Coordinate F o r m a t i o n E n e r g y ( e V ) (e) D e n s i t y o f S t a t e s ( p e r e V ) V M0 V M+1
Energy (eV) V M1 (f) Fermi Energy (eV) F o r m a t i o n E n e r g y ( e V ) V M0 V M1 G W , V B M G W , C B M D F T , V B M D F T , C B M (g) FIG. 4. (a) and (b) show the top and side view of MV-(55 |
66) in 7 × |
9) and the wave function localized at defect from top and sideview respectively. The isosurface of wave functions is plotted at 3.4 × − e/˚A . The neighbourto the vacancy is shown in yellow and the other neighbour atoms are shown in red. (e) showsthe variation of formation energy while the structure changes from MV-(55 |
66) to MV-(5 | atoms. The formation energy of MV-(5 |
9) is lower than that of MV-(55 |
66) by 330 meV perdefect. In order to estimate the barrier between the two structures, we perform a climbingimage nudged elastic band (CI-NEB) calculation as implemented in PASTA package [62].Fig. 4(e) shows how the formation energy changes as MV-(55 |
66) transforms to the morestable MV-(5 |
9) structure. The energy barrier between the two structures is found to beonly 5 meV. Due to this low barrier, we expect the defect to always be in the MV-(5 | | |
9) vacancy in 3 charge states: neutral, positive and negative.11e perform spin–polarized calculations. The neutral vacancy gives rise to an unoccupiedflat band in the band gap. The DOS calculation (Fig. 4(f)) shows a state 0.23 eV above theVBM. We plot the corresponding wave function in Fig. 4(c) and 4(d). It is clear that thewave function of this site is localized around the vacancy along a diagonal of the supercell.We can understand the origin of this state as follows: a P atom neighbouring the vacancyforms two bonds and has a dangling bond which gives rise to the flat band in the gap. Thefilled defect state lies deep within the valence band and we found that the vacancy is notstable in positive charge state. The defect can accommodate an electron and change itsstate to negatively charged state. Upon accepting an electron, the system rearranges itselfsuch that both the defect bands are in the band gap and they become degenerate (Fig.4(f)). It is to be noted that MV-(55 |
66) also induces defect states but those are within thevalence band. These states hybridize with the valence band states making the defect alwaysnegatively charged. However, this configuration is always higher in energy than negativelycharged MV-(5 | |
9) in neutral and negatively charged states.In Fig. 4(g)) the black line represents the formation energy of neutral vacancy. The vacancyhas a formation energy of 1.42 eV in the neutral state. Due to this low formation energy,vacancy defects are expected to be abundant in phosphorene [25]. The variation of formationenergy in negatively charged state with the Fermi energy is shown in Fig. 4(g) within DFT(blue dashed line) and DFT+GW (blue solid line). The CTL ε (0 / −
1) is 0.62 eV abovethe VBM within DFT and 1.06 eV above the VBM within DFT + GW. This implies thatvacancy in monolayer phosphorene behaves as deep acceptor.
IV. DEFECTS IN BILAYER PHOSPHORENEA. Interstitial
Interstitial in bilayer has a similar structure to that of interstitial in monolayer phos-phorene (Fig. 5(a) and 5(b)). The interstitial atom prefers to bond with two atoms in theouter sub-layer of bilayer phosphorene facing vacuum. The configuration with the intersti-tial atom between the layers is higher in energy by 580 meV per defect. Like the monolayer,spin–polarized calculation on interstitial in bilayer shows that it also induces two defect12 a) (b) D O S ( s t a t e s / e V / s u p e r c e ll ) I B0 I B+1
Energy (eV) I B1 (c) Fermi Energy (eV) F o r m a t i o n E n e r g y ( e V ) I B+1 I B0 I B1 G W , V B M G W , C B M D F T , V B M D F T , C B M (d) FIG. 5. (a) and (b) are the top and side view of interstitial defect with the localized wave functionin bilayer phosphorene. (c) Spin–polarized DOS of interstitial in bilayer phosphorene in positive,neutral and negative charged states. VBM are set to zero and the Fermi level are shown in maroondashed line. Red and blue lines represent two spin states. (d) CTLs calculated within both DFTand DFT+GW. Solid lines denote the formation energies within DFT+GW and the dashed linesare used for DFT. Red, black and blue are for +1, 0 and –1 charge states. Vacuum is set to zero. states in the gap. One of the defect states is filled rendering interstitial to be possible in+1, 0, –1 charge states. DOS calculations on the three charge states reveal that while thedefect states are non-degenerate in neutral state they become degenerate by accepting ordonating an electron in charged states (Fig. 5(c)). In the neutral charge state, the formationenergy is 1.45 eV which is lower than the corresponding defect in the monolayer. Fig. 5(d)shows the formation energies of interstitial in its three charge states +1, 0 and –1 withinboth DFT and DFT+GW. We find the +1 charged interstitial is always higher in energywhen the Fermi energy is in band gap in DFT. In contrast, the interstitial is more stable inpositively charged state within DFT+GW when Fermi level is below 0.43 eV with respectto VBM which marks the CTL ε (+1 / ε (0 / −
1) is at 0.86 eV (0.33 eV) aboveVBM within DFT+GW (DFT). This suggests that the interstitial in bilayer can act as bothdonor– and acceptor–type defect like the interstitial defect in monolayer.13 a)(b) E n e r g y ( e V ) M X Y 1.000.750.500.250.000.250.500.751.00 M X Y (c) (d)
Fermi Energy (eV) F o r m a t i o n E n e r g y ( e V ) V B0 V B1 G W , V B M G W , C B M H S E , V B M H S E , C B M (e) FIG. 6. (a) and (b) show the top and side view of MV-(5 |
9) in bilayer phosphorene. The defectwave function is plotted at isosurface value of 3.4 × − e/˚A . (c) and (d) show the band diagramsof the defect structure relaxed with PBE and HSE respectively. Red and blue lines denote twospin states. (e) is the formation energy plot of vacancy in neutral and negative state. B. Vacancy
In bilayer phosphorene, the vacancy can reside in two inequivalent positions. In oneconfiguration the vacancy faces the second layer while in the other it faces vacuum. Wefound the vacancy is most stable in MV-(5 |
9) structure facing vacuum (Fig. 6(a) and 6(b)).Spin–polarized calculation within GGA-PBE on this structure does not show any distinctdefect state in the band gap (Fig. 6(c)). This calculation suggests hybridization betweendefect state and states inside valence band edge. However, after a diagonal GW correction,a state emerges in the gap. This poses a problem to the diagonal GW approximation (Eqn.4), in which the self–energy matrix element is computed using the DFT wave functions as:E
QPi = E
DFTi + (cid:104) ψ DFTi | Σ GW ( { E DFT , ψ
DFT } ; E QP ) | ψ DFTi (cid:105) − (cid:104) ψ DFTi | V XC | ψ DFTi (cid:105) (4)One way to overcome this problem is to perform self–consistent GW calculation. Anotherway is to choose a better mean–field. To address the issue we start with alternative mean–field. We use hybrid functional approximation (HSE) [29] for exchange–correlation whichis an improvement over GGA. We relax the vacancy structure with HSE. For all furthercalculations, the total energies are calculated with HSE. The obtained band structure gives14 distinct state in the band gap, however the state is very close to valence bands (Fig. 6(d)).Vacancy has a formation energy of 0.97 eV in neutral state. The calculated CTL ε (0 / − G W approximation is notgoing to be adequate. We construct the GW Hamiltonian in the DFT wave function basisfollowing Eqn. 5.H ij = E DFTi δ ij + (cid:104) ψ DFTi | Σ GW ( { E DFT , ψ
DFT } ; E QP ) | ψ DFTj (cid:105) − (cid:104) ψ DFTi | V XC | ψ DFTj (cid:105) (5)Due to computational cost, we restrict the construction of Hamiltonian matrix with wavefunctions with energy within ±
500 meV of the defect state as we are interested in thedefect state and the states close to it. The self–energy matrix is evaluated within static–screening limit (static–COHSEX) [34] and we diagonalize the constructed Hamiltonian [63].The eigenfunctions of the Hamiltonian are the new wave functions with which we performed G W calculation. It should be noted that in this way, the wave functions are iterated whilethe dielectric screening is still constructed from mean–field wave functions. The quasiparticleenergies are obtained by evaluating the self–energy using the standard plasmon pole model.For the vacancy in bilayer, G W calculation shows that there is mixing between the defectstate and the valence bands close to it. After the first iteration with the updated wavefunctions, the quasiparticle energy of the defect state is 370 meV away from the VBM edge.Instead of calculating the CTL using the quasiparticle energy at the equilibrium structure,we can calculate that at any structure and adjust the relaxation energy accordingly (Fig. 7)[64]. We can denote an intermediate structure by (cid:126) R I and rewrite the Eqn. (2) in followingway: ε (q / q (cid:48) ) = (E fq [ (cid:126) R q ] − E fq [ (cid:126) R I ]) + (E fq [ (cid:126) R I ] − E fq (cid:48) [ (cid:126) R I ]) + (E fq (cid:48) [ (cid:126) R I ] − E fq (cid:48) [ (cid:126) R q (cid:48) ])q (cid:48) − q (6)= E relax , q + E QP + E relax , q (cid:48) (7)Fig. 7 shows two energy curves for two charge states of the defect. (cid:126) R I is an intermediatestructure at which the quasiparticle energy is calculated (brown straight arrow in Fig. 7).15 IG. 7. The figure shows calculation of CTL starting from an intermediate structure. The brownarrow represents the quasiparticle energy calculated at (cid:126) R I and the curved blue arrows are therelaxation energy. To obtain the CTL, the relaxation energies from the (cid:126) R I to (cid:126) R q and (cid:126) R q (cid:48) (blue curved arrows inFig. 7) are taken into account. For the vacancy calculation in bilayer phosphorene, we adoptthis scheme to calculate the CTL within DFT + GW. The intermediate structure is chosenas discussed above. We perform the GW calculation with the intermediate structure to getthe quasiparticle energy (E QP ) and calculate the relaxation energy from that structure toequilibrium neutral and negatively charged structure as the vacancy is stable in neutral andnegative state. Using this formalism (Eqn. (6)), the CTL ε (0 / −
1) is found to be at 0.24eV with respect to VBM within DFT + GW (Fig. 6(e)) which again implies that vacancyin bilayer phosphorene behaves like an acceptor.
V. CONCLUSION
We have extensively studied the formation and electronic properties of vacancy and self–interstitial defects in mono– and bi–layer phosphorene. We have taken into account thespurious electrostatic correction while studying charged defects. The defects have formationenergies between 0.9 eV – 1.6 eV in neutral state. Depending on the charge state, these de-fects can further lower their formation energies. It has been also observed that the formationenergy of defects in bilayer phosphorene is smaller than that in monolayer suggesting thatwith the increase of layer number formation energy decreases. We calculate the CTLs of thedefects within DFT and DFT+GW formalism which suggest that while vacancy behaves as16cceptor type defect, interstitial can act as both acceptor and donor type defect.
ACKNOWLEDGMENTS
We thank the Supercomputer Education and Research Centre (SERC) at IISc for pro-viding the computational resources. [1] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V.Grigorieva, and A. A. Firsov, Science , 666 (2004).[2] A. K. Geim and K. S. Novoselov, Nat. Mater , 183 (2007).[3] F. V. Kusmartsev, W. M. Wu, M. P. Pierpoint, and K. C. Yung, “Application of graphenewithin optoelectronic devices and transistors,” in Applied Spectroscopy and the Science ofNanomaterials , edited by P. Misra (Springer Singapore, 2015) pp. 191–221.[4] F. Schedin, A. K. Geim, S. V. Morozov, E. W. Hill, P. Blake, M. I. Katsnelson, and K. S.Novoselov, Nature Materials , 652 (2007).[5] F. Schwierz, Nature Nanotechnology , 487 (2010).[6] K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Phys. Rev. Lett. , 136805 (2010).[7] R. Branimir and A. Kis, Nat. Mater. , 815 (2013).[8] Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman, and M. S. Strano, Nature Nanotech-nology , 699 (2012).[9] B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, and A. Kis, Nature Nanotechnology , 147 (2011).[10] L. Li, Y. Yu, G. J. Ye, Q. Ge, X. Ou, H. Wu, D. Feng, X. H. Chen, and Y. Zhang, Nat. Nano. , 374 (2014).[11] H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tomnek, and P. D. Ye, ACS Nano , 4033(2014).[12] A. Castellanos-Gomez, L. Vicarelli, E. Prada, J. O. Island, K. L. Narasimha-Acharya, S. I.Blanter, D. J. Groenendijk, M. Buscema, G. A. Steele, J. V. Alvarez, H. W. Zandbergen, J. J.Palacios, and H. S. J. van der Zant, 2D Materials , 025001 (2014).
13] M. Buscema, D. J. Groenendijk, G. A. Steele, H. S. van der Zant, and A. Castellanos-Gomez,Nature Communications , 4651 (2014).[14] A. Brown and S. Rundqvist, Acta Crystallographica , 684 (1965).[15] F. Xia, H. Wang, and Y. Jia, Nature Communications , 4458 (2014).[16] Z. Luo, J. Maassen, Y. Deng, Y. Du, R. P. Garrelts, M. S. Lundstrom, P. D. Ye, and X. Xu,Nature Communications , 8572 (2015).[17] R. Schuster, J. Trinckauf, C. Habenicht, M. Knupfer, and B. B¨uchner, Phys. Rev. Lett. ,026404 (2015).[18] L. Liang, J. Wang, W. Lin, B. G. Sumpter, V. Meunier, and M. Pan, Nano Letters , 6400(2014).[19] M. H. Naik and M. Jain, Phys. Rev. B , 165125 (2017).[20] J. Qiao, X. Kong, Z.-X. Hu, F. Yang, and W. Ji, Nature Communications , 4475 (2014).[21] S. Das, W. Zhang, M. Demarteau, A. Hoffmann, M. Dubey, and A. Roelofs, Nano Letters , 5733 (2014).[22] W. Hu and J. Yang, The Journal of Physical Chemistry C , 20474 (2015).[23] F. Banhart, J. Kotakoski, and A. V. Krasheninnikov, ACS Nano , 26 (2011).[24] Y. Liu, F. Xu, Z. Zhang, E. S. Penev, and B. I. Yakobson, Nano Letters , 6782 (2014).[25] B. Kiraly, N. Hauptmann, A. N. Rudenko, M. I. Katsnelson, and A. A. Khajetoorians, NanoLetters , 3607 (2017).[26] A. A. Kistanov, Y. Cai, K. Zhou, S. V. Dmitriev, and Y.-W. Zhang, 2D Materials (2017).[27] K. L. Utt, P. Rivero, M. Mehboudi, E. O. Harriss, M. F. Borunda, A. A. Pacheco SanJuan,and S. Barraza-Lopez, ACS Central Science , 320 (2015).[28] J. Pei, X. Gai, J. Yang, X. Wang, Z. Yu, D.-Y. Choi, B. Luther-Davies, and Y. Lu, NatureCommunications , 10450 (2016).[29] J. Heyd, G. E. Scuseria, and M. Ernzerhof, The Journal of Chemical Physics , 8207(2003).[30] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. , 3865 (1996).[31] M. Jain, J. R. Chelikowsky, and S. G. Louie, Phys. Rev. Lett. , 216806 (2011).[32] P. De´ak, E. Khorasani, M. Lorke, M. Farzalipour-Tabriz, B. Aradi, and T. Frauenheim, Phys.Rev. B , 235304 (2019).[33] L. Hedin, Journal of Physics: Condensed Matter , R489 (1999).
34] M. S. Hybertsen and S. G. Louie, Phys. Rev. B , 5390 (1986).[35] M. Jain, J. R. Chelikowsky, and S. G. Louie, Phys. Rev. Lett. , 216803 (2011).[36] M. H. Naik and M. Jain, Phys. Rev. Materials , 084002 (2018).[37] V. Wang, Y. Kawazoe, and W. T. Geng, Phys. Rev. B , 045433 (2015).[38] Y. Guo and J. Robertson, Scientific Reports , 14165 (2015).[39] D. Wang, D. Han, X.-B. Li, N.-K. Chen, D. West, V. Meunier, S. Zhang, and H.-B. Sun,Phys. Rev. B , 155424 (2017).[40] J.-H. Yang and B. I. Yakobson, ArXiv e-prints (2017).[41] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L.Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. Gironcoli, S. Fabris, G. Fratesi, R. Gebauer,U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri,R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero,A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys. Condens. Matter , 395502 (2009).[42] D. R. Hamann, M. Schl¨uter, and C. Chiang, Phys. Rev. Lett. , 1494 (1979).[43] S. Grimme, J. Comput. Chem. , 1787 (2006).[44] H. J. Monkhorst and J. D. pack, Phys. Rev. B. , 5188 (1976).[45] J. Deslippe, G. Samsonidze, D. A. Strubbe, M. Jain, M. L. Cohen, and S. G. Louie, ComputerPhysics Communications , 1269 (2012).[46] M. S. Hybertsen and S. G. Louie, Phys. Rev. B , 5390 (1986).[47] S. Ismail-Beigi, Phys. Rev. B , 233103 (2006).[48] J. Deslippe, G. Samsonidze, M. Jain, M. L. Cohen, and S. G. Louie, Phys. Rev. B , 165124(2013).[49] F. Ferreira and R. M. Ribeiro, Phys. Rev. B , 115431 (2017).[50] V. Tran, R. Soklaski, Y. Liang, and L. Yang, Phys. Rev. B , 235319 (2014).[51] V. Tran, R. Fei, and L. Yang, 2D Materials , 044014 (2015).[52] L. Liang, J. Wang, W. Lin, B. G. Sumpter, V. Meunier, and M. Pan, Nano Letters , 6400(2014).[53] X. Wang, A. M. Jones, K. L. Seyler, V. Tran, Y. Jia, H. Zhao, H. Wang, L. Yang, X. Xu, andF. Xia, Nature Nanotechnology , 517 (2015).[54] A. Stathopoulos and J. R. McCombs, ACM Trans. Math. Softw. , 21:1 (2010).
55] L. Wu, E. Romero, and A. Stathopoulos, SIAM Journal on Scientific Computing , S248(2017).[56] C. Freysoldt, J. Neugebauer, and C. G. Van de Walle, Phys. Rev. Lett. , 016402 (2009).[57] M. H. Naik and M. Jain, Computer Physics Communications , 114 (2018).[58] J.-Y. Noh, H. Kim, and Y.-S. Kim, Phys. Rev. B , 205417 (2014).[59] H.-P. Komsa, N. Berseneva, A. V. Krasheninnikov, and R. M. Nieminen, Phys. Rev. X ,031044 (2014).[60] H.-P. Komsa, N. Berseneva, A. V. Krasheninnikov, and R. M. Nieminen, Phys. Rev. X ,039902 (2018).[61] M. Jain, J. R. Chelikowsky, and S. G. Louie, Phys. Rev. Lett , 216803 (2011).[62] S. Kundu, S. Bhattacharjee, S.-C. Lee, and M. Jain, Computer Physics Communications ,261 (2018).[63] M. Jain, J. Deslippe, G. Samsonidze, M. L. Cohen, J. R. Chelikowsky, and S. G. Louie, Phys.Rev. B , 115148 (2014).[64] T. Biswas and M. Jain, Phys. Rev. B , 144102 (2019)., 144102 (2019).