First-principles study of defect formation energies in LaO$X$S₂ ($X=$ Sb, Bi)

By $N_{1}\times N_{2}\times N_{3}$ supercell calculation, the formation energy of a point defect $D$ in charge state $q$ is evaluated as defect1 ; defect2

where $\mathbf{N}=(N_{1},N_{2},N_{3})$, $E[D^{q};\mathbf{N}]$, $E_{\mathrm{corr}}[D^{q};\mathbf{N}]$ and $E_{\mathrm{P}}$ represent the supercell size, the total energies of the supercell with the defect, its energy correction described below, and the total energy of the perfect unit cell without any defect, respectively. $n_{i}$ represents the number of removed (with a sign of $-$) or added (with a sign of $+$) atoms $i$, the chemical potential of which is denoted as $\mu_{i}$. $\epsilon_{\mathrm{VBM}}$ is the energy level of the valence-band maximum (VBM), and $\epsilon_{\mathrm{VBM}}+\Delta\epsilon_{\mathrm{F}}$ represents the Fermi level of the system.

In this study, we considered the energy correction $E_{\mathrm{corr}}[D^{q};\mathbf{N}]$ consisting of the following two terms:

The first term $E_{\mathrm{pc}}(q;\mathbf{N})$ denotes the point-charge correction Makov_Payne ; defect1 , which was evaluated through the Ewald summation for screened Coulomb interaction in a periodic supercell consisting of a single point charge $q$ in a uniform background charge $-q$. The screened Coulomb interaction is represented using the macroscopic static dielectric tensor in the way described in Ref. [diele_defect, ]. As described later in this section, the remaining finite-size error is removed by the extrapolation with respect to the supercell size.

The second term $E_{\mathrm{be}}[D^{q}]$ in Eq. (2) represents the band-edge correction, which is necessary when one uses a calculation method with a sizable band-gap error, such as the local density approximation and the generalized gradient approximation in density functional theory. We applied the band-edge correction for a shallow defect level as described in Refs. [Zunger1, ; Zunger2, ]. Suppose a point defect $D$ with a charge $q_{0}$ does not provide any carrier to the system (e.g., $q_{0}=1$ for $\mathrm{F}_{\mathrm{O}}$ that represents F^- substitution for O^2-), and then we define

where the band-edge correction to the valence band maximum $\Delta_{\mathrm{bg}}^{\mathrm{v}}$ and the conduction band minimum $\Delta_{\mathrm{bg}}^{\mathrm{c}}$ are evaluated by elaborate approximations in first-principles calculation that can provide an accurate band gap, as described in the next section. A relative energy level from the band edge, where the formation energy curves of $D^{q}$ and $D^{q_{0}}$ intersect, is kept unchanged by this correction. However, this approximation is not necessarily valid when the impurity levels are not shallow. Thus, we did not consider $E_{\mathrm{be}}[D^{q}]$ for that case. In the following sections, we mention whether $E_{\mathrm{be}}[D^{q}]$ was applied for each defect. In any case, the VBM energy in Eq. (1) is defined including this correction:

where $\epsilon_{\mathrm{VBM}}^{0}$ is an uncorrected VBM energy. This modification of $\epsilon_{\mathrm{VBM}}$ changes the origin of $\Delta\epsilon_{\mathrm{F}}$ to the corrected VBM energy. We consider defect formation energies for $\Delta\epsilon_{\mathrm{F}}$ lying between the corrected band edges.

After calculating $E[D^{q};\mathbf{N}]+E_{\mathrm{corr}}[D^{q};\mathbf{N}]$ in Eq. (1) for a set of supercells with several $\mathbf{N}$, we performed the least-squares fitting of $E[D^{q};\mathbf{N}]+E_{\mathrm{corr}}[D^{q};\mathbf{N}]$ by $c_{0}+c_{1}N^{-1}$, where $c_{0}$ and $c_{1}$ are coefficients of this fitting, and $N=N_{1}N_{2}N_{3}$ (cf. Ref. [defect1, ] for the $N^{-1}$ (or equivalently, $\Omega^{-1}$ with $\Omega$ being the supercell volume) dependence). After determining $c_{0}$ and $c_{1}$, we took the $N\to\infty$ limit, i.e., simply took $c_{0}$, to get the thermodynamic limit. By this procedure, we can get the defect formation energy at the thermodynamic limit,

