# Charged Vacancy Defects in Black Phosphorus Monolayer Phosphorene

Biswas Rijal<EMAIL_ADDRESS>Materials Science and Engineering,
University of Florida, Gainesville, Florida 32611, U.S.A  Quantum Theory
Project, University of Florida, Gainesville, Florida 32611, U.S.A  Anne Marie
Z. Tan  School of Mechanical and Aerospace Engineering, Nanyang Technological
University, Singapore 639798, Singapore  Materials Science and Engineering,
University of Florida, Gainesville, Florida 32611, U.S.A  Quantum Theory
Project, University of Florida, Gainesville, Florida 32611, U.S.A  Christoph
Freysoldt  Max Planck Institut für Eisenforschung, Düsseldorf, Germany
Richard G. Hennig<EMAIL_ADDRESS>Materials Science and Engineering,
University of Florida, Gainesville, Florida 32611, U.S.A  Quantum Theory
Project, University of Florida, Gainesville, Florida 32611, U.S.A

###### Abstract

The two-dimensional semiconductor phosphorene has attracted extensive research
interests for potential applications in optoelectronics, spintronics,
catalysis, sensors, and energy conversion. To harness phosphorene’s potential
requires a better understanding of how intrinsic defects control carrier
concentration, character, and mobility. Using density-functional theory and a
charge correction scheme to account for the appropriate boundary conditions,
we conduct a comprehensive study of the effect of structure on the formation
energy, electronic structure and charge transition level of the charged
vacancy point defects in phosphorene. We predict that the neutral vacancy
exhibits a 9-5 ring structure with a formation energy of 1.7 eV and
transitions to a negatively charged state at a Fermi level 1.04 eV above the
valence band maximum. The corresponding optical charge transitions display
sizeable Frank-Condon shifts with a large Stokes shift of 0.3 eV. Phosphorene
vacancies should become negatively charged in $n$-doped phosphorene, which
would passivate the dopants and reduce the charge carrier concentration and
mobility.

##  I Introduction

Black phosphorus is the most stable allotrope of phosphorus and was
synthesized accidentally in 1914 when white phosphorus was subjected to the
high pressure of 1.2 GPa and a temperature of 473 K Bridgman (1914); Ribeiro,
Pimenta, and de Matos (2018); Ling  _et al._ (2015). Like graphite, black
phosphorus is a layered crystal held together by weak van der Waals forces.
One hundred years after its discovery, black phosphorus was mechanically
exfoliated into a few-layered form Li  _et al._ (2014); Liu  _et al._ (2014).
This two-dimensional (2D) form of phosphorus, named phosphorene, is the only
elemental 2D material besides graphene that can be mechanically exfoliated Liu
_et al._ (2014). Phosphorene has been widely studied for various applications
in energy technology Li, Guo, and Li (2020); Batmunkh, Bat-Erdene, and Shapter
(2018), electronic devices Liu  _et al._ (2014), catalysis Ran  _et al._
(2018); Wang  _et al._ (2016), and sensors Liu  _et al._ (2017); Buscema  _et
al._ (2014), due to its high carrier mobility Liu  _et al._ (2014), optical
and electronic anisotropy Tran  _et al._ (2014); Xia, Wang, and Jia (2014),
and tunable bandgap Tran  _et al._ (2014). Fig. 1 shows that single-layer
phosphorene displays an anisotropic puckered honeycomb structure comprised of
two layers of phosphorus atoms bonded by $sp^{3}$ hybridized orbitals Rudenko
and Katsnelson (2014); Boulfelfel  _et al._ (2012).

Figure 1: (Color online) Structure of 2D phosphorene with a single vacancy
point defect. a) The top view shows the 9-5 ring (turquoise lines) formed by
relaxation of the defect structure. b) The side view shows the unrelaxed
structure, and the green circle indicates the atom that is removed to create
the vacancy.

For the use of phosphorene in electronic devices, it is paramount to determine
the formation energy of point defects and understand their effect on the
electronic properties. Obtaining this information from experiments alone is
difficult due to the high concentration of defects and impurities introduced
during synthesis and the difficulty in establishing thermodynamic equilibrium.
Computational methods, such as density-functional theory (DFT), provide a
reliable approach to predict defect properties in semiconductors and
complement experiments Van de Walle and Neugebauer (2004); Freysoldt  _et al._
(2014); Komsa  _et al._ (2014); Freysoldt and Neugebauer (2018). For DFT
calculations of charged defects in 2D materials with plane-wave basis sets,
unphysical electric fields appear in the vacuum region of supercell
calculations that require correction Komsa  _et al._ (2014, 2018); Freysoldt
and Neugebauer (2018).

