Charged Point Defects in the Flatland: Accurate Formation Energy Calculations in Two-Dimensional Materials

Impurities and defects frequently govern materials properties, with the most prominent example being the doping of bulk semiconductors where a minute amount of foreign atoms can be responsible for the operation of the electronic devices. Several computational schemes based on a supercell approach have been developed to get insights into types and equilibrium concentrations of point defects, which successfully work in bulk materials. Here, we show that many of these schemes cannot directly be applied to two-dimensional (2D) systems, as formation energies of charged point defects are dominated by large spurious electrostatic interactions between defects in inhomogeneous environments. We suggest two approaches that solve this problem and give accurate formation energies of charged defects in 2D systems in the dilute limit. Our methods, which are applicable to all kinds of charged defects in any 2D system, are benchmarked for impurities in technologically important h -BN and MoS 2 2D materials, and they are found to perform equally well for substitutional and adatom impurities.


I. INTRODUCTION
Graphene and other two-dimensional (2D) materials [1] have been found to possess many intriguing properties that can be used in a plethora of applications [2][3][4].Point defects play a central role in determining material characteristics of bulk semiconductors [5] but are equally important for 2D systems, which furthermore consist only of surface and are thus very sensitive to the environment.Moreover, because of quantum confinement of both host and defect wave functions, the role of defects may considerably differ from those in the bulk systems [6,7], affecting electronic conductivity, optical spectra, magnetic response, etc. [8][9][10][11][12][13][14].
Numerous theoretical studies focusing on defects in 2D materials have been carried out [13][14][15][16][17][18], but only very few addressed charged defects [10,[19][20][21].This is astonishing in view of the importance of charged defects in bulk semiconductors.Graphene, being a semimetal, cannot possess charged defects, but semiconducting MoS 2 , insulating h-BN, and many other 2D materials with a gap in the electronic spectrum can have charged vacancies, adatoms, and impurities.In particular, the thermodynamic equilibrium concentration of defects as a function of the Fermi-level position cannot be evaluated without taking into account their charge states.
First-principles calculations are ideal for obtaining this information [22,23].However, it is difficult to carry out such calculations using a supercell approach because of large spurious electrostatic interactions of the charged defect with its periodic images.The previously proposed schemes for overcoming such issues in low-dimensional systems involve either introducing a neutralizing charge to the system [24][25][26] or fixing the potential at the cell boundaries [27][28][29].
In this paper, we demonstrate that a straightforward extrapolation of the formation energies for charged defects in 2D materials to the dilute limit may lead to fundamentally wrong or even diverging results.A large difference in the magnitude of the interactions between the defect images for single-atomic-layer h-BN and MoS 2 sheets is demonstrated and rationalized.We then propose and benchmark several approaches for overcoming these errors and, consequently, for obtaining accurate results, thereby opening a road for calculating the properties of charged defects in these technologically important materials.

II. METHODS
Irrespective of the dimensionality of the system, the prime quantity of interest is the defect formation energy [23,30], defined for defect X at charge state q as Published by the American Physical Society under the terms of the Creative Commons Attribution 3.0 License.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
where E tot is the total energy of the supercell, n i are the differences in the number of atoms, and μ i are their chemical potentials [31].ϵ v is the valence band maximum (VBM) of the pristine host, and ϵ F is the Fermi level with respect to VBM, which is set to zero here.Defect formation energies calculated using Eq. ( 1) are very sensitive to defect concentration in the system, i.e., the employed supercell size.Normally, the goal is to assess the properties of isolated defects, which implies that the localized singleparticle wave functions at the defects do not overlap.Thus, in the dilute limit, we define the formation energy of an isolated defect , where E corr accounts for all spurious electrostatic interactions in the supercell.
Calculations of the total energy were carried out within the framework of the density-functional theory (DFT) with plane waves as a basis set and the projector-augmented wave scheme as implemented in the VASP code [32,33].For BN, the plane-wave cutoff was set to 400 eV and reciprocal space sampled with k-point mesh corresponding to 16 × 16 × 1 or higher in the primitive cell.Corresponding numbers for MoS 2 are 500 eV and 12 × 12 × 1 mesh.Exchange and correlations were treated with the Perdew-Burke-Ernzerhof (PBE) functional [34].The ionic degrees of freedom were fully optimized.We have verified that the defect states containing the added or removed charge remain within the band gap.Although the defect Kohn-Sham levels are sensitive to supercell size in bulk systems [35], this is not the case in 2D materials.Even though the size of the vacuum region affects the potential profile throughout the cell, the potential profile within the material sheet is largely unaffected except for a constant shift.

III. RESULTS AND DISCUSSION
The basic problem in the charged defect calculations in the supercell approach is illustrated in Fig. 1(a), where a charged point defect is introduced into the system along with the neutralizing background so that the total charge is zero and the total energy of the system is thus finite.Different scaling of the supercell leads to different spurious interactions of the defect with its periodic images, as well as with the homogeneously distributed background charge.The calculated formation energies without the energy correction are shown in Fig. 1(b) in the case of a carbon substituting a boron atom in the BN sheet in the þ1 charge state [C B ðþ1Þ; note that E f is negative, in agreement with previous calculations [19,36]].When the lateral size of the supercell is fixed and L z is maximized (the "standard" approach in neutral 2D system simulations), the formation energies vary dramatically and diverge at the limit of infinite vacuum.This is due to a spurious repulsive interaction between the images of the defects within the plane, while the interaction with the neutralizing background charge vanishes.When the vacuum size (i.e., the separation between the periodic images of the sheets in the z direction) is kept fixed and the lateral size is increased, the situation is equally troublesome, as evident from the comparison of the energies for 4 × 4 and 8 × 8 supercells at a fixed L z .In this case, strong interlayer repulsion dominates.Only when the system size is increased uniformly do the energies appear to converge to a finite value, corresponding to E f ∞ .Knowing E f ∞ , we can estimate the error in this type of calculation.For example, using a large 8 × 8 supercell with L z ¼ 15 Å yields 0.4 eV lower formation energy than the correct value.These problems are then directly transferred to the calculation of charge transition levels, as shown in Fig. 1(c).
The method of uniform extrapolation is intuitively simple and generally applicable, although since the asymptotic form is unknown, it can be difficult to obtain an accurate extrapolated value.Furthermore, the method may be somewhat tedious in practice when a large number of defect cases are considered.
As the first approach for correcting the induced errors a posteriori, we calculate the Madelung [37] energy within an effective medium.Using the generalization of the Ewald sum to an arbitrary dielectric tensor [38] and the effective dielectric constants over the whole supercell (from DFT), we obtain the results shown in Fig. 1(b).The energies are overestimated, and the results do not improve when moving to larger supercells.Although the correction will approach zero at the limit of an infinitely large supercell, for the employed set of supercells, it even appears that the results extrapolate to different values.Thus, this approach does not seem to be viable.
A method for correcting electrostatics for charged defects at surfaces and interfaces was introduced in Ref. [39].The method relies on estimating the electrostatic energy of a localized defect charge interacting with its periodic images and the neutralizing charge in a dielectric environment described by a macroscopic dielectric constant profile (cf.Refs.[39,40] for detailed instructions for carrying out the correction).The critical ingredient to be modeled is the profile of the spatially varying macroscopic dielectric constant, in which the atomic level details have been averaged out.However, the applicability of this method for 2D materials is unclear, as these materials lack any bulk region and only exhibit surfaces.In the context of a layer that is one (or a few) atoms thick, averaging over unit cells cannot be carried out in the direction perpendicular to the layer, as there is no longer a clear repeatable unit, and the macroscopic dielectric constant profile becomes rather poorly defined.
To extend the applicability of the method to 2D materials, we propose here a scheme for constructing the dielectric constant profile.It is clear that non-negligible screening can only occur in regions with significant charge density, and no screening should be present within the vacuum.A natural way to achieve this is to choose the dielectric constant profile to follow the charge density distribution of the system.In addition, using elementary electrostatics (see Refs. [41][42][43][44]), we can calculate the effective dielectric constant over the whole supercell from the dielectric constant profile: We then normalize the model profiles in such a way that the effective dielectric constants ε are equal to those obtained from, e.g., DFT-based perturbation-theory calculations [45].To be precise, the normalization proceeds by rescaling the polarizability χðzÞ ¼ ðεðzÞ − 1Þ=ð4πÞ, which guarantees that the vacuum dielectric constant remains 1.
Since the atomic geometries are relaxed in our study, the ion contribution to dielectric response has to be accounted for; i.e., the static dielectric constant is used throughout.This approach is aimed at retaining a model screening response that is consistent with the screening in the DFT calculation, both in the short and in the long range.The normalization was found to lead to corrections largely independent of the exact shape of the profile.A comparison over different model profiles is shown in Ref. [40].
The corrected defect formation energies for C B ðþ1Þ in h-BN are shown in Fig. 1(b).Independent of the supercell size, all energies lie at almost exactly the same energy, proving the good quality of the correction.Importantly, the corrected energies also agree with the extrapolated value.We may now extract the "correct" formation energy for the isolated C B ðþ1Þ defect in BN: −0.77 eV, indicating that the previously reported value of −1.2 eV [19] was somewhat underestimated.Although the spurious stabilization obtained with smaller vacuum sizes is now removed, the formation energy still remains negative in the limit of Fermi energy lying close to the valence band maximum.
For the cells with fixed lateral dimensions, depending on vacuum size, the formation energy can be either below or above the correct formation energy.Consequently, there exists a vacuum size L s z that will yield correct formation energy without any need for energy corrections or extrapolations, owing to cancellation of errors due to stabilizing (defect charge vs neutralizing charge) and destabilizing (defect charge vs its periodic images) interactions.With the knowledge of the correct formation energy, we may determine vacuum size yielding this particular energy, which we call "special vacuum."The values are listed in Table I, and the formation energies with these cells are shown in Fig. 2. The C B ðþ1Þ defect was used for the determination of the special vacuums, and thus for this defect, they work well by construction [Fig.2(a)].In order to investigate if these vacuums work equally well in other  [40] for results on triangular 4 carbon clusters).The formation energies are sometimes slightly underestimated and sometimes overestimated, but they are generally within 0.1 eV of the extrapolated results and thus in marked improvement over the uncorrected values.It is especially interesting to notice that the special vacuums work well even in the case of adatoms for which the dielectric environment is rather different from the case of substitutional defects.We next consider defects in a MoS 2 sheet.This system consists of three atomic layers and thus is considerably thicker.It also clearly exhibits a smaller band gap of about 2 eVand stronger screening.Because of the larger thickness of MoS 2 , the dielectric constant profile is taken to be square (smoothed by an error function at the edges), as shown in Fig. 3(a).We first consider the Mn atom substituting Mo Mn Mo ðþ1Þ, which is of interest as it has been proposed to lead to ferromagnetism [17,18].The results are shown in Fig. 3(c).Qualitatively, they are similar to the case of defects in BN, and the corrected and special vacuum energies lie within some tens of meVs.Contrary to defects in BN, even the uncorrected energies are off by only 0.15 eV because of the larger dielectric constant and lattice constants.Note, however, that the errors are expected to scale as q 2 and will become significant for higher charge states.The special vacuums for MoS 2 are also listed in Table I and were determined using the Mn Mo ðþ1Þ defect, which lies symmetrically in the middle of the MoS 2 layer.
The results for the F substitutional in the S site F S ðþ1Þ are shown in Fig. 3(d) and follow closely those for Mn Mo ðþ1Þ in the main paper.The correction does not seem to succeed to fully remove the linear dependence on L −1 z , but the results are still clearly improved.All calculated formation energies, as well as charge transition levels, are collected in Table II.
The reason for the good performance of the special vacuum scheme lies in the similarity of the electrostatic interactions over the supercell for different kinds of defects, consequently yielding similar errors (except scaled by q 2 ).Thus, it can be sufficient to estimate the error (and thus the correction) in the case of only one defect and use the same value for all defects, scaled by q 2 .For example, from Figs. 3(c) and 3(d), the formation energies for the two defects in the 6 × 6 supercell with 18 Å vacuum are 1.05 and −2.82 eV.Comparing these to the extrapolated values given in Table II, we obtain an error of about ð70 meVÞ • q 2 , independent of the defect.
In the case of two-dimensional materials, the dielectric characteristics of the environment may also significantly affect the properties of the system with defects.Extrapolation is, in principle, valid when the environment is fully included in the calculation, although the system size quickly becomes computationally untractable.Within the energy correction scheme, it is fairly straightforward to introduce any arbitrary dielectric constant for the surroundings.Not only are the formation energies and charge transition levels affected by the dielectric environment, but the band edge states are also [46,47].Extraneous complications may thus arise when attempting to calculate the charge transition level position with respect to band edges.However, the position of the charge transition levels is a rigid quantity within DFT when given with respect to the vacuum level [48][49][50], a quantity readily available from these calculations.We note that charge-preserving excitations, such as optical transitions, are likely to be only very weakly affected by the environment, and these problems are avoided.