which shall be shown in the following sections.

We used the Perdew–Burke–Ernzerhof parametrization of the generalized gradient approximation (PBE-GGA) PBE and the projector augmented wave (PAW) method paw as implemented in the Vienna Ab initio Simulation Package vasp1 ; vasp2 ; vasp3 ; vasp4 . The core electrons in the PAW potentials were taken as follows: [He] for O and F, [Ne] for S, [Kr]$4d^{10}$ for La and Sb, and [Xe]$5d^{10}4f^{14}$ for Bi. The spin-orbit coupling (SOC) was not included unless noted because a huge computational cost is required for supercell calculations including SOC. We performed spin-unpolarized calculation unless noted, also due to high computational cost. For ions with a closed-shell electronic configuration such as F^-, the spin polarization is expected to be unstable. Thus, our approximation likely works well for point defects involving such an ionic state. On the other hand, this is not always the case for some point defects, and we note that our approximation can cause some error in that case. The plane-wave cutoff energy of 500 eV and the Gaussian smearing with the smearing width of 0.15 eV were used. The structural optimization was performed until the Hellmann–Feynman force becomes less than 0.02 eVÅ^-1 on each atom. All the calculations were performed at zero temperature and zero pressure.

For calculating O₂ and F₂ molecules, we placed an isolated molecule in a $15\ \mathrm{\AA}\times 15\ \mathrm{\AA}\times 15\ \mathrm{\AA}$ cell, and only the atomic coordinates were optimized. Spin-triplet oxygen molecule was calculated using the spin-polarized calculation. For calculating several solid compounds used for evaluating the chemical potential $\mu_{i}$ of each atom $i$, their crystal structures were optimized by the computational conditions described in the previous paragraph. A sufficiently fine $\bm{k}$-mesh was taken for each material. The lattice parameters and the atomic coordinates were optimized within a constraint of a fixed space group.

For supercell calculations of LaO$X$S₂ ($X=$ Sb, Bi), we first optimized both the atomic coordinates and the lattice parameters of the perfect crystal (i.e., without any defect), assuming the space group $P$2₁/$m$ (monoclinic). This space group was realized in LaOBiS₂ SymLowPress ; SymLow ; BiS2_HP . Figure 1 shows the crystal structure of LaOBiS₂. Also for Sb compounds, this space group was experimentally reported for Ce(O,F)SbS₂ CeOFSbS2 and Ce(O,F)Sb(S,Se)₂ CeOFSbSSe2 , and was theoretically suggested to be stable for LaO$X$S₂ ($X=$ As, Sb, Bi) calc_P21m and LaOSbSe₂ Hirayama . The space group is not necessarily common in $Ln$O$XCh_{2}$ compounds, because several BiS₂ compounds, such as CeOBiS₂ SymLow , and also $Ln$OSbSe₂ ($Ln=$ La, Ce) LnOSbSe2 were experimentally reported to have the tetragonal $P$4/$nmm$ space group. Nevertheless, in this study, we can assume $P$2₁/$m$ in the calculation without loss of generality because this is a subgroup of $P$4/$nmm$, and found that the crystal structure indeed exhibits the monoclinic distortion as shown in Table 1.

After calculating the perfect crystal, we performed the defect calculation by putting an isolated defect in the supercell with different sizes. Here, we considered $3\times 3\times 2$, $4\times 4\times 2$, and $5\times 5\times 2$ supercells, which contain 180, 320, and 500 atoms, respectively. For each supercell, we used $4\times 4\times 1$, $3\times 3\times 1$, and $3\times 3\times 1$ $\bm{k}$-meshes, respectively. Supercell calculation optimized the atomic coordinates using the fixed lattice parameters optimized for the perfect crystal. This treatment is justified because we are interested in the dilute limit of a defect.

