Native defects in monolayer GaS and GaSe: electrical properties and thermodynamic stability

Structural, electronic and thermodynamic properties of native defects in GaS and GaSe monolayers are investigated by means of accurate ab-initio calculations. Based on their charge transition levels we assess the influence of the studied defects on the electrical properties of the monolayers. Specifically, we show that native defects do not behave as shallow dopants and their presence cannot account for the experimentally observed intrinsic doping. In addition, we predict that native defects are efficient compensation and recombination centers. Besides pointing out their detrimental nature, we also calculate the corresponding finite temperature formation energies and provide a window of growth conditions able to reduce the concentration of all relevant native defects.


INTRODUCTION
Ultrathin two-dimensional (2D) materials exhibit unusual characteristics that show great promise for application in next generation electronic and optoelectronic devices [1]. The quintessential 2D material is graphene [2,3], a monolayer of carbon atoms arranged in a honeycomb lattice, first exfoliated in 2004 [4] and intensively studied due to its impressive electrical, optical, magnetic and mechanical properties, mostly associated with its electronic bands exhibiting linear dispersion at the Fermi level [5][6][7][8][9]. Unfortunately, in its pristine semimetallic form, graphene does not posses a band gap and its applicability in semiconductor technology is currently limited. Thus, there is great interest in finding and understanding novel semiconducting 2D materials.
Beyond graphene, a prototypical example of 2D materials are single-layer transition metal dichalcogenides (TMDs) MX 2 (e.g., M = Mo, W; X = S, Se) which are direct semiconductors with high absorbance and exhibiting novel optical properties, including valley-selective circular dichroism and coupling of spin and valley degrees of freedom [10][11][12][13][14]. Similar to TMDs but with even more exotic properties, group III-VI metal monochalcogenides MX (M = Ga, In; X = S, Se, Te) have recently attarcted interest due to their novel electrical and optical features [15,16]. Specifically, what truly sets group III-VI metal monochalcogenides apart from the greater family of 2D materials, is the direct to indirect band gap transition observed in both experiments and simulations when their structure reach a critical number of layers [17][18][19][20][21][22][23][24]. Below this critical number of layers, the conduction band minimum (CBM) remains at the brillouin zone center (Γ-point) while the valence band maximum (VBM) moves away from the Γ towards the K -point. This process creates a nonparabolic dispersion for holes near the top of the valence band, usually referred to as "Mexican hat" dispersion [16]. Due to the large degeneracy of states associated with it, such electronic dispersion results in a large den-sity of states (DOS) near the band edge, signaling a van Hove singularity [17,18,25] which can lead to magnetic or superconducting instabilities [26][27][28]. Furthermore, the singularity in the DOS and the large number of conducting modes at the band edge can result in enhanced thermoelectric properties of these materials [29].
In particular, the ultrathin gallium chalcogenides GaS and GaSe have been found to have exceptionally large photoresponse, which makes them excellent candidates to build efficient flexible photodetectors [30][31][32][33][34], and strong nonlinear optical properties [35,36]. Moreover, these materials have been reported to exhibit potential as building blocks of photovoltaic devices [37], as channel materials in field effect transistors [38], as photocalysts for water splitting [39][40][41][42] and can transition into a topological insulator state by means of strain and oxygen functionalization [43,44]. From the point of view of their synthesis, single-layers of GaS and GaSe have been obtained by means of micromechanical cleavage technique [30,32,45,46], pulsed laser deposition (PLD) [34], molecular beam epitaxy (MBE) [47][48][49] and vapor phase transport (VPT) [31,33,[50][51][52]. The last two techniques has so far proved to be the most succesful approaches to obtain large-area high-quality monolayers of GaS and GaSe [16]. As any other material and regardless of the growth method used to synthesize them, ultrathin gallium monochalcogenides are expected to contain point defects which affect their electrical, magnetic and optical properties [53,54]. However, although the impact of defects in 2D materials is particularly strong due to the reduced dimensionaly, there is currently no clear understanding of their properties in the case of monolayer GaS and GaSe. This despite the fact that some of the most interesting properties and future applications of these materials are closely linked to the presence of defects, e.g., their topological insulator state [44] and their potential as material platform for single photon emitters [16,55,56]. In addition, on a more basic level, the origin and nature of intrinsic doping in monolayer GaS and GaSe is still under debate and native defects have been proposed as a arXiv:2110.10238v1 [cond-mat.mtrl-sci] 19 Oct 2021 possible explanation [22,49,57].
Although some valuable attemps to shed light on these issues have come from theoretical studies [58][59][60][61][62], their results are hampered by at least one of the following issues: reduced-size supercells (which causes artificially delocalized defect states), incorrect description of the dielectric properties of the monolayers, severe under-or overestimation of the band gap, incorrect or lack of correction for the electrostatic interaction between charged defects and their periodic images, and incomplete assessment of the possible charge states of the defects. Furthermore, the influence of temperature on defect formation is crucial and closely related to the growth process by which the material is prepared. Therefore, a complete and accurate inspection of the native defects, as the one available for MoS 2 [63], is still lacking in the case of GaS and GaSe. This constitutes the relevance and pertinence of the systematic theoretical study presented in this paper.
By means of accurate density funtional theory (DFT) calculations we obtain physically sound thermodynamic charge transition levels of relevant native defects in monolayer GaS and GaSe. With such information at hand we are able to reveal their electrical properties and to determine whether they are detrimental. In addition, we calculate their formation energies at experimentally relevant conditions, allowing us to predict which point defects are more likely to exist depending on the growth method. More importanly, our results enable us to propose changes in growth conditions such that harmful defects are avoided.
The paper starts with a detailed description of the methods used in Sec. II, where we provide a summary of the computational parameters and a precise account of the approach used to calculate defects formation energies. This includes the proper choices of chemical potentials, the correction used to mitigate the finite-size effects and the technique employed to obtain the dielectric constants of the monolayers. Afterwards, our results are reported in Sec. III where we present the basic structural, thermodynamical and electrical properties of all considered defects along with their formation energies at finite temperature. Furthermore, in Sec. III, we discuss the implications of our findings. Finally, in Sec. IV, we summarize our key results.

A. Computational setup and basic parameters
Density functional theory (DFT) calculations were performed using the VASP [64] simulation package with projector augmented-wave potentials (PAWs) for the effective potential associated to the nucleus and the core electrons. A converged plane-wave energy cutoff of 500 eV was applied and the convergence criterion for forces during ionic relaxation was set to 0.01 eV/Å. Both the Perdew-Burke-Enzerhof (PBE) functional with van der Waals interactions included via the empirical correction proposed by Grimme (PBE-D) and the Heyd, Scuseria and Enzerhof (HSE06) hybrid functional were used as exchange-correlation functionals. On the one hand, as reported in Refs. [65][66][67], the PBE-D approach delivers accurate structural parameters and binding energies for layered materials and it allows us to correctly describe the properties of defects lying outside of the monolayer [63]. On the other hand, although HSE06 calculations do not account for the van der Waals interactions and are more computationally demanding than PBE-D ones, they are useful to accurately assess the electronic properties of defects as they predict band gaps and defect level positions closer to experiments and beyond-DFT methods [68][69][70][71]. Throughout this paper, it is clearly pointed out when and why one or the other functional were used. For all results presented in this work spin polarization was taken into account. Since spin-orbit interaction has been found to be mostly relevant for the description of low-lying valence bands [22,29], they are not expected to be crucial for the thermodynamics of defects and they were not included in our calculations. Basic parameters of the monolayers were calculated using a converged 16 × 16 × 1 k -point grid for the unit cell with 25Å of vacuum and are reported in Table I. As can be seen there, both PBE-D and HSE06 deliver similar in-plane lattice parameters, a, and slightly underestimate the experimental room-temperature values by ∼0.5%. Regarding their electronic properties, optical band gaps measured via catholuminescence experiments are 3.4 eV and 3.3 eV for single-layer GaS and GaSe, respectively [72]. In the case of fundamental gaps, scanning tunneling spectroscopy measurements were conducted and arrived to a value of 3.5 eV in the case of monolayer GaSe [57], which implies an exciton binding energy of 0.2 eV. Such results differ from the predictions of high-level G 0 W 0 , GW 0 and Koopmans-complaint hybrid calculations, which have been reported for monolayer GaSe and overestimate the band gap and the exciton binding en-ergy [26,62,73,74]. Using the experimental results and since both GaS and GaSe exhibit similar structural and electronic properties, it is safe to assume that the exciton binding energy of both monolayers would be of the same order and we assume a fundamental band gap of 3.6 eV for the GaS case. When compared with the experimental values, the results presented in this paper using both PBE-D and HSE06 underestimate the band gaps of the monolayers under study. However, HSE06 predictions are never more than 0.5 eV away from the experimental value and, having in mind that the fundamental gap decreases when the dielectric constant of the environment increases, its predictions should be closer to experiments involving substrates. In the case of PBE-D, it considerably underestimates the band gap and in principle cannot be used to draw conclusions on the electrical properties of defects in these monolayers. However, as shown in Sec. III A, when aligned with respect to vacuum the levels predicted using PBE-D are in line with those predicted by HSE06. Furthermore, due to the reduction of the band gap as temperature rises, results obtained using PBE-D become more appropiate when trying to draw conclusions for finite temperatures. These two features combined, make it possible to use PBE-D to study the formation energies of defects at relevant growth conditions (see Sec. III B).
Finally, in order to assess the spurious interaction between defects and their periodic images (see Sec. II B), static dielectric constants were calculated using the method presented in Ref. [75], applying either density functional perturbation theory (PBE-D) or self-consistent response to finite electric fields (HSE06) to a × a × na (n = 2, 4, 8, 16) supercells, with a being the in-plane lattice constants, and analyzing the convergence of the obtained results (see the supplementary material). The obtained values for both in-plane and out-of-plane components of the dielectric constant are also reported in Table I. Unfortunately, there are no experimental values available for the dielectric properties of the monolayers under study. There are, however, previous calculations for GaSe. [62,74] Surprinsingly, as the authors of Ref. [62] themselves recognize, their conclusion is that ε ⊥ of the monolayer is close to the value of the bulk, which is clearly unphysical.

B. Defect formation energies and transition levels
Defects were studied using 6 × 6 × 6 supercells with a 3 × 3 × 1 k -point grid. This supercell size implies a large enough vacuum of ∼22Å between monolayers. All other parameters used are the converged ones described in the previous section. In order to avoid local minima and to find accurate relaxed atomic configurations for the defects, the symmetry of the initial defective supercells was broken by introducing small random displacements to all atoms.
As mentioned in the introduction, one of the aims of this work is to produce a complete study of the native defects in which realistic growth conditions are taken into account. To that end, the approach described in Refs. [68,69,71,79] and, specifically, the methodology proposed in Ref. [63] for monolayer MoS 2 were followed.
The key quantity in a defect study where temperature effects are accounted for, is the Gibbs free energy of formation G f . For a given defect D in charge state q, G f is defined as and F [bulk] are the Helmholtz free energies of the defective and pristine supercells, respectively, V f is the formation volume of the defect, n i is an integer that indicates the number of atoms of type i added or removed to create defect D and µ i is the chemical potential of the same added or removed atoms. The fifth term in Eq. 1 accounts for the energy needed to charge the system. Consequently, it contains the valence band maximum (VBM) energy, E VBM , and the Fermi energy referenced to the VBM, E F . The final term, E corr , accounts for the correction needed to eliminate the spurious electrostatic interactions between periodic images of the defects under study. At T = 0 K, the Gibbs free energy of formation reduces to the formation energy, G f → E f . Following Ref. [63], here it is assumed that the F [D q ] − F [bulk] difference in Eq. 1 cancels out most of the vibrational contributions to the Helmholtz free of the system. Furthermore, due to the large gap of monolayer GaS and GaSe, it is also assumed that finite temperature effects on the electronic contribution to the free energies are negligible. This means that F [D q ] and F [bulk] can be replaced by E 0 [D q ] and E 0 [bulk], i.e., the corresponding total energies as obtained from DFT calculations at 0 K. Regarding the defect formation volume, it is assumed to be zero as the concept is not valid for an atomically-thin material.
As it is clear from the statement expressed in Eq. 1, it is assumed that the pressure and temperature dependence of G f enters via the chemical potentials, which define the environmental conditions. Naturally, in order to produce results comparable with experiments, one is interested in realistic limits for the chemical potentials. Thus, it is assumed that one limiting condition for the chemical potentials of Ga, S and Se is that they are in equilibrium with the corresponding monolayer for all temperatures and pressures considered in our study, In addition, a further limiting condition is that µ i cannot be lower than the corresponding values for the stable elemental phases, µ ref i . In other words, this second limiting condition means that the differences ∆µ i = µ i − Above or below the lines that limit these areas, species forming the monolayers prefer their elemental phases for the corresponding temperature. Solid red and yellow lines correspond to the reference curves based on thermochemical tables [76][77][78]. Partial pressures are given in bar.
For the actual calculation of µ i , it should be kept in mind that for elemental phases the chemical potentials are proportional to the corresponding total Gibbs free energy, where N ref is the number of particles in the unit cell of interest. At temperature T and reference pressure P • , the free energy is given by where E vib is the energy of zero-point vibrations, directly obtained from DFT results, and ∆H and ∆S account for the changes in enthalpy and entropy, respectively, associated to finite temperatures. Values for ∆H and ∆S at the standard pressure of P • = 100 kPa are readily available from thermochemical tables. [76][77][78] E 0 for Ga, S and Se, was calculated using the fact that at T = 0 all these elements prefer solid phases. Specifically, Ga stabilizes in an orthorhombic structure, and both S and Se crystallize in monoclinic allotropes formed by S 8 and Se 8 rings. After obtaining all the needed values, it is possible to use Eq. 5 to focus on realistic growth conditions. As mentioned in the Sec. I, the most succesful approaches to obtain large-area highquality monolayers of GaS and GaSe are MBE [47][48][49], and VPT [31,33,[50][51][52], which usually imply growth temperatures of about 700 K for the former and 1000 K for the latter. Therefore, when discussing Gibbs formation energies of defects in monolayer GaS and GaSe, it will always refer to these two temperatures.
On the one hand, since Ga melts at 302.9 K, it is expected to be liquid at the relevant growth temperatures and we can directly apply Eq. 5 using thermochemical tables for ∆H and ∆S. [77] On the other hand, at these same temperatures both S and Se are expected to be found in their diatomic gas phase [76][77][78]. As a result, Eq. 5 cannot be directly used. Following Ref. [63], assuming ideal gas behavior, the chemical potential of S 2 and Se 2 can be written as where X = S, Se and P X2 is the partial pressure of the corresponding diatomic gas. E X2 is the total energy of the corresponding isolated diatomic molecules and E vib X2 (0K) is their energy of zero-point vibrations. Both E X2 and E vib X2 (0K) are directly obtained from DFT calculations. In the case of ∆H and ∆S, as done for Ga, a P • = 100 kPa is assumed and thermochemical tables are used [76][77][78].
The resulting limits for the chemical potentials in the monolayers under study and using PBE-D are shown in Fig.1, where the white areas are the regions for which the monolayers are stable. Above or below the lines that limit these areas Ga, S and Se prefer their elemental phases at the corresponding temperature. Specifically, if the Garich limit is exceeded, Ga will form precipitates and if Sor Se-rich limits are the ones exceeded, there will be S and Se desorption from the monolayer. In all cases, a rich limit means that the condition µ i = µ ref i has been achieved for one of the species in the monolayer. For example, in Fig.1(a), µ Ga = µ ref Ga in the Ga-rich limit and, using Eq. 2, µ Ga = µ GaS − µ ref S in the S-rich limit. A further detail shown in Fig.1, are solid red and yellow lines which correspond to the tabulated values for the lowest energy phases [76][77][78] and serve as references. At this point it is important to remark that an often overlooked issue when dealing with small molecules, e.g. diatomic S or Se, is that depending on the used computational parameters one can have relatively large errors on their calculated energies. This naturally affects the obtained chemical potentials shown in Fig.1. However, what is crucial for defect calculations is to be consistent with respect to the choice of computational methods used to obtain all the needed chemical potentials in a given study [68].
Being aware that realistic growth temperatures for the monolayers of GaS and GaSe are between 700 K and 1000 K, it is possible to see in Fig.1 that the reduction induced in the chemical potentials due to temperature effects is considerable. Thus, it is not correct to use µ i (T = 0 K) to obtain reliable Gibbs formation energies of defects at relevant finite temperature. Moreover, Fig.1 also shows that the actual values of the chemical potentials are also dependent on the S and Se partial pressures. Throughout this paper, Eq. 1 and PBE-D are applied to calculate defect Gibbs formation energies at 700 K and 1000 K, using the values for µ i shown in Fig.1 for S and Se partial pressures of 1 bar. The reason why PBE-D is used for finite temperature calculations is explained in Sec. III B.
After having determined accurate values for µ i , only the E corr term has to yet be calculated in order to take full advantage of Eq. 1. As mentioned before, this term is needed when using the supercell approach and takes care of correcting the spurious electrostatic interactions due to defect periodic images immersed in a jellium background. The existence and relevance of these spurious effect is known since the 1985 seminal work of Leslie and Gillan [81]. After a great deal of effort in this direction, currently there are several reliable schemes to calculate E corr for bulk materials [82][83][84][85]. Among those, the one proposed by Freysoldt, Neugebauer and Van de Walle (FNV) in Ref. 85 is rigorous and has been proven to be the most accurate option [86,87]. Since there is less screening in the case of a supercell containing a defective monolayer, the spurious electrostatic interactions between charged defects in such systems are expected to have a greater impact. This is the reason why during the last years there has been interest in finding accurate schemes that include the inhomogeneity of the dielectric environment natural to monolayers [88][89][90][91][92][93][94]. Unfortunately, contrary to the bulk case, there is still no detailed study of how these methods compare to each other. However, based on the success of the FNV method for bulk systems and the rigorous method behind it, its extension for two-dimensional materials (FNV-2D) is used in this paper [94]. The for- mation energy of the negatively charged S vacancy, V −1 S , in monolayer GaS at different supercell sizes is shown in Fig. 2. For each supercell size used, both the uncorrected and corrected values are presented. As can be seen there, for this very localized defect the FNV-2D is able to deliver accurate predictions even for the smallest case included in our study, i.e, the 4 × 4 × 4 supercell. Although using such small simulation box would reduce the computational cost of our study, it is important not to forget that we must also use a supercell size that guarantees the correct localization of the defects states. This is the reason why in this paper 6 × 6 × 6 supercells were used: they allow the calculation of formation energies values that are less than 0.1 eV away from dilute limit and are large enough to contain even the most delocalized defects under study, i.e., the interstitials. As an example, in Fig. 3 we present the charge density difference between the neutral and positively charged Ga interstitial in GaS in a 6 × 6 × 6 supercell, where it is clear that the hole is correctly localized within the boundaries of the used supercell.
A further defect characterization tool are charge transition levels (CTLs), defined as the Fermi level for which the formation energies of charge states q and q are equal. They must not be confused with defect Kohn-Sham states, which are called states throughout this paper. Charge transition levels for defect D at T = 0 K can be obtained as follows: The transition levels are relevant because they can be directly related to experiments in which the defects are able to fully relax after the charge transition [68].

III. RESULTS
In order to understand the properties of native defects in the monolayers under study, we start the presentation of our results by discussing their CLTs together with their basic geometries. With this first part we aim at understanding whether native defects are a possible explanation to the experimentally observed intrinsic doping on the monolayers. Since it is not central to our goal, we only discuss the electronic structure of the studied defects in a brief manner. Based on such results, in a second part we present the formation energies of the same relevant native defects at realistic growth temperatures. Our aim there is to reveal whether harmful defects are likely to form when growing monolayers using state-of-the-art methods. More importantly, in this second section we propose growth conditions that would limit the formation of the most detrimental native defects.

A. Transition levels and intrinsic doping
It may be tempting to directly compare calculated Kohn-Sham (KS) states for a given system with experiments. But it is important to keep in mind that KS states, even the ones obatined with hybrid functionals, cannot be directly associated to any experimentally significant levels [68]. The relevance of CTLs to characterize defects in semiconductors lies on the fact that they are, as mentioned in Sec. II B, one of the quantities which can be predicted via ab-initio methods and are amenable to comparison with experiment. This interesting feature of CTLs is rooted on the fact that they are calculated using Eq. 7 and total energies, which are well described by DFT. Furthermore, although defects can go through excitations in which their charge state does not change, it is more common for them to trap charge carriers. This is the process described by the CTLs.
In Fig. 4 we present the relaxed geometries of the studied native defects in monolayer GaS and in their charge neutral state. We show optimized geometries as obtained with PBE-D, which delivers a better description for out-of-plane defects, i.e. adatoms, and we do not find considerable differences with HSE results for defects inside the monolayer. Since relevant features of the optimized geometries of the native defects are common to both monolayer GaS and GaSe, we only show the corresponding GaSe configurations in the Supplementary Material.
In addition, we present the vacuum-aligned transition levels of all relevant defects for both monolayers in Fig. 5. For most defects we present the obtained CTLs using both PBE-D and HSE06. However, due to the previously discussed fact that HSE06 does not account for van der Waals interactions, for adatoms we only show the PBE-D results. Naturally, due to its more accurate description of the band gap, CTLs predicted using HSE06 are expected to be more accurate at T = 0 K and are the ones we use to draw conclusions on the electrical properties of the studied defects. Nevertheless, by showing both PBE-D and HSE06 results in the same figure, we are able to observe that both functionals provide a consistent picture of the electrical properties of the defects and predict CTLs that agree well when given with respect to vacuum. This first general conclusion is fundamental to explain why we use PBE-D for finite-temperature calculations, as we argue in Sec. III B.

Gallium interstitial
The optimized atomic configuration of the neutral Ga interstitial, Ga i , is shown in Fig. 4(a). As we can see there, at this charge state such interstitial has only a marginal effect on the positions of its surrounding atoms. Since Ga i is located in the middle of the monolayer, its optimized geometry indicates that the D 3h symmetry of the pristine monolayer is retained. The electronic structure of the neutral Ga i features three occupied defect states within the band gap, deep in the case of the GaS monolayer and close to the conduction band in the GaSe system. These findings suggest that the Ga i may be prone to donate electrons, which is actually the case. As we can see in Figs. 5(a) and 5(b), with respect to the VBM this defect has (+1/0) = 2.49 eV and (+1/0) = 2.31 eV for GaS and GaSe, respectively. This means that for both materials the Ga i is a deep donor, with ionization energies larger than 0.65 eV. Furthermore, also for both materials, as the defect donates the electron it exhibits a Jahn-Teller distortion and the symmetry is broken to C s .
FIG. 6. Formation energies of native defects in monolayer GaS at 700 K for the (a) Ga-rich limit and the (b) S-rich limit, and at 1000 K for the (c) Ga-rich limit and the (d) S-rich limit. Results obtained with PBE-D and not given with respect to vacuum.

Antisites
Contrary to the case of the Ga i , the antisites in monolayer GaS and GaSe induce significant changes in the positions of their surrounding atoms. Their relaxed configurations in the neutral state and for the GaS system are shown in Fig. 4(b) and Fig. 4(c). For the corresponding relaxed structures in monolayer GaSe please refer to the supplementary material. For both monolayers and for all relevant charge states the antisites reduce the D 3h symmetry to a C 3v case. We do not find extra symmetry breaking when these defects become charged. As they trap or donate electrons, the antisites in the monolayers under study induce localized, but symmetric, displacements in their surrounding atoms.
In the specific cases of the Ga S and Ga Se , their electronic structure is characterized by three defects states in the gap, one occupied and two unoccupied. Such features suggest that these antisites can, in principle, be amphoteric. The obtained CTLs for Ga S and Ga Se are presented in Figs. 5(a) and 5(b), and detailed values of the CTLs are summarized in Table II. As we can see in both the figures and the table, the amphoteric character of the Ga S and Ga Se is confirmed by the calculated CTLs. In addition, based in the position of the CTLs, we conclude that these antisites can behave as both deep acceptors or deep donors.
Regarding the S Ga and Se Ga antisites, their electronic structure is dominated by three defect states in the gap, one occupied and very close to the VBM and two unoc-cupied and deep in the gap. Due to these features, as in the case of Ga S and Ga Se , one would expect that both S Ga and Se Ga antisites have an amphoteric behavior. The calculated CTLs, shown in Figs. 5(a) and 5(b) and summarized in Table II, confirm their amphoteric character. The main difference with respect to the case of Ga S and Ga Se , is that both S Ga and Se Ga exhibit also an extremely deep (−1/ − 2), that would be stable only under strong n-type doping. Based on the position of the CTLs, we can predict that both S Ga and Se Ga behave as deep acceptors or deep donors.  Formation energies of native defects in monolayer GaSe at 700 K for the (a) Ga-rich limit and the (b) Se-rich limit, and at 1000 K for the (c) Ga-rich limit and the (d) Se-rich limit.

Chalcogenide interstitials
Analogous to the case of the Ga i and as can be seen in Fig. 4(d), the neutral chalcogenide interstitials S i and Se i do not induce large changes in the positions of their surrounding atoms and retain the D 3h symmetry of the pristine monolayers. However, contrary to what we found for Ga i , associated to the charging process of both S and Se interstitials, we only observe localized but symmetric distortions, i.e., the D 3h symmetry is kept even in their charged states. Regarding their electronic structure, the chalcogenide interstitials feature three defect states in the band gap, two of which are occupied. Among these defects states, only one, induced by Se i , is located close to the VBM and the rest are deep within the band gap. These results suggest again an amphoteric tendency, which is confirmed by the calculated CTLs shown in Figs. 5(a) and 5(b) and summarized in Table III. Based on the position of the CTLs, we can predict that both S i and Se i behave as deep acceptors or deep donors.

Vacancies
In the case of the gallium vacancy V Ga , we do observe large relaxations associated to the presence of the defect, Fig. 4(e). Specifically, the neighbouring gallium atom relaxes to the middle of the monolayer. Due to this change, the neutral V Ga retains the D 3h symmetry of the nondefective monolayers. Such symmetric configuration is  also observed when the V Ga traps an extra electron, a process that only induces a symmetric outward "breathing" distortion affecting the atoms surrounding the vacancy. The electronic structure of this defect is characterized by the appearance of three defect states in the band gap: one just above the VBM, occupied (unoccupied) in its spin-up (spin-down) channel, and two unoccupied and deep in the gap. In principle, with such features one would expect a tendency towards trapping extra electrons for the V Ga . As we can see in Figs. 5(a) and 5(b), and in Table IV this defect is a deep acceptor with ionization energies larger than 0.7 eV for both monolayers. Chalcogenide vacancies, V S and V Se , break the symmetry of the system to C 3v as can be seen in Fig. 4(f). Although relaxations are observed when these vacancies change their charge state, the symmetry is not lowered further. For S and Se vacancies the electronic structure is dominated by three defect states, all of which are deep in the gap and only one is occupied. In principle, these anti-sites could then be either positively or negatively charged. However, the obtained CTLs presented Figs. 5(a) and 5(b), and in Table IV show that such charging processes would happen only under extreme conditions. Our assesment is that V S and V S are extremely deep amphoteric defects.

Adatoms
As mentioned at the beginning of this section, conclusions regarding the electrical properties of intrinsic defects in monolayer GaS and GaSe should be obtained from CTLs calculated using hybrids. Unfortunately, since HSE06 does not account for van der Waals interactions, we cannot use it to correctly predict the relaxed configuration of out-of-plane defects. Thus, it cannot deliver reliable results for adatoms. This is the reason why for these defects, Figs. 5(a) and 5(b) only show PBE-D results. In the case of gallium adatoms, whose relaxed configuration is directly on top of the center of the hexagons as shown in Fig. 4(g), PBE-D predicts (+1/0) = 1.88 eV and (+1/0) = 1.67 eV for monolayer GaS and GaSe, respectively. In addition, PBE-D predicts that S and Se adatoms prefer to be located on top of chalcogenide sites, see Fig. 4(h), and are always neutral. The key to draw relevant conclusions regarding the electrical activity of adatoms from these PBE-D results is that, despite suffering from the band gap problem, CTLs obtained using PBE-D agree well with HSE06 ones when aligned with respect to vacuum. This means that opening the gap would only result in gallium adatoms being deep donors, even deeper compared to what PBE-D predicts, and the chalcogenide adatoms are expected to be electrically inactive as they are neutral for all doping conditions.

Native defects and intrinsic doping
The central goal of our research is to shed light on the currently unclear origin and nature of intrinsic doping in monolayer GaS and GaSe. Experimental studies have measured shifts in the Fermi level and have proposed defects as a possible explanation for such behavior [22,49,57]. However, our results do not agree with this hypothesis. As shown in Figs. 5(a) and 5(b), no native defect exhibits a shallow nature. Specifically, in the monolayers under study all ionization energies are above 0.6 eV. This means that native defects in these systems cannot be easily ionized and are not efficient dopants. Therefore, we conclude that the intrinsic doping observed in experiments should be caused by charge transfer with the substrates or by the presence of shallow extrinsic dopants.
The fact that native defects cannot explain intrinsic doping in monolayer GaS and GaSe does not mean that they are irrelevant. Actually, it is quite the opposite situation. Being deep defects, they can act as carrier traps and degrade the efficiency of any device which uses these monolayers. In addition, being stable in different charge states, the studied native defects can act as compensating centers and, as a consequence, make it difficult to engineer the doping conditions of the monolayers. We dedicate the next section to elucidate which growth conditions can lower the concentrations of harmful defects.

B. Formation energies at realistic growth conditions
We now use the methodology described in Sec. II B and tackle the calculation of formation energies at finite temperatures. To that end, we keep in mind that high-quality monolayers of GaS and GaSe are grown at about 700 K via MBE [47][48][49], and at about 1000 K via VPT [31,33,[50][51][52]. Thus, in this section we present our results for these two experimentally relevant temperatures. In Sec. II B we clearly explained that in our approach, based on Ref. [63], the temperature dependence entered via the chemical potentials. However, to be closer to reality, the temperature dependence of the formation energies should also enter via changes in the band gap, which reduces as temperature rises. For bulk GaS and GaSe it has been reported that the band gap decreases by about 0.45 meV/K and 0.51 meV/K, respectively [95,96]. Assuming a rate of the same magnitude for the monolayers means that at 1000 K the experimental band gap would be reduced to E g = 3.15 eV for monolayer GaS and to E g = 2.99 eV for monolayer GaSe. Actually, assuming an analogous behavior to the one observed for MoS 2 , for which the rate of change of the band gap increases from 0.4 meV/K in bulk [97,98] to 0.7 meV/K in the monolayers [99], one would expect that the experimental band gaps of monolayer GaS and GaSe could be reduced further to E g 2.8 eV for monolayer GaS and to E g 2.6 eV, respectively. Based on the values reported in Table I and keeping in mind that the fundamental gap decreases due to the presence of substrates, we can safely say that at high temperatures PBE-D should deliver band gaps, defect formation energies and CTLs more in line with experiments. In addition, since it takes van der Waals interactions into account, PBE-D allows us to correctly describe the physics of adatoms. Due to these two relevant advantages, we adopted PBE-D for our finite-temperature calculations. The resulting formation energies for both monolayers at 700 K and 1000 K are presented in Figs. 6 and 7. In there, chemical potentials are obtained using the methodology described in Sec. II B and from Fig. 1. Ideally, one would use the finite temperature formation energies and translate them to concentrations [68]. The usual procedure to do so implies the self-consistent calculation of the intrinsic Fermi level [63,68]. Although simple and clear, this method cannot be directly applied to monolayers due to the strong influence of substrates, i.e., the Fermi level in the monolayer is likely to be far from the intrinsic one. Therefore, in the following, we use the calculated formation energies and differences between them to qualitatively draw conclusions about the relative abundance of specific defects.
For Ga-rich conditions at 700 K, Figs. 6(a) and 7(a), the native defects that are most likely to occur are chalcogenide vacancies and gallium adatoms. A little bit less common, but expected to occur in considerable quantities, are gallium vacancies and the Ga S /Ga Se antisites. For the same Ga-rich conditions but at 1000 K, Figs. 6(c) and 7(c), we predict a more or less equal defect landscape, with the only difference of slightly larger formation energies for chalcogenide vacancies, gallium adatoms and Ga S /Ga Se antisites, and slightly lower formation energies for gallium vacancies. The consequence of these findings can be understood based on the discussion presented in Sec. III A. On the one hand, S and Se vacancies are expected to be electrically harmless. On the other hand, gallium adatoms, gallium vacancies and the Ga S /Ga Se antisites exhibit deep CTLs and are expected to be harmful. Specifically, Ga S /Ga Se have relatively high formation energies and are not expected to play a relevant role as compensation center. However, they do feature a deep CTLs very close to the center of the band gap in both monolayers and are expected to be very efficient charge carrier traps. In the case of gallium adatoms (deep donors) and gallium vacancies (deep acceptors), our results indicate that both are expected to be efficient compensation centers. Naturally, under Ga-rich conditions, gallium adatoms are more likely to occur than gallium vacancies, causing an asymmetric behavior with respect to doping. In other words, when grown under Ga-rich conditions, monolayer GaS and GaSe are easier to dope n-type than p-type.
In the case of chalcogenide-rich conditions at 700 K, Figs. 6(b) and 7(b), the native defects that are most likely to occur are S and Se adatoms, gallium vacancies and gallium adatoms. Due to their slightly higher, but still low, formation energies we predict considerable amounts of chalcogenide vacancies and S Ga /Se Ga antisites. As expected, at 1000 K and the same chalcogenide-rich conditions, Figs. 6(d) and 7(d), we predict an essentially identical defect landscape. At such temperature chalcogenide adatoms, gallium vacancies and S Ga /Se Ga antisites exhibit slightly larger formation energies, and gallium adatoms have slightly lower formation energies. Based on the discussion presented in Sec. III A we know that S and Se adatoms are expected to be electrically harmless and that gallium adatoms, gallium vacancies and S Ga /Se Ga antisites induce deep CTLs and are expected to be harmful. Analogous to what we found for Ga S /Ga Se antisites in Ga-rich conditions, when samples are grown under chalcogenide-rich conditions, S Ga /Se Ga antisites have relatively high formation energies and are not expected to be efficient compensation centers but very efficient charge carrier traps. Regarding compensation in samples grown under chalcogenide-rich conditions, we predict a strongly asymmetric process driven by the presence of gallium adatoms (deep donors) and gallium vacancies (deep acceptors), and heavily dominated by the latter. Thus, our results show that monolayers GaS and GaSe grown under chalcogenide-rich conditions will be far easier to dope p-type than n-type. The reason for the strong asymmetry is caused by the large difference in formation energies of the two defect types.
At this point we are certain that, although some native defects are harmless, most of them are harmful and should be avoided. To this end we can use a strong tendency hidden in Figs. 6 and 7: harmful defects observed at Ga-rich conditions have higher formation energies in chalcogenide-rich conditions, and vice versa. Therefore, we propose that growing monolayer GaS and GaSe at intermediate conditions should increase (decrease) the formation energies (concentration) of all harmful defects at the same time. As an example of the approach, in Fig. 8 we show the formation energies of native defects in monolayer GaSe at 700 k for intermediate growth conditions, µ Ga = −3.84 eV and µ Se = −4.3 eV (which are certainly within the stability region of the monolayer, see Fig. 1). In there we can see how such small change in the chemical potentials leads to a situation where the formation energies of all harmful defects (gallium/chalcogenide vacancies and antisites) are simultaneously increased above ∼ 1 eV. Since most high-quality monolayer GaS and GaSe grown via MBE or VPT methods are actually grown in chalcogenide-rich conditions, our proposal to experimen-talists is to slightly increase the gallium concentration or decrease the S/Se partial pressures.

IV. CONCLUSIONS
We have performed first-principles calculations to assess the electrical and thermodynamic properties of native defects in GaS and GaSe monolayers. To that end we calculated their charge transition levels and formation energies, using state-of-the-art schemes to account for finite-size effects and the temperature dependence of the chemical potentials. Based on the position of the obtained charge transition levels, we conclude that among the studied native defects there are no shallow dopants and the experimentally observed intrinsic doping should then be caused by charge transfer to/from the substrates or by the presence of shallow extrinsic dopants. In addition, the position of the charge transition levels also allowed us to predict that most native defects are efficient compensation and recombination centers. Thus, native defects should be avoided in order to grow technologically reliable GaS and GaSe monolayers. By studying relevant formation energies at realistic growth temperatures, we are able to prove that by growing GaS and GaSe monolayers at intermediate conditions, far from the chalcogenide-or Ga-rich regimes, it should be possible to simultaneously decrease the concentration of all relevant harmful native defects. In other words, when growing monolayer GaS and GaSe under S/Se-rich conditions, our findings translate to a simple experimental directive: to slightly increase the gallium concentration or to slightly decrease the S/Se partial pressures.