Several computational and experimental studies investigated neutral point
defects in phosphorene Cai  _et al._ (2016); Babar and Kabir (2016); Ziletti
_et al._ (2015); Kiraly  _et al._ (2017). However, only a few address the
charged defects in phosphorene Zhan  _et al._ (2019); Guo and Robertson
(2015); Gaberle and Shluger (2018) even though in field-effect transistors,
these native point defects can become charged and reduce the charge carrier
mobility Liu  _et al._ (2014). The previous DFT studies disagree on the charge
transition levels (CTLs) of the phosphorene vacancy. Guo and Robertson Guo and
Robertson (2015) reported a [+1/$-1$] CTL 0.24 eV above the valence band
maximum (VBM), while Gaberle and Shluger Gaberle and Shluger (2018) predicted
a [0/$-1$] CTL at 0.55 eV. In the previous studies the effect of defect
configuration on the electronic structure has not been thoroughly
investigated. Therefore, there is a need to conduct a comprehensive study of
the effect of the structure and apply recently developed correction schemes
Freysoldt and Neugebauer (2018) to accurately calculate and validate the
formation energy and charge transition levels of vacancies in phosphorene.

The prevalent approach Cohen  _et al._ (1975) for defect calculations embeds
the defect into periodic supercells of increasing size to reduce the defect-
defect interactions and converge the defect properties to the dilute limit Van
de Walle and Neugebauer (2004); Shim  _et al._ (2005); Leslie and Gillan
(1985). However, additional care must be taken when performing such
calculations for charged defects due to spurious electrostatic interactions
between the defect and its periodic images and the homogeneous compensating
background charge.

In this work, we calculate the vacancy formation energy in phosphorene
considering the $q=-1,0,+1$ charge states. We employ the Freysoldt-Neugebauer
correction method Freysoldt and Neugebauer (2018) to account for the spurious
electrostatic interaction between the defect charge and its images and the
compensating background charge. We show that this method enables calculations
of defect formation energies to within 0.1 eV accuracy in moderately sized
simulation cells of 96 and 140 atoms. We discuss the electronic character of
the defect state in terms of the relevant orbitals and the projected density
of states. We predict a vacancy formation energy of 1.7 eV for the neutral
vacancy and a [0,$-1$] CTL at 1.04 eV above the valence band maximum (VBM).
This deep acceptor level exhibits a strong electron-lattice coupling with a
large Stoke shift of 0.3 eV that could be measured by photoluminescence. The
results show that intrinsic vacancy defects in phosphorene lead to deep
acceptor levels that can compensate n-type dopants, trap charges, and scatter
electrons, thereby reducing the carrier concentration and mobility.

##  II Method

###  II.1 Density Functional Theory

We perform DFT Hohenberg and Kohn (1964); Kohn and Sham (1965) calculations to
characterize the stability, charge transition levels, and density of states of
various structures of the vacancy in phosphorene. We use the plane-wave code
VASP Kresse and Hafner (1993, 1994); Kresse and Furthmüller (1996a, b) with
the projector augmented wave method Blöchl (1994). For the exchange-
correlation functional, we employ the Perdew-Burke-Ernzerhof (PBE) Perdew,
Burke, and Ernzerhof (1996) generalized gradient approximation (GGA) and the
strongly constrained and appropriately normed (SCAN) meta-GGA approximation
Sun, Ruzsinszky, and Perdew (2015). A plane-wave basis cutoff energy of 520 eV
and a $k$-point density of 1000 per reciprocal atom for the Brillouin-zone
integration with a Monkhorst-Pack mesh Monkhorst and Pack (1976) ensure energy
convergence to within 1 meV/atom. We perform spin-polarized calculations and
relax the atomic structures until the energy difference between subsequent
interactions was smaller than 0.001 meV for the simulation cell.

First, we relax both lattice vectors and ionic positions of the pristine
rectangular four-atom unit cell of single-layer phosphorene obtained from
MaterialsWeb.org Ashton  _et al._ (2017) to determine the equilibrium lattice
parameters, shown in Tab. 1. Then, from the relaxed unit cell, we create
supercells of size $4\times 3$, $6\times 4$, and $5\times 7$ with 48, 96, and
140 atoms, respectively. For the defect calculations, we remove one phosphorus
atom per supercell, fix the lattice parameters to the equilibrium lattice
constant of the pristine material, and relax only the ionic positions. This
constraint enables the calculation of the elastic dipole of the defect, i.e.,
the derivative of the defect energy with respect to strain, and hence the
calculation of the formation energy for small changes in the lattice constant
if desired. To determine the energy convergence with respect to vacuum
spacing, we investigate interlayer spacings of 10, 20, and 30 Å.

Table 1: Calculated lattice parameters $a$ and $b$, fundamental bandgap $E_{g}$, and averaged permittivity, $\epsilon_{r}$, of phosphorene, calculated with the PBE and SCAN functionals compared to literature and experimental values. |  $a$ (Å) |  $b$ (Å) |  $E_{g}$ (eV) | $\epsilon_{r}$
---|---|---|---|---
PBE | 3.29 | 4.62 | 0.91 | 15.44
SCAN | 3.28 | 4.60 | 1.28 | 12.50
HSE06 |  |  | 1.6 |
PBE Wang, Kawazoe, and Geng (2015) | 3.30 | 4.61 | 0.91 |
HSE06 Wang, Kawazoe, and Geng (2015) | 3.30 | 4.50 | 1.56 |
$G_{0}W_{0}$ Ferreira and Ribeiro (2017) |  |  | 2.06 |
Exp. Liang  _et al._ (2014); Wang  _et al._ (2015) |  |  | 2, 2.2 |