For calculating the density of states (DOS), we used $24\times 24\times 4$ and $5\times 5\times 1$ $\bm{k}$-meshes for the primitive cell and the $4\times 4\times 2$ supercell, respectively. We performed these DOS calculations using the fixed charge density obtained in the self-consistent-field (SCF) calculation.

We calculated the macroscopic static dielectric tensor based on the density functional perturbation theory to evaluate the point-charge correction. We considered both the electronic and the ionic contributions of the dielectric tensor. The local field effect for the Hartree and the exchange-correlation potentials was included. The derivative of the cell-periodic part of the Kohn–Sham orbitals were calculated using the finite-difference method using a $12\times 12\times 2$ $\bm{k}$-mesh.

For evaluating the band-edge correction, we performed band-structure calculation using the HSE06 hybrid functional HSE06 including SOC. Here, we adopt HSE06 because it reproduces the experimental band gap of LaOBiS₂ well: the direct band gap at the X ($=\mathbf{a}^{*}/2$) point calculated by HSE06$+$SOC is 0.97 eV, while the experimental optical gap is 1.0 eV LaOBiS2_gap . Band-edge correction was evaluated as an energy difference $\Delta_{\mathrm{bg}}=\epsilon_{\mathrm{HSE06,\ SOC}}-\epsilon_{\mathrm{PBE,\ no\ SOC}}$ for the highest-occupied (for $\Delta_{\mathrm{bg}}^{\mathrm{v}}$) and lowest-unoccupied (for $\Delta_{\mathrm{bg}}^{\mathrm{c}}$) states at the X point, where $\epsilon_{\mathrm{HSE06,\ SOC}}$ and $\epsilon_{\mathrm{PBE,\ no\ SOC}}$ are the orbital energies obtained by HSE06 $+$ SOC and PBE without SOC calculations, respectively. For the HSE06$+$SOC calculation, we used the crystal structures optimized by PBE calculation without SOC for a fair comparison of the orbital energies. We used a $6\times 6\times 1$ $\bm{k}$-mesh for evaluating the band-edge correction. Obtained band-edge corrections are shown in Table 2.

While the band-edge correction was applied in the way described in the previous paragraph, the calculation error coming from PBE still remains in our calculation. For example, it is often said that GGA tends to delocalize the electronic state owing to the self-interaction error, which might result in an inaccuracy of the stability for the local structural change and a localized electronic state around it. Such an inaccuracy would be partially removed by using more sophisticate approximations such as the hybrid functionals for structural optimization, which is an important future issue.

To evaluate the chemical potentials of atoms, $\mu_{i}$ in Eq. (1), we determined possible sets of compounds that can coexist in chemical equilibrium in the following way. Note that the chemical potentials of atoms are not uniquely determined because we can consider different experimental environments such as sulfur-rich and oxygen-poor ones. In this section, we shall show all the possible equilibriums and then choose a representative one for which we shall present defect formation energies from the next section.

We first calculated the total energies of several solids consisting of La, O, F, Sb, Bi, and S atoms. Table 3 in Appendix A shows a list of all the compounds considered here. Next, we found out all possible subsets of the compounds, where every compound $\mathrm{La}_{v}\mathrm{O}_{w}\mathrm{F}_{x}X_{y}\mathrm{S}_{z}$ included in a subset satisfies

while any compound not included in this subset satisfies

Here, the chemical potential of an atom is denoted as $\mu[\mathrm{A}]$ for atom A. Each subset consists of five compounds, which enables consistent determination of the atomic chemical potentials using Eq. (6). Note that this procedure was done separately for $X=$ Sb and Bi. These conditions represent the chemical equilibrium where five compounds included in a subset coexist while other compounds are unstable to exist. For $X=$ Sb, we assumed that LaOSbS₂ is included in every subset. On the other hand, for $X=$ Bi, we found that we cannot include LaOBiS₂ because LaOBiS₂ is slightly unstable with respect to the following reaction,

