Native point defects from stoichiometry-linked chemical potentials in cubic boron arsenide

TL;DR

This paper develops a framework linking crystal stoichiometry to chemical potentials, improving defect formation energy calculations in semiconductors like cubic boron arsenide.

Contribution

It introduces a canonical approach that relates chemical potentials directly to composition, avoiding nonphysical results from ad hoc assumptions.

Findings

Antisites are the main defects in off-stoichiometric BAs.

BAs antisites and vacancies dominate in stoichiometric BAs.

The new framework ensures positive formation energies across Fermi levels.

Abstract

The presence of a point defect typically breaks the stoichiometry in a semiconductor. For example, a vacancy on an A-site in an AB compound makes the crystal B-rich. As the stoichiometry changes, so do the chemical potentials. While the prevalent first-principles methods have provided significant insight into characters of point defects in a transparent manner, the crucial connection between crystal stoichiometry and chemical potentials is usually not made. However, ad hoc choices for chemical potentials can lead to nonphysical negative formation energies in some Fermi level ranges, along with questions about charge balance. Herein, we formulate a canonical framework describing how the chemical potential of each element is directly linked to the composition of the crystal under (off-)stoichiometric conditions instead of the ad hoc assumption that the chemical potential is the elemental limit under a certain growth condition. Consequently, the chemical potential changes with the Fermi level within the band gap, and the formation energies are positive. Using such an approach, we present $ab$ $initio$ results for native point defects in BAs, a semiconductor with ultra-high room temperature thermal conductivity. We find that antisites are the constitutional defects in off-stoichiometric material, while B$_\mathrm{As}$ antisites and B vacancies dominate in the stoichiometric material. We further discuss the thermodynamic equilibrium and charge neutrality point in BAs in light of our stoichiometry-determined chemical potentials. As discussed, our work offers a more applicable and accessible approach to tackle defect formation energies in semiconductors, especially the ones with wide gap where negative formation energies are commonly seen.

Tables

Table 1. Table 1: Summary of zero-temperature formation energies E f subscript 𝐸 𝑓 E_{f} and Arrhenius-fitted formation enthalpies H f subscript 𝐻 𝑓 H_{f} at mid gap for native point defects in cubic BAs, along with their relative differences Δ ​ E rel Δ subscript 𝐸 rel \Delta E_{\mathrm{rel}} and Arrhenius-fitted formation entropies S f subscript 𝑆 𝑓 S_{f} .

Defects		$E_f\left(\mathrm{eV}\right)$	$H_f\left(\mathrm{eV}\right)$	$\Delta E_{\mathrm{rel}}$(%)	$S_f$($k_B$)
	$V_{\mathrm{As}}^0$	4.79	4.773$\pm$0.001	0.42	1.39
	$V_B^0$	3.36	3.628$\pm$0.003	8.04	11.8
B-rich	$B_{\mathrm{As}}^{-1}$	0	0	0	$-$
	${\mathrm{As}}_B^{2+}$	5.36	5.702$\pm$0.004	6.74	10.8
	$B_i^{+}$	4.20	4.312$\pm$0.01	2.62	0.426
	${\mathrm{As}}_i^0$	8.45	8.635$\pm$0.001	2.25	10.5
	$V_{\mathrm{As}}^0$	7.47	7.624$\pm$0.001	2.01	6.79
	$V_B^0$	0.68	0.7763$\pm$0.0014	14.1	6.46
As-rich	$B_{\mathrm{As}}^{-1}$	5.35	5.702$\pm$0.004	6.54	10.08
	${\mathrm{As}}_B^{2+}$	0	0	0	$-$
	$B_i^{+}$	6.88	7.163$\pm$0.011	4.07	5.83
	${\mathrm{As}}_i^0$	5.77	5.784$\pm$0.001	0.17	5.14
	$V_{\mathrm{As}}^0$	5.98	5.982$\pm$0.001	0	5.12
	$V_B^0$	2.37	2.418$\pm$0.003	2.04	8.15
Stoichi-	$B_{\mathrm{As}}^{-1}$	2.37	2.418$\pm$0.003	2.04	7.45
ometric	${\mathrm{As}}_B^{2+}$	3.04	3.464$\pm$0.018	13.9	2.55
	$B_i^{+}$	5.29	5.521$\pm$0.01	4.35	4.15
	${\mathrm{As}}_i^0$	7.50	7.442$\pm$0.001	0.8	5.05

Equations

$${\color[rgb]{0,0,0}E_{f}(\mathrm{V_{A}})=E\mathrm{(V_{A})}-E_{\rm perf}+\mu_{\rm A}+Q(E_{F}+E_{\mathrm{v}})},$$

$${\mu }_{\mathrm{A}}+{\mu }_{\mathrm{B}}=E(\mathrm{AB}),$$

$$C_X(T)=C_{\mathrm{s}}{\theta }_X\exp \left[-F_X^f/\left(k_{\mathrm{B}}T\right)\right],$$

$$\frac{N_{\mathrm{As}}}{V}$$

$$\frac{N_{\mathrm{B}}}{V}$$

$$C_{{\mathrm{V}}_{\mathrm{As}}}+C_{{\mathrm{B}}_{\mathrm{i}}}+C_{{\mathrm{B}}_{\mathrm{As}}}-C_{{\mathrm{V}}_{\mathrm{B}}}-C_{\mathrm{A}{\mathrm{s}}_{\mathrm{i}}}-C_{\mathrm{A}{\mathrm{s}}_{\mathrm{B}}}=2\delta \frac{N}{V}.$$

$$\Delta n_D(E)=n_D(E)-p_D(E)=\sum _X{Q}_X{C}_X$$

$$E_{F_D}=E_{F_i}+k_{B}T\ln \left(\frac{\Delta n_D+n_i}{n_i}\right)$$

$$m_{n,\mathrm{dos}}^{\ast }=6^{2/3}\times (m_l{m}_t{m}_t{)}^{1/3}=2.19m_0.$$

$$E_{\mathrm{i}}$$

$$n_i(T)=\sqrt{N_c{N}_v}\exp \left(-\frac{E_g}{2k_{B}T}\right).$$

$$N_c(T)$$

$$N_v(T)$$

$$n_i(T)=8.85\times 10^{15}{\left(\frac{T}{\mathrm{K}}\right)}^{3/2}\exp \left(-\frac{9351}{T/\mathrm{K}}\right){\mathrm{cm}}^{-3}.$$

$$S_f^{\mathrm{conf}}=k_{\mathrm{B}}[c-c\ln (c)+c\ln ({\theta }_X)],$$

$$F(p,T)=H(p)+\frac{1}{4\pi }\int g\left[\omega (p)\right]h\omega (p)d\omega +k_{B}T\int g\left[\omega (p)\right]\ln \left[1-\exp \left(-\frac{h\omega (p)}{2\pi k_{B}T}\right)\right]d\omega$$

$$F_{\mathrm{As},\mathrm{vib}}(T)=0.0252-7.7436\times 10^{-7}{T}^2+3.9500\times 10^{-10}{T}^3-8.392\times 10^{-14}{T}^4$$

$$F_{\mathrm{B},\mathrm{vib}}(T)=0.1336-9.3379\times 10^{-8}{T}^2-6.4218\times 10^{-11}{T}^3+2.4277\times 10^{-14}{T}^4$$

$$F_{\mathrm{BAs},\mathrm{vib},512}(T)=2\left(0.0828-3.1397\times 10^{-7}{T}^2+7.1954\times 10^{-11}{T}^3-6.0088\times 10^{-15}{T}^4\right)$$