Table 1 shows the lattice parameters, fundamental bandgap, and permittivity
(static dielectric constant) for phosphorene calculated using the PBE and SCAN
functionals. The lattice parameters match well for both functionals. However,
the calculated band gaps underestimate the experimental fundamental bandgap by
about 1 eV, typical for semi-local exchange-correlation functionals. The meta-
GGA functional SCAN somewhat improves the bandgap, and the hybrid functional
HSE06 Heyd, Scuseria, and Ernzerhof (2003) yields a value closest to
experiment. The experimental optical bandgap of phosphorene is 1.3 eV due to a
strong exciton binding energy of 0.9 eV Wang  _et al._ (2015). To accurately
reproduce the fundamental and optical gap of phosphorene requires
computationally demanding many-body $G_{0}W_{0}$ and Bethe-Salpeter
calculations Ferreira and Ribeiro (2017), beyond the scope of this work.

###  II.2 Defect Formation Energy and Charge Transition Level

The formation energies of point defects in different charge states in
semiconductors determine their charge transition levels (CTLs) and equilibrium
concentrations as a function of the Fermi level Zhang and Northrup (1991); Van
de Walle and Neugebauer (2004). The defect formation energy,
$E_{\mathrm{f}}[X^{q}]$, of a defect $X$ with charge $q$ is given by

$\begin{split}E_{\mathrm{f}}[X^{q}]=&E_{\rm{tot}}[X^{q}]-E_{\rm{tot}}[\textrm{host}]-\\\
&\sum_{i}n_{i}\mu_{i}+q\left(\epsilon_{F}+\epsilon_{v}\right)+E_{\mathrm{corr}},\end{split}$
(1)

where $E_{\textrm{tot}}[X^{q}]$ and $E_{\textrm{tot}}[\textrm{host}]$ are the
total energies of the supercell containing the defect $X$ and the pristine
host structure, respectively, $n_{i}$ denotes the number of atoms of type $i$
added or removed by the creation of the defect, and $\mu_{i}$ denotes their
chemical potentials. For the single-component material phosphorene, the
chemical potential $\mu_{P}$ is the energy per atom of pristine phosphorene.
For the charged defects, the Fermi level, $\epsilon_{F}$, represents the
energy of the electron reservoir in the material, and we reference
$\epsilon_{F}$ to the VBM of the host material, $\epsilon_{v}$. The term
$E_{\mathrm{corr}}$ comprises the correction for the spurious electrostatic
interactions between the charged defect and its periodic images and the
homogeneous compensating background charge Leslie and Gillan (1985); Makov and
Payne (1995); Rozzi  _et al._ (2006); Freysoldt, Neugebauer, and Van de Walle
(2009).

In 2D materials, the compensating charge background generates a quadratic
potential across the vacuum region, which causes an unphysical linear
divergence of the monolayer energy with increasing vacuum spacing Komsa  _et
al._ (2014, 2018); Freysoldt and Neugebauer (2018). We utilize the Freysoldt-
Neugebauer correction, which is applied in postprocessing and neither requires
additional defect calculations nor relies on fitted parameters. The correction
energy is obtained from a Gaussian approximation for the local defect charge
density and aligns the long-range potential calculated for the defect
potential with the model potential Freysoldt, Neugebauer, and Van de Walle
(2009, 2011); Freysoldt and Neugebauer (2018). The method requires as input
the electrostatic potential of the host and defect simulation cell and the
permittivity. We calculate the permittivity using DFT, assuming that the
monolayer is dielectrically isotropic, i.e., the in-plane
($\epsilon_{\parallel}$) and out-of-plane components ($\epsilon_{\bot}$) are
equal. The computed permittivity tensor ($\epsilon_{\textrm{sc}}$) includes
contributions from both the slab and the vacuum spacing in the unit cell, with
the in-plane components acting in parallel and the out-of-plane component
acting in series, resulting in the averaged permittivity for the monolayer of
Freysoldt  _et al._ (2008); Tan, Freysoldt, and Hennig (2020a)

$\epsilon^{\textrm{slab}}=\frac{\epsilon^{\textrm{sc}}_{\parallel}-1}{1-(\epsilon^{\textrm{sc}}_{\bot})^{-1}}.$
(2)

Next, we calculate the thermodynamic charge transition level [$q/q^{\prime}$]
which is defined as the Fermi energy at which the defect formation energy of
the charge state $q$ ($E_{\textrm{f}}(X^{q})$) equals that of the charge state
$q^{\prime}$ ([$E_{\textrm{f}}(X^{q^{\prime}})$]),

$[q/q^{\prime}]=\frac{E_{\mathrm{f}}(X^{q};\epsilon_{F}=0)-E_{\mathrm{f}}(X^{q^{\prime}};\epsilon_{F}=0)}{q^{\prime}-q},$
(3)