by 20 meV. This slight instability might be a calculation error or perhaps means that LaOBiS₂ is metastable at zero temperature. It is difficult to discuss such a small energy difference and is outside the scope of this study. Here, instead of including LaOBiS₂ into a subset, we just chose the chemical equilibriums where the instability of LaOBiS₂ is small, by discarding the subsets with $\Delta_{\mathrm{LaOBiS}_{2}}\equiv E[\mathrm{LaOBiS}_{2}]-(\mu[\mathrm{La}]+\mu[\mathrm{O}]+\mu[\mathrm{Bi}]+2\mu[\mathrm{S}])>100$ meV. For both $X=$ Sb and Bi, we also imposed that every subset includes only one fluorine compound. Thanks to this constraint, we can immediately obtain the chemical potentials of La, O, $X$, and S atoms satisfying the equilibrium condition without fluorine, by simply eliminating the fluorine compound from each subset. This procedure yields the same chemical potentials of La, O, $X$, and S atoms regardless of the existence of fluorine, which is convenient to compare defect formation energies in chemical equilibrium with and without fluorine.

All possible chemical potentials of atoms, to say, all possible chemical equilibriums, obtained here are listed on Tables 4 and 5 in Appendix A, for LaOSbS₂ and LaOBiS₂, respectively. Figure 2 summarizes these chemical potentials. For this plot, we defined $\Delta\mu_{i}\equiv\mu_{i}-\mu_{i}^{0}$ for each atom $i$, where $\mu_{i}^{0}$ is the total energy per atom for La (solid), O₂ molecule, F₂ molecule, Sb (solid), Bi (solid), and $\alpha$-S (solid), for La, O, F, Sb, Bi, and S, respectively. As a representative set, we chose the A1 (LaF₃, La₂S₃, La₂O₂S₂, La₂O₂S, and LaOSbS₂) on Table 4 for LaOSbS₂ and the B1 (LaF₃, La₂S₃, La₂O₂S₂, La₂O₂S, and Bi₂S₃) on Table 5 for LaOBiS₂, the chemical potentials for which are shown by red bold circles in Fig. 2. As shown in Fig. 2, these two chemical equilibriums provide similar values of chemical potentials for LaOSbS₂ and LaOBiS₂. While we use these chemical potentials in the following subsections, we again note that the chemical potentials of atoms correspond to the experimental environment, which can be controlled by the experimental setup. Thus, calculated defect formation energies presented hereafter can vary by differences in chemical potentials shown in Fig. 2.

Here we briefly mention the consistency between theoretically estimated chemical equilibriums and actual experimental environments for synthesis. It was experimentally reported that La₂O₂S impurity is often found in the synthesis of LaOBiS₂ (e.g., Ref. [LaOBiS2_gap, ]). LaF₃ is used for introducing fluorine into LaOBiS₂ and is often found as an impurity phase (e.g., Ref. [BiS2_2, ]). Also, Bi₂OS₂, Bi₂S₃, and Bi impurity phases are found in the synthesis of Bi₂OS₂, which is a sibling compound of LaOBiS₂, where La is replaced with Bi in Bi₂OS₂ LaOBiS2_gap . Given these facts, it seems that our theoretical estimate to some extent represents experimental environments for synthesis.

Formation energies of the anion point defects are shown in Fig. 3. In the figure, point defect species are denoted with the Kröger-Vink notation without showing the electronic charge since it varies against the chemical potential. Lines with the most stable $q$ are shown for each chemical potential. For each line, a value of charge $q$ is shown in the figure. Here, the slope of a line in the figure equals $q$ (see Eq. (1)). The equatorial and apical sulfur atoms, S^eq and S^ap, respectively, are defined in Fig. 1. We applied the correction $E_{\mathrm{be}}[D^{q}]$ except V${}_{\mathrm{S}^{\mathrm{eq}}}$, where the impurity energy level is not shallow, as we shall see later in this section.