$$F_{\mathrm{BAs},\mathrm{vib},64}(T)=2\left(0.0820-3.2050\times 10^{-7}{T}^2+7.8671\times 10^{-11}{T}^3-7.8746\times 10^{-15}{T}^4\right)$$

$$\frac{N_{\mathrm{As}}}{V}$$

$$\frac{N_{\mathrm{B}}}{V}$$

$$C_{V_{As}}+C_{B_i}+C_{B_{As}}-C_{V_B}-C_{A{s}_i}-C_{A{s}_B}=2\delta \frac{N}{V}.$$

$$m{\mu }_{\mathrm{A}}+n{\mu }_{\mathrm{B}}+l{\mu }_{\mathrm{C}}$$

$$C_{{\mathrm{V}}_{\mathrm{A}}}+C_{{\mathrm{B}}_{\mathrm{i}}}+C_{{\mathrm{B}}_{\mathrm{A}}}-C_{{\mathrm{V}}_{\mathrm{B}}}-C_{{\mathrm{A}}_{\mathrm{i}}}-C_{{\mathrm{A}}_{\mathrm{B}}}$$

$$C_{{\mathrm{V}}_{\mathrm{A}}}+C_{{\mathrm{C}}_{\mathrm{i}}}+C_{{\mathrm{C}}_{\mathrm{A}}}-C_{{\mathrm{V}}_{\mathrm{C}}}-C_{{\mathrm{A}}_{\mathrm{i}}}-C_{{\mathrm{A}}_{\mathrm{C}}}$$

$$C_{{\mathrm{V}}_{\mathrm{B}}}+C_{{\mathrm{C}}_{\mathrm{i}}}+C_{{\mathrm{C}}_{\mathrm{B}}}-C_{{\mathrm{V}}_{\mathrm{C}}}-C_{{\mathrm{B}}_{\mathrm{i}}}-C_{{\mathrm{B}}_{\mathrm{C}}}$$

Full text

I Introduction

Point defects are the only thermodynamically stable crystal defects at finite temperatures. For functional, optical, and electronic materials, point defects influence a material’s properties arguably more than any other factor. While defects or impurities can serve as dopants to achieve either $p$- or $n$-type conduction, they are mostly detrimental to the precise control of optical and transport properties. As a result, extensive studies from both experiment and theory have been put forward to better understand the nature of both native point defects and impurities, and thus improve control of their formation in semiconductors. While experimental observation of point defects is especially difficult and typically limited to indirect probing,JOHNSON201717; stoddard_duscher_windl_rozgonyi_2004; Asel18 their study has been one of the central points of computational materials science. Especially first principles, or ab initio calculations have been shown to be a powerful tool to identify the characteristics of point defects in functional materials, starting in the last two decades of the previous century such as reported in Refs. Scheffler88; laks1991role; Neuge94; windl1999first to pick a few. Several review chapters and articles have also been published, summarizing a gamut of topics from the calculation of point defect energetics and dynamics to their combination with process simulations.Estreicher00; Windl2004diffusion; Windl2006; Freysoldt2014first

Generally speaking, based on DFT total energies for defective and perfect cells, the defect formation energies can be calculated using (for the example of an A vacancy in a stoichiometric AB compound, VA)

[TABLE]

where Q is the charge state, $E_{F}$ the Fermi level with $0<E_{F}<E_{\rm gap}$, $E_{\rm v}$ the energy of the top of the valence band, $\mu_{\rm A}$ the chemical potential of species A, and $E_{\rm perf}$ the total energy of the perfect supercell. $E\mathrm{(V_{A})}$ is the total energy of a relaxed AB supercell with A vacancy. The same framework would hold if instead of total energies free energies were used. Although the DFT the energies in Eq. 1 can be corrected with a plenitude of schemes for e.g. finite size effects on energiesWindl2004diffusion and shallow dopant levels,Windl1998first artifacts from the density functional used,Rampi12 charge interactions between periodic images,Makov95 and other artifacts,Freysoldt2014first the understanding of the chemical potentials in Eq. 1 has yet to be well converged.

The currently prevalent heat of formation (HoF) method focuses on the ultimate limits where one of the chemical potentials is aligned with the elemental solids.Freysoldt2014first While this method provides a straightforward way to determine defect formation energies, it does not explore the intrinsic connection between the stoichiometry and the chemical potentials because of the binary choice of chemical potentials along with the intermediate and often unconsidered thermodynamic ensemble situated between grand canonical and canonical.Freysoldt2014first Without proper understanding, the HoF approach may have pitfalls that lead to defect energies and concentrations in unphysical limits, may leave charge neutrality typically unexplored, and may lead to focusing on the wrong prevalent point defects.

An aspect within work based on the HoF approach that can be easily overlooked is the connection between the Fermi level and charge neutrality point (CNP). Although explicitly defined recently in a reviewFreysoldt2014first after decades of HoF use, few publications based on HoF discuss this aspect, which when neglected makes it difficult to interpret the traditional HoF results in a physical way in a gapped material. More specifically, in HoF, which is usually explored at zero temperature for a given chemical potential, the slope in the formation-energy-versus-Fermi-level plots corresponds to the point defect’s charge. At zero temperature, charge neutrality, a fundamental requirement for point defects in any material with a band gap,Callister can be identified by finding the lowest-energy intersection point between the line for a point defect with positive slope and the line for a point defect with negative slope with the same absolute values of slope. There are no previous studies that address the CNP question from the HoF results, including discussing them in light of finite temperatures. When only considering the CNP, it can be several eV away from mid gap in large-gap materials such as Al2O3 as a prototypical example Choi13, while the defect formation energies are so high that the defect concentrations will not shift the Fermi level significantly from mid gap in the material (for details see SM Sec. S.1). When instead focusing on mid gap, the questions about charge neutrality and resulting electron and hole concentrations remain unanswered.

While thus the presently prevailing HoF method, whose approach was originally developed for compounds in contact with elemental solids,Scheffler88 focuses on the elemental limits of the constituent elements, the approach prposed by Zhang and Northrup in 1991Northrup91 focuses on the difference between the chemical potentials $\Delta\mu=\mu_{\rm A}-\mu_{\rm B}$ under the requirement that the AB compound is stable with respect to the elemental chemical potentials, with no explicit values for $\mu_{\rm A}$ and $\mu_{\rm B}$ assumed in the input or received in the output. This allows to connect the defect concentrations to the overall stoichiometry through a quantitative link to chemical potential and Fermi level and thus eliminates the above discussed problems. However, the somewhat reverse nature and the fact that the chemical potentials of the constituent elements are never explicitly determined may have kept it from being the standard approached used today. In addition, the formulation in terms of $\Delta\mu$ does not allow a straightforward extension to systems with more than two components. Here we will take this approach to a new level, where the desired stoichiometry is used as an input and where the explicit chemical potentials for all constituents can be determined.