where, $E_{\mathrm{f}}(X^{q},\epsilon_{F}=0)$ is the defect formation energy
of $X^{q}$ at the VBM. The thermodynamic or adiabatic transition levels
correspond to sufficiently slow processes in which the system fully relaxes to
its new ground state. Hence the formation energies in Eq. (3) are evaluated
for the relaxed charged defects. For instantaneous transitions occurring on
time scales much shorter than lattice relaxations, such as in optical
excitations, the atomic configuration is frozen. These optical transition
levels are obtained from the formation energies for the different charge
states with the same atomic configuration.Freysoldt  _et al._ (2014) In
general, specific finite-size corrections for optical (vertical) transitions
from charged initial states are required to account for the frozen-in
screening response of the ions Gake  _et al._ (2020); Falletta, Wiktor, and
Pasquarello (2020). In phosphorene, however, the screening is purely
electronic since the Born effective charges vanish for symmetry reasons. The
corrections therefore reduce to the difference of the standard charge
corrections according to the initial and final charge state, respectively.

##  III Results and Discussion

###  III.1 Defect structure

The relaxation of the neutral vacancy defect leads to two different
configurations depending on the initial symmetry. Starting from the pristine
phosphorene supercell with one atom removed and preserving the mirror symmetry
leads to the 5-5-6-6 ring structure Hu and Yang (2015); Babar and Kabir
(2016). However, when breaking the symmetry by perturbing a neighboring atom
to the vacancy, the structure relaxes to the lower energy 9-5 ring
configuration Babar and Kabir (2016) shown in Fig. 1(a). The observed
distortion for the neutral vacancy in phosphorene is similar to the
reconstruction surrounding the vacancy in graphene, which leads to the
formation of a new bond saturating two of the three dangling bonds with a
single dangling bond remaining for geometric reasons Robertson  _et al._
(2013); Banhart, Kotakoski, and Krasheninnikov (2011) and giving rise to a
five- and nine-membered ring structure. However, Guo and Robertson Guo and
Robertson (2015) reported that the negative vacancy in phosphorene does not
display the 9-5 ring structure of the neutral and the positive vacancy. We did
not observe such discrepancy and find that the 9-5 ring structure is the
lowest in energy. We note that if the symmetry is not broken during the
relaxation process, the structure relaxes to the higher energy metastable
5-5-6-6 structure Guo and Robertson (2015). The discrepancy may be due to the
use of different functional as they use screened exchange (sX) hybrid density
functional. Clark and Robertson (2010)

The calculated vacancy defect formation energy for the 5-5-6-6 structure is
1.95 eV, whereas the 9-5 structure has a lower energy of 1.59 eV (PBE). The
calculated defect formation energies for both structures match those reported
in the literature Hu and Yang (2015); Huang  _et al._ (2018); Babar and Kabir
(2016). The new bonds which form the 9-5 rings have bond length of 2.33 Å,
which is larger than the bond length in pristine phosphorene of 2.24 Å. The
relaxed structures of the neutral and +1 charged P vacancy are essentially
identical. However, the ring structure around the negatively charged defect is
further distorted; the 4 bonds on the left side of the ring shrink while the
ones on the right increase in length. The changes in bond lengths are less
than 0.1 Å, hence, the 9-5 ring structure is maintained.

Table 2: Calculated defect formation energy for neutral vacancy defect with the two ring structure. The ring structure depends upon the broken and the unbroken symmetry Structure | PBE | SCAN
---|---|---
9-5 | 1.64 eV | 1.72 eV
5-5-6-6 | 1.95 eV | 2.28 eV

Next, we calculate the formation energies of the charged vacancy defects and
the charge transition levels and characterize the electronic structure of the
lower energy 9-5 ring structure.

###  III.2 Energetics

We compare the corrected defect formation energies for the vacancy in 96 and
140-atom simulation cells with interlayer vacuum spacings of 10, 20, and 30 Å.
For the neutral defect in the 96 and 140-atom cells, the formation energies
changes by less than 10 meV when the vacuum spacing increases from 20 to 30 Å.
For the charged cases, the energy correction removes the linear divergence of
the uncorrected formation energy with increasing vacuum spacing. The corrected
formation energies of the charged vacancy in the 96 and 140-atom cells differ
by 70 and 5 meV, respectively, for vacuum spacings of 20 and 30 Å.

The corrected formation energies for the neutral and negatively charged
vacancy show a similar convergence with respect to supercell size, indicating
that the energy correction removes the in-plane electrostatic interactions
between the charged defects. The formation energy of the neutral vacancy
changed from 1.64 eV for the 96 atom cell to 1.59 eV for the 140 atom cell.
The formation energy for the positively charged vacancy exhibited a stronger
dependence on the in-plane cell size due to the defect state being below the
VBM, resulting in the excess charge being delocalized in the 2D material,
which is not accounted for by the correction scheme. Overall, we estimate that
the defect formation energies are converged to within about 0.1 eV for the
supercells of 96 atoms and vacuum spacings of 20 Å and use this supercells for
the characterization of the vacancy.