In Fig. 3, we can see that S_O and O_S are relatively unstable compared with V_O and V_S. For $X=$ Bi, the possible formation of O and S vacancies is consistent with the experimentally observed n-type carriers in the mother compound without fluorine doping. Note that here we adopted the S-rich environment, i.e., a high chemical potential $\mu[\mathrm{S}]$, as shown in Fig. 2, but V_S can be lowered by about 0.6 eV for $X=$ Bi when one adopts the S-poor environment as shown in Fig. 12(c). It is also noteworthy that the oxygen vacancy was experimentally observed in LaOBiSSe LaOBiSSe_VO . For $X=$ Sb, V_O and V${}_{\mathrm{S}^{\mathrm{ap}}}$ become much more unstable than $X=$ Bi, while V${}_{\mathrm{S}^{\mathrm{eq}}}$ becomes much more stable at the conduction bottom with $q=0$ charge. We adopt the S-rich environment here, but the S-poor environment can lower the V_S formation energy by about 0.5 eV as shown in Fig. 12(a), which makes the formation energy of V${}_{\mathrm{S}^{\mathrm{eq}}}$ negative. Also, the formation energy of V_O can be almost zero by adopting the O-poor environment.

A negative formation energy for the whole energy range of the Fermi level within the band gap is troublesome from the viewpoint of the materials stability, because a large number of the point defects will be created in that case. We got such a low defect formation energy in some cases, such as V_O in Figs. 3(b) and 12(c), V_S in Fig. 12(a). One possible reason that makes this situation is the calculation errors, such as the finite-size error and the band-gap inaccuracy, both of which can in general result in an error in the formation energy of an order of 0.1 eV, even if one carefully pays attention to these issues and applies the correction terms. Another possibility is that some sets of the chemical potentials of atoms in Tables 4 and 5 are not appropriate for materials synthesis. Considering the limitation of the accuracy of our calculation method, it is difficult to go into details further. Nevertheless, we again emphasize the consistency between our calculation and experiments where a small amount of electron carriers is observed. If a defect that can introduce the electron carriers into the system has a negative formation energy for some Fermi level (e.g., near the valence band top), the Fermi level is automatically pushed up by the introduced electron carriers. Then, the Fermi level becomes close to the conduction band bottom, while the formation energy should be positive at the conduction band edge for the materials stability. Our calculation results well represent this trend.

Among the point defects investigated in this section, V${}_{\mathrm{S}^{\mathrm{eq}}}$ is interesting in the sense that this defect offers an in-gap impurity state accompanied by a change of the (local) crystal structure. Figure 4 presents the partial DOS and the crystal structure of the $X$S₂ plane that includes the $\mathrm{S}^{\mathrm{eq}}$ vacancy. We found that the crystal structures shown in Figs. 4(b)(f) are stable for $q=2$, while those shown in Figs. 4(d)(h) are stable for $q=0,1$. The characteristic change here is the formation of the $X_{2}$ dimer in the $q=0,1$ structures as shown in Figs. 4(d)(h). Here, the DOS peaks shown with arrows in Figs. 4(a)(c)(e)(g) have a relatively large component of the $X_{2}$ dimer. This tendency is clear for Fig. 4(g) while it is not clear for Fig. 4(e), because the Bi-Bi distance surrounded by the broken line in Fig. 4(f) (4.41 Å) is much longer than that in Fig. 4(h) (3.42 Å). A similar situation happens for the Sb compound, where the Sb-Sb distance is 4.08 Å in Fig. 4(b) and 3.04 Å in Fig. 4(d). The bond formation is usually caused by the electron occupation of the bonding orbitals that have low eigenenergies. In fact, the impurity levels with a relatively large $X_{2}$ dimer components denoted with arrows are occupied in Figs. 4(c)(g) ($q=0$), while it is not the case for the $q=2$ states as shown in Figs. 4(a)(e). The valence electron density for $X=$ Sb with $q=0,2$ shown in Figs. 4(i)(j) also supports the bond formation within such an $X_{2}$ dimer for the $q=0$ case.

It is also noteworthy that the $X$-S network is strongly disarranged in $X=$ Sb with $q=0$ as shown in Fig. 4(d). This notable structural change might be related to the instability of in-plane atoms in BiS₂ compounds observed in experimental studies BiS2_review_thermo2 ; BiS2_atomic_distortion , which seems more significant in the case of $X=$ Sb.