As demonstration system, we choose cubic BAs, where we study the structure and energetics of its native point defects. BAs was initially predicted to be the best thermal conductor with an ultrahigh lattice room-temperature thermal conductivity of 2000 W/(mK) in 2013 lindsay2013first, and a few years later new theoretical work suggested a room temperature value of 1400 W/(mK) lindsay17. Meanwhile, a number of groups have succeeded in synthesizing boron arsenide crystals of increasingly better quality li2018high; tian2018unusual; kang2018experimental; tian2018seeded; kim2016thermal and high thermal conductivities around 1300 W/(mK) at room temperature have been found. Although by now experiment and theory converge on a value for the thermal conductivity limit in high-quality crystals, the question about the most prevalent native point defects in BAs is still unsettled. Several theoretical attempts have been made to unravel the nature of the dominant native point defects, using their connection to off-stoichiometry and their detrimental effect on the thermal conductivity measured in BAs lv2015experimental; protik2016ab; zheng2018antisite. However, there has been no convergence on the constitutional point defects to date, i.e. the point defects responsible for off-stoichiometric compositions.mishra2012native Initial work has focused on vacancies, which were chosen in an ad hoc fashion,lv2015experimental; protik2016ab and subsequently, Ref. zheng2018antisite suggested that BAs-AsB antisite pairs should be the lowest-energy defects. However, the stoichiometric composition of the pair, which in Fig. 2 from Ref. zheng2018antisite is the lowest-energy defect, violates the basic assumption of As-richness and indicates a flaw in the calculation: by adding the antisite pair as the only point-defect cluster out of the full ensemble of paired point defects, the antisite pair is not competing with any other cluster. Nevertheless, point defects very frequently have the tendency to cluster to minimize the number of mis-bonds, which gives the clusters lower energy than the sum of their constituent point defects and sometimes even lower than single point defects. The mismatch between the stoichiometry of the antisite pair and the stoichiometry premise thus can be taken as an indication that if a more complete set of clusters would have been considered, there should be another cluster that is As-rich which has a lower energy at the CNP. Taking cues from other tetrahedral compound semiconductors, point defect complexes confirmed by experiments are indeed often more complex than a simple antisite pair, such as the VCCSi-SiCCSi and SiC(C${}_{\rm Si})_{2}$ complexes identified by EPR as predominant annealing products in SiC in Ref. pinheiro2004silicon, which, translated to BAs, could be compatible with both the TEM images and the stoichiometry assumptions, while the predicted antisite pair from Ref. zheng2018antisite contradicts the underlying stoichiometry. This discrepancy suggests a methodological deficiency after the initial assumption of As-rich chemical potentials. Therefore, a theoretically sound answer to the question about the energetics of point defects in BAs is still outstanding, especially with a comprehensive examination of their connection to stoichiometry, which we will provide here. Since our goal is to have a method that shows better consistency between assumptions and results regarding chemical potentials and stoichiometry, we will restrict our examination to the case of single point defects, due to the potential misinterpretations resulting from comparing incomplete sets of defect complexes such as clusters to single point defects as we just discussed.

In this paper, we extend the method proposed in Ref. Northrup91 and reformulate a full canonical approach to explicitly calculate the chemical potentials self-consistently at both 0 K and finite temperatures. By using the stoichiometry as input variable and allowing the chemical potentials to be Fermi-level dependent, we treat both stoichiometric and off-stoichiometric conditions within the same framework. We also show that the off-stoichiometric limit saturates to “constitutional point defects” that have been defined within the grand-thermodynamic-potential approach developed within the field of intermetallics and alloys towards the end of the previous centuryMayer95; foiles_daw_1987; Schott97; mishra2012native; hagen1998point. We discuss how our approach can be used to determine physical chemical potentials, carrier concentrations and charge neutrality points in gapped materials.

We then present a systematic study of native point defects in cubic BAs, including vacancies, antisites and interstitials. Our results include the fully relaxed structures, their minimum-energy charge states, and their formation energies for the B-rich, As-rich and stoichiometric cases at 0 K and finite temperatures. For the latter, we calculate the fully temperature dependent free energy of formation within the quasi-harmonic approximation in the dilute concentration limit where the configurational entropy contribution is shown to be much smaller than that of vibrational entropy, and thus negligible in altering the equilibrium point defect concentrations. Our results at finite temperatures show that it is the boron vacancy and both antisites that play the most important role. In addition, we give a detailed discussion of CNP and thermodynamic equilibrium in BAs in light of these calculations, and determine quantitatively the $p$-type charge carrier concentration induced from self-doping in the frequently observed B-rich samples. These together validate our predictions as well as the underlying methodology.

II Methods

II.1 Theoretical framework

Based on the DFT total energies for defective and perfect cells as calculated in Sec. II.2, the defect formation energies can be calculated using Eq. 1. Moreover, we account for the vibrational contribution to the free energies of formation at finite temperatures, $F_{f}(T)$, which we obtain using the quasi-harmonic approximation baroni2001phonons based on supercell $\Gamma$-point calculations as described in previous work luo2009first and summarized in Sec. S2 of the Supplementary Material (SM). Free energies of formation of a certain point defect can be determined in analogy to Eq. 1 by replacing $E$ ($E_{f}$) with $F$ ($F_{f}$).

For determining the necessary chemical potentials from Eq. 1, we propose a canonical approach where they are determined by the stoichiometry of the crystal, starting from the method discussed in Ref. Northrup91. Here, the two constraints necessary to resolve $\mu_{\rm B}$ and $\mu_{\rm As}$ are the following. First, their sum is the total or free energy per formula unit of perfect BAs determined from ab initio calculations,

[TABLE]

where $E$(AB) is the total energy per AB formula unit, which is also utilized in the HoF approach and is commensurate with the limits defined in Ref. zhang1991chemical. However, unlike the HoF approach where either $\mu_{\rm B}$ and $\mu_{\rm As}$ is assumed equal to the elemental limit, the second constraint in our approach is that the overall stoichiometry is balanced by point defect concentrations, such that the ratio of the As and B number densities equals a predefined stoichiometry.

In a binary system with a given stoichiometry B1+δA1-δ, the concentrations of all defects are practically calculated from the free energy of formation ,Mayer95; mishra2012native

[TABLE]

where $C_{\mathrm{s}}$ is the concentration of lattice sites in the perfect material, $\theta_{X}$ is the multiplicity of the configuration (which is different from 1 e.g. for dumbbell interstitials which can have different orientations), and $k_{\mathrm{B}}$ is the Boltzmann constant. These allow us to calculate the concentrations of As and B atoms in BAs by

[TABLE]

where $V$ is the volume of the cell and $N$ is number of B or As atoms. For a given stoichiometry B1+δAs1-δ, the difference between the B and As concentrations needs to be $2\delta$, resulting in the equation

[TABLE]

For the total defect concentrations incorporating all charge states that is needed for Eq. 15, the concentrations in Eq. 14 need to be calculated for all charge states and summed up into a total defect concentration. For larger energy differences between the charge states and at sufficiently low temperature, that procedure can be simplified by only calculating the concentration for the lowest-energy point defect which will dominate.

Replacing the concentrations in Eq. 15 with the corresponding expressions from Eq. 3 with the free energies of formation following the form of Eq. 1 creates a second equation in addition to Eq. 2. The combination of Eqs. 2 and 15 allows solving for the chemical potentials. For systems with more than two components, Eq. 15 becomes a system of equations for pairwise stoichiometries. For the example of a ternary compound AmBnCl, one ends up with two equations for $\delta_{\rm AB}$ and $\delta_{\rm AC}$ (which already fix the BC ratio too), where the three chemical potentials can then be solved in combination with the counterpart of Eq. 2. The detailed equations for this example are given in Sec. S.3. in the SM. Since the free energies depend on temperature and Eq. 1 adds the Fermi level into the equation, the chemical potential $\mu(T,E_{F},\delta)$ as a function of $T$ and $E_{F}$ for the stoichiometric case ($\delta=0$) or a given non-stoichiometry $\delta$ can be determined.