IV. CONCLUSIONS
To summarize, we showed that the energetics of charged point defects in 2D materials, when studied using firstprinciples calculations and the supercell approach, may be strongly affected by spurious electrostatic interactions, giving rise to large errors in defect formation energies in the dilute limit.The origin of these errors was rationalized, and their magnitude was illustrated by the examples of charged defects in technologically important 2D h-BN and MoS 2 materials.We further presented approaches for overcoming these errors, which are based on an a posteriori energy correction scheme or on choosing a "special" amount of vacuum between the periodic images of the 2D sheets.Coupled with the other recent advances in defect studies [26,[50][51][52], this work provides a way for accurately calculating charged defect formation energies, equilibrium concentrations, and defect level positions in the rapidly growing family of 2D semiconducting and insulating materials.z .Extrapolations to infinitely large supercells and formation energies using special vacuums (cf.Table I) are also shown.TABLE II.Formation energies and charge transition levels for all considered defects in 2D h-BN and MoS 2 (in eV).For formation energies of charged defects, the Fermi level is set to the valence band maximum.Charge transition levels are referenced to the valence band maximum.Substitutional defect formation energies are taken from the extrapolation, whereas adatom energies are taken from the special vacuum calculation with a 6 × 6 supercell.

FIG. 1 .
FIG. 1.(a) Schematic presentation of a supercell geometry for charged defect calculation in 2D materials by the example of positively charged carbon substitutional impurity C B ðþ1Þ in an h-BN sheet.The different types of interactions between defects are denoted by arrows.(b) Formation energy of C B ðþ1Þ as a function of the inverse layer separation for two different fixed lateral sizes (L x ¼ L y ¼ 4a and 8a, where a is the lattice constant) and for a system where all dimensions of the supercell are scaled uniformly (L x ¼ L y ∼ L z , starting from L x ¼ L y ¼ 4a, L z ¼ 10 Å).Both uncorrected (dashed lines) and corrected (solid lines) energies are shown.Application of Madelung correction within an effective medium is shown with a dash-dotted line.Extrapolations to an infinitely large supercell are also presented (dotted lines).(c) Charge transition levels εð0= þ 1Þ for the same system sizes as in (b).The scale of the y axis is limited by the PBE band gap of BN (4.67 eV).

FIG. 2 .
FIG. 2. (a) A schematic showing the determination of the special vacuum from the crossing point of the known E f ∞ with the curve 4 × 4 supercells for the C B ðþ1Þ defect in h-BN.The results obtained using the special vacuums (from Table I at indicated supercell sizes) are shown with triangles (black).(b) The uncorrected and corrected formation energies for the C N ð−1Þ defect (notation as in Fig. 1) together with the special vacuum results.(c, d) Special vacuum formation energies for C, B, and N adatoms on h-BN at the (c) þ1 and (d) −1 charge states.The respective atomic geometries are visualized on top of each panel.

FIG. 3 .
FIG. 3. (a) Model dielectric constant profile for MoS 2 and comparison to plane-averaged charge density from DFT. Dotted horizontal lines denote the macroscopic dielectric constant of bulk 2H-MoS 2 .(b) Model and DFT potentials, and their difference for the Mn Mo ðþ1Þ defect.(c, d) Uncorrected (dashed lines) and corrected (solid lines) formation energies for Mn Mo ðþ1Þ (c) and F S ðþ1Þ (d) defects as a function of the inverse supercell size in the z direction L −1z .Extrapolations to infinitely large supercells and formation energies using special vacuums (cf.TableI) are also shown.