Figure 2: (Color online) Defect formation energies of the charged vacancy
defects in phosphorene as a function of Fermi level calculated with the PBE
and SCAN functional. The resulting thermodynamic and optical charge transition
levels are compared with the valence and conduction band edges from PBE, SCAN,
and HSE06 using the vacuum level as a common reference. The vacancy only shows
a [0/$-1$] charge transition inside the bandgap, close to the conduction band
minimum. The differences in energy between the thermodynamic and the optical
transition levels corresponding to absorption $0\rightarrow-1$ and emission
$-1\rightarrow 0$ indicates large Franck-Condon shifts for the deep acceptor
level of the vacancy.

Fig. 2 shows the formation energies of the charged vacancy obtained with the
PBE and SCAN functional and compares the resulting thermodynamic and optical
charge transition levels to the band edges obtained with the PBE, SCAN, and
HSE06 functional using vacuum level as the common reference. The vacancy
exhibits a [0/$-1$] charge transition inside the bandgap at 0.65 and 1.04 eV
above the VBM for the PBE and SCAN functional, respectively. The [0/$-1$]
charge transition occurs close to the conduction band minimum (CBM),
indicating that the vacancy in phosphorene acts as a deep acceptor. The [+1/0]
charge transition falls well inside the valence band region, indicating that
positively charged vacancies do not occur in monolayer phosphorene. For the
6-6-5-5 ring structure, which is not the ground state, the charge transition
levels are outside the band gap region.

Figure 3: (Color online) DOS for (a) the neutral and (b) the negatively
charged vacancy projected onto the $s$ and $p$ orbitals of the atoms nearest
to the vacancy defect calculated with the SCAN functional. The grey shaded
regions indicate the DOS of pristine phosphorene. The energies are plotted
with respect to the vacuum level. The dashed line indicate the Fermi level.
(c)–(f) Isosurfaces of the charge density for the defect levels labeled A, B,
and C (top and bottom view of the slab). (g,h) The top and bottom view of the
simulated STM image for peak B in the density of states. The STM images are
nearly identical for the defect level in the bandgap of the neutral and
negatively charged vacancy.

Optical transitions occur on time scales too fast for structural relaxations.
Therefore, the energy difference between the neutral and negatively charged
vacancy without relaxations provides an estimate of the optical CTLs for
absorption [$0\rightarrow-1$] of $-4.1$ eV and emission [$-1\rightarrow 0$] of
$-4.4$ eV relative to the vacuum level. The large observed shift in the
optical CTL relative to the thermodynamic CTL in Fig. 2 indicates sizable
Franck-Condon shifts and a combined Stokes shift of 0.3 eV for the phosphorene
vacancy that could be measured by photoluminescence.

The intrinsic phosphorene vacancy defect could be neutral or negatively
charged depending upon the Fermi level position. This could lead to two
possible excitation and emission processes. When the defect is negatively
charged, such as in $n$-type doped phosphorene, an electron could be excited
from the filled defect state to the CBM absorbing a photon of 0.38/0.57 eV
energy (SCAN/HSE06), and the emission would occur when an electron jumps from
the CBM to the defect state with the emitted photon energy of 0.09/0.28 eV
(SCAN/HSE06). Similarly, when the defect is charge neutral, such as in undoped
or $p$-doped phosphorene, an electron could be excited from the VBM to the
empty defect state, creating a hole in the VBM and absorbing 0.84/1.03 eV
energy (SCAN/HSE06), and the emission would occur when an electron moves from
the defect state to fill the VBM hole emitting a photon of 1.13/1.32 eV energy
(SCAN/HSE06). The SCAN functional predicts that this entire process takes
place in the infrared and the microwave regime depending upon the nature of
the defect charge. The more accurate band edges from the HSE06 functional
indicate that the emission and absorption occur in the infrared region.

With the Stoke’s shift in the infrared region, phosphorene could have
applications in biosensing and solar concentrators. Doped graphene quantum
dots based fluorescence biosensor have shown remarkable biosensing ability
that have potential applications in biomedical applications Kalkal  _et al._
(2020). Also, tandem solar concentrators with two different absorption onsets,
one in the near infrared region and the other in the visible region have
demonstrated improved power conversion efficiency for collecting sunlight Wu,
Li, and Klimov (2018).

###  III.3 Density of States

To characterize the electronic structure of the vacancy, we calculate the
electronic density of state (DOS) and charge density. Fig. 3(a,b) shows the
DOS projected onto the $s$ and $p$ orbitals of the P atoms neighboring the
vacancy for the negatively charged and neutral defect, illustrating that the
defect state predominantly exhibits $p_{x}$ and $p_{y}$ orbital character. For
the neutral vacancy, the singly occupied defect level is spin-polarized. The
occupied state (A) occurs inside the valence band while the empty state (B)
remains in the bandgap. For the negative vacancy, the defect level (C) becomes
doubly occupied and is inside the bandgap. Fig. 3(c-e) illustrates the charge
density of the defect orbitals corresponding to peaks A, B, and C. The defect
charge density is strongly localized on the $p$ states of the P atoms
surrounding the vacancy and extends in the $x$ and $y$ directions, validating
the DOS results. The DOS is shown with respect to the vacuum level.
Calculating the vacuum level for the neutral case is straightforward; however,
care must be taken when calculating the vacuum energy for the charged cases.
First we obtain the vacuum level in the pristine supercell. Then we subtract
the contribution from the electrostatic artifact obtained from FNV correction
scheme to remove the unphysical quadratic potential and account for the
alignment term due to the background charges.