In the limit of zero temperature where the Boltzmann factor becomes numerically inaccessible, the equation system collapses to the requirement that the formation energy of the lowest-energy stoichiometry balancing defects are equal.foiles_daw_1987 Specifically, one pair of defects resulting in B-rich and As-rich sample respectively should have the same formation energy, for example $E_{f}({\rm B_{As}})=E_{f}({\rm V_{B}})$ or $E_{f}({\rm B_{As}})=E_{f}({\rm As_{B}})$, or $E_{f}({\rm B_{As}})=E_{f}({\rm As_{i}})$, to name a few possible cases. The solution of the chemical potentials from these postulation must satisfy that all other defects are “irrelavant”, meaning that they have higher formation energies than the two stoichiometry-balancing defects. As we will show later, the 0 K chemical potential obtained in this way agrees with the value extrapolated from finite temperatures very well. One can already see that, at both 0 K and finite temperatures, the chemical potentials are obtained from the lattice stoichiometry in a self-consistent way, and all point defects should have non-negative formation energies, and thus physical equilibrium concentrations in our approach.

While Eq. 15 can be solved in principle for any off-stoichiometry $\delta$, as we show for the example of BAs in Sec. S.4 and Fig. S4 in the SM, the formation energy ($E_{f}$) of one of the point defects quickly saturates and approaches zero with increasing off-stoichiometry $\delta$, indicating that the chemical potential does not change upon alternating $\delta$ anymore. This then means that this zero-formation-energy defect is the one linked to the off-stoichiometry. This connection has been discussed before for point-defect calculations for intermetallic compounds in the 1980’s and 1990’s,Mayer95; foiles_daw_1987; Schott97; hagen1998point where the term “constitutional defect” was used for the zero-formation-energy point defect responsible for the off-stoichiometry, and the zero-formation energy was defined through the grand thermodynamic potential. Practically, this means that, for the example of a BAs crystal, one of the potential defects that can cause a B-rich sample, i.e. a Bi interstitial, BAs antisite, or VAs vacancy, will have zero formation energy, while all other defects will have higher (positive) formation energies. If for an assumed constitutional defect, another defect is found to have negative formation energy, that defect needs to be the constitutional defect instead. For the example of postulating that the BAs is the “constitutional defect” for B-rich BAs, the chemical potentials can be calculated by combining Eq. 2 with the requirement that $E_{f}({\rm B_{As}})=0$ and substituting for $E_{f}({\rm B_{As}})$ the equivalent expression from Eq. 1, $E({\rm B_{As}})-E_{\rm perf}-\mu_{\rm B}+\mu_{\rm As}+Q(E_{F}+E_{v}).$

Beyond the constitutional limit, our canonical approach is identical to the grand canonical approach, except that the chemical potential does not necessarily have to be the same with the elemental chemical potential. In fact, it should not be equal because the chemical potential in a compound is always lower than that in the elemental limit, otherwise the compound will never form. Still, comparing the constitutional chemical potentials with the elemental limits gives additional information about the stability of stoichiometric and off-stoichiometric compounds, as we discuss in the SM in Sec. S.3 and Fig. S3.

While the condition of charge neutrality is discussed in the review of the HoF approach from 2014Freysoldt2014first but typically not used in the reviewed body of work, it is a central point in the present work due to the Fermi-level dependence of our chemical potentials. While discussing the significance of charge neutrality in both approaches in detail in SM Sec. S.1, we start here from the typical condition that the CNP is determined in an undoped crystal by adding up the charges on the defects weighted by their concentration $C_{X}$, calculated by Eq. 3, for each possible Fermi level ($E_{F}$) in the band gap, and equating them to the defect-induced carrier population,

[TABLE]

where charges $Q_{X}$ on donor and acceptor point defects $X$ have opposite signs, such as $Q_{{\rm V}_{\mathrm{B}}^{+1}}=+1$ and $Q_{{\rm V}_{\mathrm{As}}^{-2}}=-2$. The resulting defect charge concentration, which we label $\Delta n_{D}(E)$, describes the electron or hole excess injected into the crystal from unbalanced charges on the defects. For extrinsic semiconductors, the extrinsic carrier concentration would have to be added into the equation, and electron vs. hole population would be subject to the mass action law.

The charge neutrality point in BAs is thus defined by the condition that the excess charge carrier density $n(E)-p(E)=0$ within the condition of stoichiometry balancing at a given temperature. For Fermi levels below the CNP, the equilibrium point defects add extra hole and above the CNP, extra electron concentration, which we express by $\Delta n_{D}$ (positive for excess electrons, negative for holes). Fermi level and defect-induced carrier concentration can then be determined by standard semiconductor physics by combining Eq. 6 with

[TABLE]

and solved self-consistently.

II.2 Computational details

The calculations were performed within density-functional theory (DFT) using the plane wave projector augmented-wave (PAW) method blochl1994projector as implemented in the VASP code kresse1993ab; kresse1994ab. A cutoff energy of 560 eV and a BZ sampling grid of 6$\times$6$\times$6 Monkhorst-Pack k-point monkhorst1976special in a 64-atom super cell were used in this work. Monopole correction $e^{2}Q^{2}\alpha/(a_{0}\varepsilon$) ($\alpha$ is the Madelung constant, a0 the lattice constant of the cubic cell and $\varepsilon$ the reported dielectric constant touat2006dynamical) is applied to calculate energies of charged systems windl1999first; Windl2004diffusion; Windl1998first. While structural relaxation was performed with the Perdew-Burke-Ernzerhof (PBE) functional perdew1996generalized, the Heyd-Scuseria-Ernzerhof (HSE) hybrid functional heyd2003hybrid was employed for bandstructure calculations to get the correct band gap and for the final total-energy calculations to maintain maximum accuracy. This procedure (PBE relaxation plus HSE final energy) is computationally efficient but compares favorably to full HSE-relaxation, as we tested for the example of BAs antisite, the constitutional defect for B-rich BAs, where the calculated “creation energy” Daw01, i.e. the difference between the defect cell and the perfect supercell (without compensating chemical potentials) results in 20.60 eV with our procedure vs. 20.55 eV for full HSE relaxation, a difference of 0.24%. The fraction of exact exchange in the HSE calculations was chosen as 0.16 to obtain an indirect zero-temperature band gap of 1.60 eV, since typical semiconductors have band gaps that shrink between 0 K and room temperature by approximately 0.1 eV varshni1967temperature and the measured room temperature gap for BAs is $\sim$1.50 eV chu1974preparation. The quasiharmonic calculations of the vibrational free energy contributions were performed as outlined in Ref. luo2009first as summarized in the SM, Section S.2.

III Results and Discussion

III.1 Intrinsic Fermi level and carrier concentration

Figure 1 shows band structure and density of states of perfect BAs.The effective mass of holes $m^{*}_{p}$ is obtained from an average of heavy hole and light hole and found to be 1.02 $m_{0}$. Similar to that in silicon, the effective mass of the anisotropic conduction band minima is characterized by a longitudinal mass ($m_{l}=1.11m_{0}$) along the corresponding equivalent (100) direction and two transverse masses ($m_{t}=0.513m_{0}$) in the plane perpendicular to the longitudinal direction.van2004principles With the sixfold degenerate conduction band minimum,van2004principles

[TABLE]

The intrinsic Fermi level (0.75 eV above the valence band edge) indicated in Fig. 1b has been calculated for a temperature of 1100 K, which is the temperature at which the BAs crystals of Lv *et al. *have been synthesized lv2015experimental, using