Formation energies of the cation point defects are shown in Fig. 5. We do not show V_La here because we found that its formation energy is too high in the overall range of the chemical potential (e.g., $>4.5$ eV for $X=$ Bi). We applied the correction $E_{\mathrm{be}}[D^{q}]$ except V_X, where the impurity energy level is not shallow, as we shall see later in this section. All the defect formation energies shown in Fig. 5 are positive even at the conduction band bottom.

Here, we briefly discuss how the choice of the atomic chemical potentials can change our conclusion. In Fig. 5, we used an $X$-poor environment (see Fig. 2 and Tables 4 and 5). If one uses an $X$-rich environment, $\mu[\mathrm{Sb}]$ and $\mu[\mathrm{Bi}]$ can be increased by 0.68 and 0.86 eV, respectively, which destabilizes V_X and La_X and stabilizes $X_{\mathrm{La}}$. The energy diagrams for the $X$-rich conditions are shown in Figs. 13(a)–(d). As for the La chemical potential, we can adopt an La-rich environment but it also leads to a high $\mu[X]$ ($X$-rich, see A3 and A5 in Table 4 and B2 in Table 5), which almost does not change the defect formation energies of $X_{\mathrm{La}}$ and La_X, as shown in Figs. 13(a)(c). In total, $X_{\mathrm{La}}$ can exhibit slightly negative formation energy only for an La-poorer and $X$-richer environment, as shown in Figs. 13(b)(d). Such a negative formation energy for the Fermi level in the whole energy region of the band gap is troublesome for the materials stability as discussed in the previous section. This situation is possibly due to the calculation error, such as the finite-size error and the band-gap inaccuracy, or means that that chemical environment is not appropriate for materials synthesis. However, even when the system has a small number of $X_{\mathrm{La}}$ defects, because of the following three reasons, we do not consider that transport properties are much affected by this defect: $X_{\mathrm{La}}$ is a defect in the insulating LaO layer rather than the conducting $X$S₂ layer, this defect introduces no carriers, and we do not find any in-gap impurity state for this defect.

For the V_X defect that shows a change of $q$ in the in-gap energy region, we present the partial DOS for $X=$ Bi with $q=-1,-3$ in Fig. 6. There are very sharp DOS peaks that mainly consist of S^ap just next to the vacancy site. The energy level of this impurity state can lie within the band gap as shown in Fig. 6(b), which induces the change of $q$ in Fig. 5. To say, $q=-3$ gives electron occupation of such defect states as shown in Fig. 6(b).

Formation energies of the point defects of $X_{\mathrm{S}}$ or S_X are shown in Fig. 7. These defects involve an anion-cation exchange, by which we found that many in-gap states take place. Therefore, we did not apply the correction $E_{\mathrm{be}}[D^{q}]$ for all the point defects investigated in this section.

We can see that the defect formation energies of $X_{\mathrm{S}}$ are all high at the conduction band bottom, where the chemical potential is considered to lie in these n-type materials. Thus, we do not go into details about these defects. On the other hand, the defect formation energy of ${\mathrm{S}}_{X}$ is relatively low at the conduction band bottom. Figure 8 presents the partial DOS and the crystal structure of the $X$S₂ plane that includes the S_X defect. In this case, the impurity level shown in Fig. 8(c)(g) mainly consists of the S₂ dimer surrounded by the black broken lines in Fig. 8(d)(h), respectively. In Fig. 8(d)(h), the S₂ dimer is more separate from the surrounding atoms than that shown in Fig. 8(b)(f), as manifested by the bond angle of in-plane sulfurs. The $q=-1$ charge state allows the occupation of the impurity levels, which seems to stabilize the local structures containing the S₂ dimer more separate from the surrounding atoms. It is also noteworthy that the Sb-S network is again strongly disarranged in Fig. 8(b), as we have seen in Figs. 4(b)(d). This also suggests that the SbS₂ plane is prone to be disordered.