Comparing defect orbitals obtained from DFT with scanning tunneling microscopy
(STM) $dI/dV$ images can help verify a defect’s type and nature Kozhakhmetov
_et al._ (2020). Fig. 3(f,g) shows the simulated STM images from the top and
the bottom of the monolayer obtained from the projected charge density between
0.25-0.4 eV above the VBM that corresponds to the peak B in the density of
states. The top view shows a characteristic asymmetry in the simulated STM
image for the defect state inside the bandgap of the neutral vacancy with the
5-9 ring structure that can help experimental validation. The STM images are
nearly identical for the negatively charged vacancy. The image shows that the
vacancy defect state is dominated by a $p$ orbital localized on the
neighboring P atom, similar to its charge density.

Importantly, the energy of the peaks for the defect states in the DOS differ
from the CTL’s since the DOS and band structure illustrate the one-electron
Kohn-Sham energies while the true CTL’s correspond to total energy differences
of different charge states that include the many-electron contributions to the
energy from exchange and correlation and the charge correction. Therefore, the
DOS cannot accurately predict the position of CTL’s.

The choice of exchange-correlation functionals affects the position of the
band edges and CTL’s. In this study, PBE and SCAN predict similar defect
formation energies leading to the qualitatively same conclusions. Similarly,
we do not expect the usually more accurate HSE06 functional to change the
conclusion of our results. In our previous study of the vacancy in 2D MoS2,
the HSE06 functional predicts a larger bandgap than SCAN but retains the
[0/$-1$] CTL at the same position relative to the vacuum level.Tan, Freysoldt,
and Hennig (2020b) If assuming a similar behavior in phosphorene, the [0/$-1$]
CTL would occur about 1.2 eV above the VBM.

We compare our structure and charge transition level with two previous
studies. Guo and Roberson Guo and Robertson (2015) reported a [+1/$-1$] CTL
within the bandgap. However, the applied charge correction Lany and Zunger
(2008); Persson  _et al._ (2005) did not account for the varying dielectric
profile and large vacuum spacing specific to 2D materials. Gaberle and Shluger
Gaberle and Shluger (2018) reported a [0/$-1$] CTL at 0.55 eV above the VBM,
which closely agrees with our result of 0.7/1.04 eV for PBE/SCAN. For charge
correction, they extrapolate the defect formation energy to the infinite
supercell. Guidon, Hutter, and VandeVondele (2009) The difference in the CTL
energy may be due to the choice PBE0-TC-LRC hybrid functional and the charge
correction method.

##  IV Conclusions

We computed the relaxed strucure and defect formation energies for charged
single vacancies in monolayer phosphorene, utilizing the charge correction
scheme developed by Freysoldt and Neugebauer to restore the appropriate
electrostatic boundary conditions. We find that the symmetry broken 9-5 ring
configuration is the ground state structure for the neutral and charged
vacancies. The neutral vacancy has a formation energy of 1.7 eV and undergoes
a charge transition to the negative state at a Fermi level of 1.04 eV above
the valence band maximum within the bandgap. Hence, phosphorene’s single
vacancy acts as a deep acceptor that could passivate dopants in $n$-type
phosphorene and reduce the carrier concentration and mobility. We predict that
the defect level in the bandgap exhibits $p$-orbital character that can be
revealed by scanning tunneling microscopy. The optical charge transition
levels display sizeable Frank-Condon shifts and a large Stokes shift of more
than 0.3 eV, providing possible opportunities for biosensing and photovoltaics
applications.

###### Acknowledgements.

This work was supported by the National Science Foundation under Grant Nos.
DMR-1748464 and OAC-1740251 and the 2DCC-MIP under Grant No. DMR-1539916.
Computational resources were provided by the University of Florida Research
Computing Center. Part of the research was performed while the authors visited
the Institute for Pure and Applied Mathematics (IPAM), which is supported by
the National Science Foundation (Grant No. DMS-1440415).