[TABLE]

where $E_{\mathrm{gap}}$ denotes the band gap between conduction band minimum and valence band maximum , $k_{\mathrm{B}}$ the Boltzmann constant and $T$ the temperature. With the same ingredients, the intrinsic carrier concentration can be determined,

[TABLE]

With $E_{g}=1.611$ eV, $m^{*}_{n}=2.19m_{0}$, $m^{*}_{p}=1.02m_{0}$, and

[TABLE]

we get

[TABLE]

$n_{i}$ has values of $6.56\times 10^{16}$ cm*-3* at 1100 K and $1.33\times 10^{6}$ cm*-3* at 300 K.

III.2 Structure and energetics of point defects

Figure 2 shows the most stable structures of the various defects in their minimum-energy charge state at mid gap as labeled. The perfect cell is relaxed to have a lattice constant of 4.821 Å compared to the experimental value of 4.7802 Å reported in Ref. lv2015experimental, with a B-As bond length of 2.088 Å. The arrows in Figs. 2a-d indicate changes in bond length as compared to perfect BAs, where we find that the lattice distortion from antisites is significantly larger than from vacancies. While B vacancies and antisites contract the lattice, their As counterparts expand it. For interstitials, the (100) dumbbell with identical species was found to be the most stable case for both boron and arsenic as compared to tetrahedral, hexagonal and (110) dumbbell interstitials. In Figs. 2e and f, the arrows represent the bond length between the two dumbbell atoms. The relaxations were performed without enforcing crystallographic symmetry to allow for symmetry breaking (e.g. Jahn-Teller distortions or off-center dumbbells), which however was found to be very small.

The formation energies of all defects as a function of Fermi level relative to the top of the valence band $E_{\rm VBM}$ in both stoichiometric and off-stoichiometric cases are shown in Fig. 3a-d. We find B-antisites to be the constitutional defects in the B-rich case for the entire Fermi-level range (Fig. 3a), while for the As-rich case, the As-antiste dominates up to $\sim 1.2$ eV, while for high Fermi levels, VB takes over (Fig. 3b). The stoichiometry-balancing defects are independent of temperature except for the close vicinity of the valence band edge at low temperatures. At 1100 K, which is about the annealing temperature in Ref. lv2015experimental, BAs and VB balance each other to maintain stoichiometry over the entire Fermi-level range (Fig. 3d), while the zero-temperature calculations predict stoichiometry balancing between AsB and BAs instead of VB in the severe $p$-type region (Fig. 3c).

III.3 Charge neutrality

In order to explore the relation between charge neutrality, Fermi level and non-stoichiometry in our results, we start from a short discussion of the charge neutrality point (CNP). The term CNP is historically used to indicate the “pinning” of the Fermi level at interfaces of a semiconductor with other materials Spicer79 and determines the change of the Fermi level with increasing defect concentrations.Brudnyi95 Thus, it is the essential link between point defect concentrations and measurable electronic properties and is explained in more detail in SI Sec. S.1.

Within the HoF approach, where the chemical potentials are Fermi-level independent, the charge of a defect can be determined by its slope in the plot of formation energy vs. Fermi level. This allows to identify the zero-temperature CNP from the lowest-energy intersection between two formation energy curves with slopes of the same magnitude but opposite signs. This is important, since the formation energies at the CNP are the in principle only sensible ones as opposed to other Fermi level positions, where the formation energies can even be negative (see Sec. S.1). This method would not work within our framework, where the chemical potentials are Fermi-level dependent, while the favorable charge state of a certain defect is the one leading to the lowest formation energy. Thus, the slope of the formation energy curves is not related to the defect charge. For example, the BAs antisite is neutral at $E-E_{\rm VBM}<0.265$ eV (Fig. 3e), but the formation energy does not have a zero slope at the same range (Fig. 3c). Therefore, the charges need to be separately recorded and reported. Figures 3e-f show the charge states for the native point defects at both 0 K and 1100 K, where the ionization levels $\varepsilon(q_{1}/q_{2})$ for the different defects can be identified from the steps in the charge vs. Fermi level plots.

In the spirit of stoichiometry balancing, the CNP is determined quantitatively by Eq. 6 using the defect concentrations as defined in Eq. 3. Figure 4a shows the electrons or holes injected from all possible native point defects in stoichiometric BAs at 1100 K, unbalanced by the defect charges. The excess charge density $|n-p|$ will diverge at the CNP, indicated by the grey line in Fig. 4a. For stoichiometric BAs, as calculated at 1100 K in Sec. III.1, the intrinsic Fermi level (Eq. 9) is 0.75 eV and the intrinsic carrier concentration (Eq. III.1) is $6.56\times 10^{16}$ cm*-3*, while the excess charge concentration is $|n-p|\sim 6.9\times 10^{9}{\rm cm}^{-3}$, seven orders of magnitude lower. With Eq. 7, the Fermi level stays almost in the middle of the gap as it should, since the overall Fermi energy in a bulk semiconductor is a thermodynamic average and the excess carriers from the charged defects could barely shift it significantly. This situation can be very different in non-stoichiometric BAs, where the concentration of the constitutional defects can be in the percent range, potentially resulting in a significant shift of the Fermi level. The connection between stoichiometry, CNP, and resulting Fermi level is thus a straightforward outcome of our approach, while in traditional HoF work it is usually not possible due to the undefined relationship between stoichiometry and chemical potentials. In the following, we will use the CNP in the traditional sense, where it refers to the charge-balancing Fermi level only from the native point defects.Spicer79

Next, we compare the CNP from stoichiometry balancing to that determined from HoF plots. Due to the instability of defect concentrations at extremely low temperatures, we extrapolate the CNP from a series of finite-temperature results. With that, we find it to be $0.40$ eV above the valence band edge at 0 K (Fig. 4b). Previous work at 0 K has also found it to be in the $p$-type regime below mid gap at $\sim 0.6$ eV above the valence band edge.zheng2018antisite Possible reasons for the slight discrepancy include that Ref. zheng2018antisite assumes $\mu_{\rm As}$ to be at its elemental limit, the inclusion of zero-point energies in our results which influence the defect concentration balancing, and the different HSE mixing parameters in the HSE functional, 0.25 in Ref. zheng2018antisite vs. 0.16 here (see Sec. II.2).