Since here we considered the $X$-poor and S-rich environment (see Fig. 2), the ${\mathrm{S}}_{X}$ defect cannot be further stabilized by a different choice of the atomic chemical potentials. The $X_{\mathrm{S}}$ defects can be stabilized to some extent, but still has positive formation energy even if one considers the $X$-rich and S-poor environment, as shown in Fig. 14.

Finally, to see the effects of fluorine doping, we calculated fluorine-related point defects as shown in Figs. 9(a)–(b) for LaOSbS₂ and LaOBiS₂, respectively. We applied the correction $E_{\mathrm{be}}[D^{q}]$ for all the defects considered here. For interstitial fluorine doping, F_i, Figure 10 shows the most stable position of a fluorine atom we found. Note that we also found that the defect formation energy of F_i for the case when one places a fluorine atom between two $X$S₂ layers is almost similar. The energy difference between that case and the most stable case is at most 0.15 eV for $X=$ Sb and 0.02 eV for $X=$ Bi.

It is clearly presented that, for both compounds, the substitutional doping for oxygen, F_O, is the most stable point defect when fluorine is introduced into the crystal. This situation is unchanged by using different sets of chemical potentials of atoms (see Fig. 15). Its charge state is basically F_O^∙ (i.e., $q=+1$) when the chemical potential lies in the band gap, which is naturally understood by considering that O^2- is replaced with F^-. This is consistent with many experimental studies of LaOBiS₂ where fluorine is doped into the crystal to introduce electron carriers BiS2_2 . The defect formation energy of F_O is higher in $X=$ Sb, which is consistent with the previous theoretical study Hirayama . Even so, it is noteworthy that the defect formation energy of F_O is negative both for LaOSbS₂ and for LaOBiS₂. In fact, an experimental study found that lattice parameters are changed and the electrical conductivity is to some extent increased by fluorine doping for LaOSbSe₂ LnOSbSe2 , which suggests that fluorine was successfully doped into the crystal in experiments. In the following subsection, we discuss why carrier control is still difficult in Sb compounds, even though the fluorine substitutional doping seems to be realized.

We note that the defect formation energy of F_O is close to zero at the conduction band bottom, in particular for $X=$ Sb, which can be positive by using different chemical potentials of atoms, as shown in Fig. 15(b). In addition, the calculation results, of course, should be different for different compounds (i.e., different lanthanoid and chalcogen elements). In this sense, our calculation results suggest that fluorine substitutional doping for oxygen can be energetically unfavorable in some $X=$ Sb compounds or some experimental setup.

Here, we summarize the important differences between $X=$ Sb and Bi compounds: (i) V_O has a lower formation energy in $X=$ Bi than in $X=$ Sb. (ii) In-plane structural instability seems to result in the relatively stable V${}_{\mathrm{S}^{\mathrm{eq}}}$ for $X=$ Sb (and S_X). The in-plane $X$-S bonding network is more strongly disarranged for $X=$ Sb. (iii) F_O has a lower formation energy in $X=$ Bi than in $X=$ Sb.

Both (i) and (iii) suggest that the electron carriers are much more difficult to introduce because of higher formation energies of corresponding point defects in $X=$ Sb compounds than in $X=$ Bi compounds. Even when one succeeds in introducing the electron carriers, the strong in-plane structural instability in $X=$ Sb can be an obstacle for realizing high electrical conductivity, as suggested by (ii). While we have only investigated the point defects in this study, several possible higher-order defects, such as a twin defect, can occur in these compounds. It is an important future problem to clarify the origins of the relatively low electrical conductivity in $X=$ Sb compounds from this viewpoint.

Our calculation suggests that the anion vacancy (V_O or V_S) can be relatively stable in some chemical environment. This observation means that oxygen-rich and sulfur-rich environment can be beneficial to synthesize the crystal with a low concentration of defects. It also means that one can possibly control the electron carrier concentration by changing the chemical potentials of oxygen and sulfur. When one would like to introduce a large amount of electron carriers, the fluorine substitutional doping F_O is shown to be effective as is well known in experimental studies of BiS₂ compounds.