## References

  * Bridgman (1914) P. Bridgman, J. Am. Chem. Soc. 36, 1344 (1914).
  * Ribeiro, Pimenta, and de Matos (2018) H. B. Ribeiro, M. A. Pimenta, and C. J. de Matos, J. Raman Spectrosc. 49, 76 (2018).
  * Ling  _et al._ (2015) X. Ling, H. Wang, S. Huang, F. Xia, and M. S. Dresselhaus, Proc. Natl. Acad. Sci. 112, 4523 (2015).
  * Li  _et al._ (2014) L. Li, Y. Yu, G. J. Ye, Q. Ge, X. Ou, H. Wu, D. Feng, X. H. Chen, and Y. Zhang, Nat. Nanotechnol. 9, 372 (2014).
  * Liu  _et al._ (2014) H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tománek, and P. D. Ye, ACS Nano 8, 4033 (2014).
  * Li, Guo, and Li (2020) J. Li, C. Guo, and C. Li, ChemSusChem 13, 1047 (2020).
  * Batmunkh, Bat-Erdene, and Shapter (2018) M. Batmunkh, M. Bat-Erdene, and J. G. Shapter, Adv. Funct. Mater. 8, 1701832 (2018).
  * Ran  _et al._ (2018) J. Ran, W. Guo, H. Wang, B. Zhu, J. Yu, and S.-Z. Qiao, Adv. Mater. 30, 1800128 (2018).
  * Wang  _et al._ (2016) X. Wang, G. Sun, N. Li, and P. Chen, Chem. Soc. Rev. 45, 2239 (2016).
  * Liu  _et al._ (2017) X. Liu, T. Ma, N. Pinna, and J. Zhang, Adv. Funct. Mater. 27, 1702168 (2017).
  * Buscema  _et al._ (2014) M. Buscema, D. J. Groenendijk, S. I. Blanter, G. A. Steele, H. S. Van Der Zant, and A. Castellanos-Gomez, Nano Lett. 14, 3347 (2014).
  * Tran  _et al._ (2014) V. Tran, R. Soklaski, Y. Liang, and L. Yang, Phys. Rev. B 89, 235319 (2014).
  * Xia, Wang, and Jia (2014) F. Xia, H. Wang, and Y. Jia, Nat. Commun. 5, 1 (2014).
  * Rudenko and Katsnelson (2014) A. N. Rudenko and M. I. Katsnelson, Phys. Rev. B 89, 201408 (2014).
  * Boulfelfel  _et al._ (2012) S. E. Boulfelfel, G. Seifert, Y. Grin, and S. Leoni, Phys. Rev. B 85, 014110 (2012).
  * Van de Walle and Neugebauer (2004) C. G. Van de Walle and J. Neugebauer, J. Appl. Phys. 95, 3851 (2004).
  * Freysoldt  _et al._ (2014) C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti, and C. G. Van de Walle, Rev. Mod. Phys. 86, 253 (2014).
  * Komsa  _et al._ (2014) H.-P. Komsa, N. Berseneva, A. V. Krasheninnikov, and R. M. Nieminen, Phys. Rev. X 4, 031044 (2014).
  * Freysoldt and Neugebauer (2018) C. Freysoldt and J. Neugebauer, Phys. Rev. B 97, 205425 (2018).
  * Komsa  _et al._ (2018) H.-P. Komsa, N. Berseneva, A. V. Krasheninnikov, and R. M. Nieminen, Phys. Rev. X. 8, 039902 (2018).
  * Cai  _et al._ (2016) Y. Cai, Q. Ke, G. Zhang, B. I. Yakobson, and Y.-W. Zhang, J. Am. Chem. Soc. 138, 10199 (2016).
  * Babar and Kabir (2016) R. Babar and M. Kabir, J. Phys. Chem. C 120, 14991 (2016).
  * Ziletti  _et al._ (2015) A. Ziletti, A. Carvalho, D. K. Campbell, D. F. Coker, and A. C. Neto, Phys. Rev. Lett. 114, 046801 (2015).
  * Kiraly  _et al._ (2017) B. Kiraly, N. Hauptmann, A. N. Rudenko, M. I. Katsnelson, and A. A. Khajetoorians, Nano Lett. 17, 3607 (2017).
  * Zhan  _et al._ (2019) F. Zhan, W. Xu, R. Zou, J. Yang, J. Fan, X. Wu, and R. Wang, J. Phys. Chem. C 123, 27080 (2019).
  * Guo and Robertson (2015) Y. Guo and J. Robertson, Sci. Rep. 5, 1 (2015).
  * Gaberle and Shluger (2018) J. Gaberle and A. L. Shluger, Nanoscale 10, 19536 (2018).
  * Cohen  _et al._ (1975) M. L. Cohen, M. Schlüter, J. R. Chelikowsky, and S. G. Louie, Phys. Rev. B 12, 5575 (1975).
  * Shim  _et al._ (2005) J. Shim, E.-K. Lee, Y. Lee, and R. M. Nieminen, Phys. Rev. B 71, 035206 (2005).
  * Leslie and Gillan (1985) M. Leslie and N. Gillan, J. Phys. C: Solid State Phys. 18, 973 (1985).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Hafner (1994) G. Kresse and J. Hafner, Phys. Rev. B 49, 14251 (1994).
  * Kresse and Furthmüller (1996a) G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996a).
  * Kresse and Furthmüller (1996b) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996b).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Perdew, Burke, and Ernzerhof (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Sun, Ruzsinszky, and Perdew (2015) J. Sun, A. Ruzsinszky, and J. P. Perdew, Phys. Rev. Lett. 115, 036402 (2015).
  * Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  * Ashton  _et al._ (2017) M. Ashton, J. Paul, S. B. Sinnott, and R. G. Hennig, Phys. Rev. Lett. 118, 106101 (2017).
  * Wang, Kawazoe, and Geng (2015) V. Wang, Y. Kawazoe, and W. Geng, Phys. Rev. B 91, 045433 (2015).
  * Ferreira and Ribeiro (2017) F. Ferreira and R. M. Ribeiro, Phys. Rev. B 96, 115431 (2017).
  * Liang  _et al._ (2014) L. Liang, J. Wang, W. Lin, B. G. Sumpter, V. Meunier, and M. Pan, Nano Lett. 14, 6400 (2014).
  * Wang  _et al._ (2015) X. Wang, A. M. Jones, K. L. Seyler, V. Tran, Y. Jia, H. Zhao, H. Wang, L. Yang, X. Xu, and F. Xia, Nat. Nanotechnol. 10, 517 (2015).
  * Heyd, Scuseria, and Ernzerhof (2003) J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 118, 8207 (2003).
  * Zhang and Northrup (1991) S. Zhang and J. E. Northrup, Phys. Rev. Lett. 67, 2339 (1991).
  * Makov and Payne (1995) G. Makov and M. Payne, Phys. Rev. B 51, 4014 (1995).
  * Rozzi  _et al._ (2006) C. A. Rozzi, D. Varsano, A. Marini, E. K. Gross, and A. Rubio, Phys. Rev. B 73, 205119 (2006).
  * Freysoldt, Neugebauer, and Van de Walle (2009) C. Freysoldt, J. Neugebauer, and C. G. Van de Walle, Phys. Rev. Lett. 102, 016402 (2009).
  * Freysoldt, Neugebauer, and Van de Walle (2011) C. Freysoldt, J. Neugebauer, and C. G. Van de Walle, Phys. Status Solidi (b) 248, 1067 (2011).
  * Freysoldt  _et al._ (2008) C. Freysoldt, P. Eggert, P. Rinke, A. Schindlmayr, and M. Scheffler, Phys. Rev. B 77, 235428 (2008).
  * Tan, Freysoldt, and Hennig (2020a) A. M. Z. Tan, C. Freysoldt, and R. G. Hennig, Phys. Rev. Mater. 4, 064004 (2020a).
  * Gake  _et al._ (2020) T. Gake, Y. Kumagai, C. Freysoldt, and F. Oba, Phys. Rev. B 101, 020102(R) (2020).
  * Falletta, Wiktor, and Pasquarello (2020) S. Falletta, J. Wiktor, and A. Pasquarello, Phys. Rev. B 102, 041115 (2020).
  * Hu and Yang (2015) W. Hu and J. Yang, J. Phys. Chem. C 119, 20474 (2015).
  * Robertson  _et al._ (2013) A. W. Robertson, B. Montanari, K. He, C. S. Allen, Y. A. Wu, N. M. Harrison, A. I. Kirkland, and J. H. Warner, ACS Nano 7, 4495 (2013).
  * Banhart, Kotakoski, and Krasheninnikov (2011) F. Banhart, J. Kotakoski, and A. V. Krasheninnikov, ACS Nano 5, 26 (2011).
  * Clark and Robertson (2010) S. J. Clark and J. Robertson, Phys. Rev. B 82, 085208 (2010).
  * Huang  _et al._ (2018) L.-Y. Huang, X. Zhang, M. Zhang, and G. Lu, Phys. Rev. Materials 2, 054003 (2018).
  * Kalkal  _et al._ (2020) A. Kalkal, R. Pradhan, S. Kadian, G. Manik, and G. Packirisamy, ACS Appl. Bio Mater. 3, 4922 (2020).
  * Wu, Li, and Klimov (2018) K. Wu, H. Li, and V. I. Klimov, Nat. Photonics 12, 105 (2018).
  * Kozhakhmetov  _et al._ (2020) A. Kozhakhmetov, B. Schuler, A. M. Z. Tan, K. A. Cochrane, J. R. Nasr, H. El-Sherif, A. Bansal, A. Vera, V. Bojan, J. M. Redwing, N. Bassim, S. Das, R. G. Hennig, A. Weber-Bargioni, and J. A. Robinson, Adv. Mater. 32, 2005159 (2020).
  * Tan, Freysoldt, and Hennig (2020b) A. M. Z. Tan, C. Freysoldt, and R. G. Hennig, Phys. Rev. Mater. 4, 114002 (2020b).
  * Lany and Zunger (2008) S. Lany and A. Zunger, Phys. Rev. B 78, 235104 (2008).
  * Persson  _et al._ (2005) C. Persson, Y.-J. Zhao, S. Lany, and A. Zunger, Phys. Rev. B 72, 035211 (2005).
  * Guidon, Hutter, and VandeVondele (2009) M. Guidon, J. Hutter, and J. VandeVondele, J. Chem. Theory Comput. 5, 3010 (2009).