Further, we compare the intrinsic-defect-doping in the different stoichiometry limits with the intrinsic carrier concentration and extrinsic-impurity-doping. The atomic concentration of B/As in a stoichiometric sample is $\sim 3.67\times 10^{22}$ cm*-3*, where Eq. III.1 calculates the intrinsic carrier concentration to be $6.56\times 10^{16}{\rm cm}^{-3}$ at 1100 K and the thermal point defects are predicted to produce $n$-type doping with an electron concentration of $6.9\times 10^{9}$ cm*-3* at 1100 K, seven orders of magnitude smaller than the intrinsic carrier concentration, which then leaves the Fermi level at almost mid gap as expected for a stoichiometric semiconductor. However, for the heavily B-rich sample from Ref. lv2015experimental, assuming all off-stoichiometry is caused by the B${}_{\rm As}^{-1}$ antisite as constitutional defect, the 2.8% As-deficiency corresponds to a B${}_{\rm As}^{-1}$ concentration of approximately $5.14\times 10^{20}$ cm*-3* acting as acceptors. In the non-degenerate region, the carrier concentration can be estimated by $N_{A}/\{1+4{\rm exp}[(E_{A}-E_{v})/k_{B}T]\}$ for holes, where $N_{A}=5.14\times 10^{20}$ cm*-3* is now the acceptor concentration and $E_{A}=0.265$ eV is the ionization energy (Fig. 3e). Thus the hole concentration from self-doping is $\sim 4.45\times 10^{15}$ cm*-3* at 300 K, and $\sim 7.68\times 10^{18}$ cm*-3* at 1100 K. In combination with the low enthalpy of formation of B12As2, which puts the BAs phase above the convex hull and makes it much less likely to grow an As-rich sample (SM Sec. S.5), this leads to the conclusion that standard growth techniques most likely produce p-type boron rich samples. Although the hole concentration from the B${}_{\rm As}^{-1}$ antisites alone are not as high as the magnitude of $10^{19}$ cm*-3* from earlier works which suggests that Si/C impurities are indeed dominant in the $p$-type conduction, chae2018point; lyons2018impurity; meng2020pressure our results indicate that even without impurities, the as-grown samples would still be $p$-type, although to a lesser degree. Our canonical native point defect formation energies thus enable a quantitative analysis on the origin of holes from the native vs. extrinsic defects, thus allowing the assessment of the effect of sample purity vs. sample stoichiometry. To this end, our canonical approach opens up new opportunities to quantitatively assess the self-doping in non-stoichiometric samples, especially important in materials where intrinsic point defects can have a smaller ionization energy.

III.4 Free energies and thermodynamic stability

Since we have calculated the full free energies of formation, we now plot the various defect concentrations as a function of (inverse) temperature in an Arrhenius plot as shown in Fig. 5 for intrinsic BAs at mid gap in order to examine if zero-temperature formation energies are a viable way to estimate defect concentrations and thus can reasonably well substitute for the formation enthalpies $H_{f}$ from the full free energies of formation (independently of how the chemical potentials are determined). We choose to examine them at mid gap because as discussed earlier, the intrinsic carriers dominate in the chosen temperature range and leave the overall Fermi level close to mid gap, with only small shifts from the charged defects as long as we assume small non-stoichiometries. All lines are in very good approximation linear, and with Ff = Hf $-$ TSf, we can use a linear fit to determine formation enthalpy and entropy for each defect, shown in Table 1. These formation enthalpies contain contributions from zero-point energies and compensate for non-linearities from the vibrational entropy contribution, and thus cannot a priori be expected to be identical to the formation energies calculated from zero-temperature total energies. The number of digits shown in the enthalpy column indicates the degree of linearity of the fit, and the error bars shown limit the 95% confidence interval. The concentrations follow an Arrhenius temperature dependence to an excellent degree, although that does not necessarily have to be the case for the type of calculations performed here as previously found stoddard2005ab; grabowski2009ab; Grabowski2010formation; glensk2014breakdown. We approximate the configurational entropy per atom in the dilute limit by Freysoldt2014first

[TABLE]

where $c=n/N$ with $n$ being the number of point defects of a specific type and $N$ is the number of lattice sites; $\theta_{X}$ is again the configurational degree of freedom. For the most relevant defects, i.e. antisites and vacancies, $\theta_{X}=1$. Even if we use the simulated supercell to account for defect concentration (1 defect in 32 lattice sites thus $c\sim 3\%$), which is usually exaggerated, the configurational entropy is about $S_{f}^{\rm conf}=0.1396k_{\rm B}$, corresponding to $-0.012$ eV at 1000 K. This contribution is negligible compared to the vibrational entropy listed in Table 1.

While the agreement between calculated formation energies and fitted formation enthalpies is in general good with an average deviation of 5%, the considerable variations in the (vibrational) formation entropies shown in Table 1 cause some larger deviations up to 14%. While the sequence and relative importance of nearly all point defects is not dependent on temperature, there is one exception: the AsB antisite dominates close to zero temperature in proximity to the valence band edge over VB, while VB dominates otherwise in most of the relevant temperature and Fermi level range as the stoichiometry balancing defect. In combination with the significant variations in formation entropy in Table 1 and thus different temperature dependence, our results suggest that a calculation of free energies is indeed an important step to determine the dominant defects at finite and especially high annealing temperatures, while the calculation of zero-temperature formation energies give a reasonable estimate for the general trends and for identifying the most relevant defect types. While this was done for the specific choice of chemical potentials proposed in our paper (stoichiometric chemical potential from the canonical assumption for $N_{\rm B}=N_{\rm As}$; B- and As-rich cases from the assumption that the respective constitutional defects have zero formation energies), this comparison could have been done the same way for HoF-determined chemical potentials.

The combination of Figs. 3a and 5a shows that for B-rich material, the highest concentration of thermal, non-constitutional defects is below 108 cm*-3*, and thus all thermal defects should be negligible when calculating materials properties. This means the only relevant defect in B-rich material, which is the most relevant case for comparison to experiment as we just showed, is the boron antisite. The situation is different for As-rich BAs (Figs. 3b and 5b), where boron vacancies can reach concentrations in the 1020 cm*-3* range, which may require including them into transport and carrier-concentration calculations in addition to the AsB antisites, while other defects still have negligible concentrations. For stoichiometric BAs (Figs. 3c-d and  5c), we find that VB and BAs balance the stoichiometry, with the equilibrium concentration of B vacancies twice that of BAs antisites, since one BAs needs two B vacancies to result in perfect stoichiometry. VB and BAs concentrations have relatively low concentrations of $\sim$1010 cm*-3* at elevated temperatures and significantly less around room temperature, giving them negligible relevance for electronic properties and thermal conduction in stoichiometric BAs and in agreement with the recent findings of close-to-ideal thermal conductivity in apparently stoichiometric BAs.kang2018experimental

At this point, it is important to note that the chemical potentials numerically solved in a self-consistent manner as a function of stoichiometry, temperature and Fermi level $\mu(E_{F},T,\delta)$ actually fulfill the stability criterion and do not exceed the boundaries set by the elemental limits (for more details see SM Sec. S.3). It is well known that increasing the temperature will facilitate defect formation, which in our temperature-dependent approach manifests itself in the fact that the temperature-dependent elemental limits for the chemical potentials form a wider range at elevated temperatures, where the numerical solution for the chemical potential stays well within the elemental limits for the entire Fermi level range (Fig. S3 in SI). Importantly, the expansion of the chemical potential range at higher temperatures agrees with the phase stability analysis regarding the BAs and the B12As6 phases on a convex hull (Fig. S5 in SI). Therefore, this method with stoichiometry-linked chemical potentials can also be useful for tackling for example solubility ranges and phase stability problems in alloys or intermetallic compounds.

IV Conclusions

In summary, we have presented a new approach to calculate defect formation energies in non-metals based on stoichiometry-linked chemical potentials and shown results for the energetics of native point defects in cubic boron arsenide. We have demonstrated how, under the conditions of thermodynamic equilibrium and overall stoichiometry, chemical potentials are functions of temperature, Fermi level, and stoichiometry, and have demonstrated how they can be determined for our demonstration system BAs. With those, we have calculated the formation energies and equilibrium concentrations of the native point defects, which included vacancies, antisites and interstitials.

We find that antisites are the constitutional defects in B-rich and to a large degree in As-rich material, while BAs antisites and B vacancies balance to keep stoichiometry over most conditions, while AsB replaces VB in extremely $p$-type material at low temperatures. In addition, full free energies of formation were calculated for a range of temperatures, and we find that considerable differences up to 14% exist between Arrhenius fitted formation energies and 0 K total-energy values.

This work takes an important step further to practically include charge neutrality in a gapped material from the charged point defects in the analysis, which to a very large degree has been left ambiguous in most of the HoF literature. We have shown the critical significance of sensible charge neutrality analysis for our demonstration system, where we find that the as-grown B-rich samples are already intrinsically $p$-type, independent of the external doping through Si or C found to dominate doping in previous work.chae2018point; lyons2018impurity; meng2020pressure

While the widely used HoF approach seems to be at the first glance transparent and straightforward, some care is required when interpreting the results. If that care is taken, the HoF is a quick way to get an idea about dominant point defects and its general influence on the Fermi level, albeit not in a quantitative way. Our stoichiometry-balancing approach allows to correlate chemical potentials with stoichiometry in a physically rigorous way for a real material, and therefore allows a quantitative insight into point defect populations instead of the qualitative picture that is gained from common HoF calculations, which in turn enables a quantitative connection to materials properties in semiconductors. Finally, we would like to point out that this method, when collapsed to a zero-gap system, directly connects to previous approaches for point-defect calculations in alloys and intermetallic compounds,Mayer95; foiles_daw_1987; Schott97; mishra2012native; hagen1998point where it allows a quantitative evaluation of solubilities and point defects, especially the constitutional ones.

V Data Availability

The data are available from the corresponding author on request.

VI acknowledgments

This work was funded by the Air Force Office of Scientific Research under award number FA9550-14-1-0322. Y. W. also acknowledges partial funding from the Center for Emergent Materials, an NSF MRSEC (DMR 1420451). All high-performance computations were performed on machines of the Ohio Supercomputer Center, grant number PAS0072.

VII Supplementary material

See supplemental material for detailed discussions on charge neutrality point, quasiharmonic approximation, chemical potential saturation, as well as thermodynamic equilibrium and phase stability of BAs.

Table of Contents

Section S.1 Charge Neutrality Point

Figure S1. Fermi-level dependent formation energies of point defects in Al2O3 and BAs from previous work.

Section S.2 Quasiharmonic Approximation

Figure S2. Convergence test for free energies in BAs.

Section S.3 Thermodynamic Equilibrium of BAs with Equilibrium Point Defects.

Figure S3. Chemical potentials for As in stoichiometric BAs at 0 K and finite temperature.

Section S.4 Transition from canonical to grand canonical at off-stoichiometric limit

Figure S4. Stoichiometry (As/B atomic ratio) dependent defect formation energies in BAs.

Section S.5 Phase Stability of Cubic Stoichiometric BAs

Figure S5. Convex hull plot of the formation energy of the phases in the B-As system.

VIII Charge Neutrality Point From Heat-of-formation Method

Upon irradiation, the Fermi level of semiconductors starts to shift due to the point defects created, until it eventually stabilizes at a “limiting” value or neutrality level, from where it does not change with further irradiation Brudnyi95. Also, the CNP can pin the Fermi level at interfaces, which can lead to unexpected band line-ups Spicer79. Thus, the charge neutrality point (CNP) is a crucial ingredient for the interpretation of point defect energies in binary compounds. As an example for the meaning of the CNP in point-defect calculations, Fig. 6(a) shows the Fermi-level dependent formation energy of native point defects in Al2O3, redrawn after the results from Ref. Choi13. The formation energies were calculated within the heat-of-formation approach for the Al-rich case, i.e. the chemical potential of Al, $\mu_{\rm Al}$, was set equal to the atomic energy in elemental aluminum, and $\mu_{\rm O}=\left[E_{\rm tot}({\rm Al}_{2}{\rm O}_{3})-2\mu_{\rm Al}\right]/3$.

We have added into the figure a black dashed line at zero formation energy. As can be seen, a lot of the calculated formation energies are in part significantly below zero, which defies the definition of formation energy in a stable crystal and thus is unphysical. While it is usually not discussed explicitly, there is actually only one Fermi-level position in this graph where the results within this approach are sensible, which is at the Fermi level value that defines the CNP, indicated by the red dashed line in the figure. At that Fermi level value, the line for the lowest-energy point defect with positive slope intersects the lowest-energy negative-slope line with the same absolute value of slope. Since the slope in this fixed-chemical-potential approach corresponds to the point defect?s charge, at this point, positively and negatively charged point defects have the same formation energy and thus the same concentration, compensate each other?s charge and thus are used to define what Fermi level is sensible for the chosen chemical potential. Here, Ali (slope $+3$) and VAl (slope $-3$) intersect at 4.9 eV, which thus defines the CNP. All the regions with the unphysical negative formation energies are thus irrelevant. At the CNP, the stable point defect (lowest energy) is V${}_{\rm O}^{0}$ with zero charge. If thus V${}_{\rm O}^{0}$ is the predominant point defect, the crystal will be Al-rich without losing charge neutrality, which is the premise of these results, and thus we end up finally with a piece of physical and valuable information and have identified the constitutional point defect and its charge. One could pick a different set of chemical potentials between the elemental limits and end up with a different CNP Fermi level. More specifically, in Fig. 6(b) for the O-rich case, VO (slope +2) intersects with VAs (slope -2) at around 2 eV, where the defect having the lowest formation energy is the uncharged oxygen interstitial, corresponding to an O-rich crystal. In a similar fashion, as long as only single point defects are considered, Ref. zheng2018antisite finds for the As-rich case of BAs, that at the charge neutrality point, which is 0.5-0.6 eV above the valence band edge, the neutral boron vacancy has the lowest energy and thus is constitutional defect, resulting in an As-rich crystal in agreement with the premise (Fig. 6(c)). However, the ad hoc addition of the antisite pair to the complete set of single point defects now suddenly seems to suggest that the antisite pair would be the constitutional defect. Since however the antisite pair preserves the perfect stoichiometry, a contradiction is now created to the assumption of an As-rich crystal, suggesting a methodological problem by singling out one ad hoc cluster to throw into the mix. Independent of that, all three examples show that chemical potentials are clearly coupled to the Fermi level once a stoichiometry is chosen, also in the common heat-of-formation approach. While the CNP determined in this way can explain qualitatively the effect of Fermi-level pinning described above, it does not allow for a quantitative connection between stoichiometry and Fermi level, and does not connect to the well-known fact that a crystal with very few point defects has its Fermi level at mid gap.

IX Quasiharmonic Approximation

As described in Ref. luo2009first, we used the quasiharmonic approximation (QHA) with VASP input for energies and phonon frequencies to calculate the vibrational contribution to the free energy, which approximates the free energy by

[TABLE]

where $\omega$ denotes phonon frequencies, $g$ is the phonon density of states (PDOS), the second term on the right hand side is the zero-point energy, and the third term is the temperature-dependent contribution. In order to calculate the PDOS, we determine the normal modes in the defect supercell and equivalent perfect supercell using density-functional perturbation theory within VASP with the same settings as for the PBE calculations. The PDOS is calculated from a sum over the calculated $\Gamma$-point phonons, where each of them is broadened by a Gaussian of width 0.5 THz. As we have shown in Ref. luo2009first, this rather fast method agrees typically within a few percent with a free-energy calculation from the PDOS calculated from a full integration by tetrahedron method over the Brillouin zone. This is expected since the $\Gamma$-point is a “special” k-point as defined first by Chadi and Cohen Chadi73 that is also frequently used to approximate electronic Brillouin-zone integration, and the accuracy of such single or multiple special points has been discussed in the childhood of DFT calculations extensively by e.g. Baldereschi baldereschi73 or Monkhorst and Pack monkhorst1976special, respectively. Since all structures are relaxed to zero pressure, the enthalpy term $H(p)$ is equal to the total energy. As additional approximation for computational efficiency, the effect of thermal expansion on the phonon frequencies is neglected.

While in the calculations in the paper the actual numerical results for the specific temperatures are used, for reference, we report a fit of the vibrational contribution here based on a fourth-order polynomial without linear term, which fits the actual DFT results over the temperature range within 0.01 eV per atom. As reference and for cross-validation, we report here results for solid As (relaxed minimum-energy $R\bar{3}m$ structure from MaterialsProject.org within a $6\times 6\times 2$ 432-atom supercell), B (relaxed minimum-energy $R\bar{3}m$ structure from MaterialsProject.org within a $3\times 3\times 1$ 324-atom supercell), we find

[TABLE]

and

[TABLE]

respectively. For BAs in zinc blende structure, we perform a convergence test with respect to supercell size. The results from a $4\times 4\times 4$ 512-atom supercell,

[TABLE]

are nearly identical to those from a $2\times 2\times 2$ 64-atom supercell,

[TABLE]

as also can be seen in Fig. 7.

X Thermodynamic Equilibrium of BAs with Point Defects

In the main text, for finite temperatures, the chemical potentials of As and B are determined from the requirement that thermal point defects maintain stoichiometry, i.e. that the system is closed and no atoms are exchanged with the environment when point defects are created. Within this approach, point defects are then coupled to stoichiometry, since they are the only way non-stoichiometry in a crystal such as B1+δAs1-δ can be accommodated for $\delta\neq 0$. The point defects that create the specific stoichiometry with the lowest energy penalty are called the “constitutional defects.”

The concept of constitutional point defects and coupling point defects to the specific stoichiometry is based on the grand thermodynamic potential and is in use since about three decades, where first applications targeted intermetallics Mayer95; Schott97; hagen1998point with later extension to insulating systems mishra2012native. The ground state of an ordered binary compound that is off-stoichiometric, such as B0.52As0.48, is defined by the system with the lowest energy. Since it is off-stoichiometric, the “perfect” off-stoichiometric system needs to have constitutional point defects to accommodate the off-stoichiometry, for example antisites or vacancies. Assuming they don?t interact and thus no second phase precipitates, their formation energy with respect to the “perfect” non-stoichiometric system (which they define) has to be zero, otherwise they would not be the constitutional point defects. For perfectly stoichiometric crystals, there are no constitutional point defects, but only thermally activated ones according to their formation energies, whose concentration is typically many orders of magnitude lower than that of the constitutional point defects, since measurable off-stoichiometry is typically only created by defects with formation energies of at most a few tenths of an eV.

In practice, the concentrations of all defects $X$ are expressed in the form of Eq. (2) in the main manuscript, from which the concentrations of As and B atoms in BAs can be calculated by mishra2012native

[TABLE]

For B1+δAs1-δ, the difference between the B and As concentrations needs to be $2\delta$, resulting in the equation

[TABLE]

Replacing the concentrations in Eq. (15) with the corresponding expressions from Eq. (2) with the free energies of formation following the form of Eq. (1) creates a second equation in addition to $\mu_{\rm B}+\mu_{\rm As}=F({\rm BAs})$ as described in the main text that allows solving for the chemical potentials. Since the free energies depend on temperature and Eq. (1) adds the Fermi level into the equation, the chemical potential $\mu(T,E_{F},\delta)$ as a function of $T$ and $E_{F}$ for stoichiometry ($\delta=0$) or a given non-stoichiometry $\delta$ can be determined. For a ternary system AmBnCl, there will be four conditions following the notation in the main text

[TABLE]

where only three in four equations are linearly independent, allowing us to solve the chemical potentials for A, B, and C respectively.

Next, we show that the chemical potentials generated this way fulfill the stability criterion and do not exceed the boundaries set by the B-rich and As-rich limits, i.e. the elemental chemical potentials of B and As. Figure 8 shows the As chemical potential as a function of Fermi level position as determined from the requirement of stoichiometry in comparison to the limits of As-rich and B-rich material, all calculated within QHA as described in Sec. IX. The figure shows that for elevated temperatures, the chemical potential is within the limits, indicating that the solution found from the stoichiometry requirement indeed fulfills the condition of thermodynamic stability for all Fermi level positions. At 0 K, the allowed range for the chemical potential of As ($\mu_{\mathrm{As}}$) is much narrower than that at higher temperatures, and for extreme-$p$ type conditions with Fermi levels $<0.25$ eV, stoichiometric BAs is unstable. More restrictive than that, Zheng et al. zheng2018antisite even argued that stoichiometric BAs at zero temperature is not a thermodynamically stable phase at all. To shine more light on this question, we examine the thermodynamic stability of stoichiometric BAs in Sec. XII. We further note that, elemental As has a sublimation point of $\sim 887$ K. In fact, the upper limit of $\mu_{\rm As}$ should be the minimum value within all possible phases. Here we take the limit of $\mu_{\rm As}$ in its solid phase, calculated by the quasiharmonic approximation, because the chemical potential in the gas phase varies significantly upon the partial pressure and so on.

XI transition from canonical to grand canonical at Off-stoichiometric limit

As discussed in Sec. IIA Theoretical framework in the main text, our canonical approach allows the chemical potentials, thus the defect formation energies, in a compound to be determined from a specific lattice stoichiometry B1+δAs1-δ. Here we show in Fig. 9 that, in BAs, the self-consistent solution for the chemical potentials do not vary upon the atomic ratio of As/B and the formation energies of the lowest energy defect approaches zero when the off-stoichiometry $\delta$ is on the order of $\pm 10^{-6}$. These lowest formation energy defects are called “constitutional defects”. Therefore, for an experimentally detectable composition, the off-stoichiometry limit thus aligns with the grand canonical approach practically. However, instead of being equal to the elemental limit, the chemical potentials are determined by the condition that the constitutional defect has zero formation energy while others have positive values.

XII Phase Stability of Cubic Stoichiometric BAs

In this section we discuss the phase stability of BAs with respect to its left and right neighbors in the phase diagram, elemental B and B6As, by constructing the convex hull at different temperatures from free energies of formation. To consider the effect of thermal expansion, the temperature-dependent lattice constants of the four phases, i.e. elemental B tsagareishvili1986thermal, elemental As fischer1978debye, pure BAs Chen2019, and pure B6As Whiteley2013 are taken from experiments to include thermal expansion as lowest-order anharmonic effect into the QHA framework as described in Sec. IX.

Figure 10 shows the convex hull plot Hautier2014 for the formation energies of the BAs and B6As phases. The entropy contribution increases the stability of the BAs phase while decreasing that of the B6As phase. By calculating the energy difference between the free energies of formation of BAs and the convex hull, the beneficial effect of elevated temperature can be clearly shown (Fig. 10b). Only energies up to 700 K are shown because As has a sublimation point $\sim 887$ K. Extrapolating from the values up to 700 K with a second-order polynomial fit estimates that the BAs phase will be located on the convex hull at temperatures around 1900 K. This estimate does not consider the excess As:B ratio used in synthesis, which may modify the thermodynamic equilibrium to lower temperatures and enable growth of stoichiometric BAs. This is consistent with the fact that at elevated temperatures, the chemical potential range allowed for stoichiometric BAs is much wider than that at 0 K, as shown in Fig. 8.