Towards intermediate-band photovoltaic absorbers: theoretical insights on the incorporation of Ti and Nb in In2S3

Creation of a partially filled intermediate band in a photovoltaic absorber material is an appealing concept for increasing the quantum efficiency of solar cells. Recently, we showed that formation of a partially filled intermediate band through doping a host semiconductor with a transition metal dopant is hindered by the strongly correlated nature of d-electrons and the antecedent Jahn–Teller distortion, as we have previously reported. In present work, we take a step forward and study the delocalization of a filled (valence-like) intermediate band throughout the lattice: a case study of Ti- and Nb-doped In2S3. By means of hybrid density functional calculations, we present extensive analysis on structural properties and interactions leading to electronic characteristics of Ti- and Nb-doped In2S3. We find that Nb creates an occupied doublet, which can become delocalized onto the crystal at high but feasible concentrations (around 2.5 at% and above). As a consequence, doping In2S3 with adequately high concentrations of Nb allows the subgap intermediate band to conduction band absorption, which leads to higher photocurrent densities compared to pure In2S3. Ti on the other hand forms an occupied singlet intermediate band, which remains strongly localized even at high concentration of 5 at%.


INTRODUCTION
The idea of multi-transition solar cells was first introduced by Wolf in 1959 1 , 5 years after manufacturing of the first silicon p-n junction photocell. For a better utilization of the solar spectrum, he suggested to introduce deep impurity levels within the absorber band gap, e.g., by means of transition metal doping 2 . These levels would facilitate sequential absorption of lower-energy photons from the valence band (VB) to empty defect levels and from the filled defect levels to the conduction band (CB). However, since the intermediate impurity levels trigger the occurrence of Shockley-Read-Hall (SRH) non-radiative recombinations 3,4 , absorbed photons do not necessarily produce electron-hole pairs, and therefore, the idea did not evolve further until some decades later. In 1997, Luque and Martí 5 demonstrated theoretically that outperforming the Shockley-Queisser (SQ) efficiency limit 6 of solar cells is achievable by introducing a partially filled intermediate band (IB) instead of an isolated defect level inside the band gap. The core of their idea is that through formation of partially filled IBs, two photons with sub-band gap energies can excite electrons and generate electron-hole pairs. Furthermore, Luque et al. 7 argued that the non-radiative recombination can be suppressed if wavefunctions of the impurities overlap. Such condition can be achieved at high dopant concentrations.
As a result of the two-step photon absorption associated with the presence of a partially filled IB, the efficiency of a solar cell can approach 63.1% instead of 40.7% being the SQ limit for single junction solar cells 6 . In order to design an intermediate-band solar cell (IBSC) with such high efficiency, one has to take into account both the magnitude of the band gap and the position of the introduced IB. For example, the maximum efficiency of 63.1% corresponds to a host semiconductor with a band gap of 1.95 eV and a partially filled IB that divides it into 1.24 eV and 0.71 eV sub-gaps 5 .
In recent years, several attempts have been made to find doped semiconductors, which can form an IB. Chalcopyrite compounds belonging to the Cu(In,Ga)(S,Se) 2 (CIGS) family were the first thin film materials investigated as candidates for formation of IBs 8,9 . Besides, indium thiospinels are another class of semiconductor hosts with well-known thin-film deposition techniques, and recognized potentials for the development of an IBSC. Among these materials, β−In 2 S 3 represents the most promising candidate for hosting intermediate bands when doped with transition metals 10 , since its band gap of 2.03 eV 11 is very close to that of the ideal IB semiconductor. Furthermore, its advantage as possible semiconductor host for the development of an IBSC when compared to CIGS is that the thiospinels have octahedral cation sites, providing the opportunity of including transition metals in a more stable six-fold environment. β−In 2 S 3 is an n-type semiconductor that crystallizes in a defective spinel-like structure 11,12 . Its optoelectronic properties have already made it an attractive material for photovoltaic applications, namely, as buffer layer in CIGS-based thin-film solar cells [13][14][15] .
The study of impurities in semiconductors requires the application of relatively large supercells. In this respect, (semi) local functionals with moderate computational cost became the main workhorse for investigation of doped semiconductors. Nevertheless, the theoretical interpretations of physics of such systems were affected by two well-known limitations of the (semi) local functionals: the delocalization and the band gap errors. The local density and generalized gradient approximations to density functional theory (DFT) lead to considerable errors in describing nonlocality of the screened exchange interaction and consequently the electrostatic self-interaction. The consequence of this is that an electron interacts with itself, which artificially favors its delocalization. This results in particularly sizable errors for strongly localized states, e.g., the d states of transition metals (TM) 16 . Furthermore, the band gap underestimation by (semi)local functionals leads to inaccurate description of electronic properties and thermodynamic charge transition levels of such systems. One can largely surmount the above-mentioned problems by using hybrid density functionals with a tuned mixing parameter or alternatively a tuned range separation parameter [17][18][19] . In case of GaN:Mn, for instance, standard DFT predicts partially filled threefold degenerate states inside the band gap, whereas the rangeseparated hybrid functional shows splitting of t 2 states into an unoccupied singlet and an occupied doublet, in excellent agreement with GW 0 calculations and experiments. We reported a similar scenario for In 2 S 3 :V in our previous work 20 . Although the electronic and thermodynamic properties produced by hybrid functional are precise for many TM-doped systems, but it does not fulfill the generalized Koopmans theorem for every TM-doped semiconductor 21,22 . In such cases (e.g., Cr-doped SiC 21 ), the exact Hartree-Fock exchange energy does not fully cancel the selfinteraction error and hence the generalized koopmans' condition (gKC) 23,24 is not satisfied. In order to address this issue properly, we examine here the fulfillment of the gKC by our applied hybrid functional (HSE06). We find that for both Ti and Nb donors, the HSE06 functional shows a negligible deviation from the gKC. Hence, our obtained results are self-interaction-free.
The aim of this work is to present a careful, detailed theoretical picture of the electronic structure of Ti-and Nb-doped In 2 S 3 and to describe under what circumstances a higher absorption of solar spectrum can be achieved in this material. This paper is arranged as follows. Results section concentrates on inclusion of Ti and Nb (which we label as 'TM' hereafter) on In tetrahedral and octahedral substitutional sites and provides a comparison between their formation energies. Then, the electronic structure of TM-substituted In 2 S 3 is discussed. Thereafter, the effect of TM concentration on the formation of an impurity band and mixing enthalpies of the solutes are discussed. A complete description of the computational setup and the theory are given in methods section.

RESULTS
Formation energy of TMs in In 2 S 3 There are three nonequivalent crystallographic environments for In and S atoms in the β-phase of In 2 S 3 . The In substitutional sites in the crystalline matrix of In 2 S 3 are tetrahedral 8e (In 8e ), octahedral 8c (In 8c ), and octahedral 16h (In 16h ) sites 12 , corresponding to Wyckoff notations of space-group 141. Both, In 8c and In 16h sites are six-fold coordinated by S and are different in their associated In-S bond lengths. In 8e occupies two-thirds of the available tetrahedral cation sites, while the remaining one-third tetrahedral sites are unoccupied. Sulfur atoms are arranged on distorted cubic closed packed lattices, where the S positions are distinct due to the different occupation of the surrounding In sites. The atomic structure of the β−In 2 S 3 unit cell is illustrated in Supplementary Fig. 1. Figure 1 shows the calculated formation energy of TMs on 8e tetrahedral and 16h octahedral sites of In for In-rich and S-rich growth atmospheres. In this set of calculations, one TM was added to the crystalline matrix of In 2 S 3 . The formation energies are depicted from E F = 0 eV (corresponding to p-type material) to E F = 2.1 eV (corresponding to n-type material), covering the whole range of possible Fermi energies. The points, where the line slop changes represent transitions of the stable charge state of the impurity. Since the formation energies of TMs on In 8c and In 16h octahedral sites are close, here, we only show formation energies on In 16h sites. Here, our goal is to identify the most favorable defect configuration and the charged data are shown for completeness. In order to discuss Ti and Nb as charged defects, one needs to consider the possible charge compensation mechanisms. Since the concentration of TMs are quite high, the charge compensation by means of intrinsic defects is extremely unlikely. Hence, oxidation of Ti and Nb would require presence of abundant acceptor dopants in the system. Since In 2 S 3 is an n-type material, Ti and Nb can be essentially considered as solutes, i.e., in neutral charge state. Therefore, the only relevant finite size corrections would be due to the elastic interactions, which we safely ignored because the considered concentrations are comparable to what is experimentally achievable. Clearly, the calculated energies can not be used to predict equilibrium concentrations.
For both, In-and S-rich samples, TMs prefer to occupy In octahedral sites rather than tetrahedral sites. Due to its very low formation energy (0.06 eV at CBM), Ti In 16h can exist in In-rich samples in very high concentrations, whereas the formation energy of Ti In 16h in S-rich samples is one order of magnitude larger (0.59 eV at CBM). The formation energy of Nb In 16h is 0.68 and 1.22 eV for In-rich and Srich conditions, respectively. Intuitively, one would expect a lower formation energy of Nb/Ti under In-poor conditions. However, one has to consider that here the thermodynamic constraints are on the possible boundary phases of the solutes and the substituted elements, which leads to the coupling of their chemical potentials Eqs. (5), (6), and (7). Therefore, with the increase of the chemical potential of In under In-rich conditions, the chemical potential of the solute increases as well. The origins of the different formation energies of Ti and Nb will be discussed in subsection "Heat of solution of Ti and Nb in In 2 S 3 ". Electronic properties of In 2 S 3 alloyed with TMs The calculated partial and total densities of states (DOS) of pure In 2 S 3 is shown in Fig. 2. The calculated electronic band structure of β−In 2 S 3 with hybrid functional is given in ref. 25 . The DOS at the uppermost part of the valence band is mainly derived from the S 3p states mixing with In 4d states. Due to the small contribution of S 3s states, they are not shown here. Figure 3 shows HSE spinpolarized densities of states of Ti-and Nb-substituted In 2 S 3 , where the energy level of 0 eV corresponds to the highest occupied state. Here, four Ti or Nb ions were either accommodated in four In 16h or In 8c octahedral sites, corresponding to 5% concentration of these solutes. The subsequent incorporation sites are separated from each other by about 8.6 Å, which is the largest distance achievable when using the chosen supercell. It is found that the energetic position and shape of the DOS of the impurity bands are different when the solute TMs substitute In 16h and In 8c sites. As mentioned earlier, there are three nonequivalent crystallographic environment for S atoms in the structure of β−In 2 S 3 . Here, we refer to them as S1, S2, and S3. S1 atoms are four-fold coordinated with two In 8c , one In 16h , and one In 8e next neighbors; S2 atoms are three-fold coordinated with two In 16h and one In 8c next neighbors; and S3 atoms are four-fold coordinated with three In 16h and one In 8e next neighbors. TMs on In 16h sites are bonding with one S1, two S2, and three S3. While the ones on In 8c sites are bonding with four S1 and two S2 atoms. Therefore, the lattice distortions caused by TMs on In 8c or In 16h sites are different.
In 2008, Palacios et al. 26 studied the electronic structure of Vand Ti-doped In 2 S 3 using DFT. For Vanadium, they found a narrow half-filled intermediate band isolated from the VB and CB of In 2 S 3 . Recently, we revisited V-doped In 2 S 3 20 , using hybrid functional and explained that the erroneous DFT-predicted partially filled IB was caused by neglecting the strongly correlated nature of delectrons and the antecedent Jahn-Teller effect. In fact, trivalent V ions do not form any IB at all. And only if V ions get reduced, then, the divalent V introduce a totally-filled triplet IB within the band gap of In 2 S 3 20 . For Ti, they showed that when Ti replaces In octahedral site (the different In octahedral sites were not distinguished in their study), the GGA approximation predicts an overlap between the d-bands of Ti and the CB of In 2 S 3 , which yields to an IB positioned at the bottom of the conduction band. They assumed that correcting the GGA band gap would shift the position of the CBM upward and the incorporation of Ti would eventually lead to a partially filled IB inside the gap. As illustrated in our previous work 20 , going from semi-local PBE to HSE functional, the band gap opens from 0.84 eV to 2.03 eV, and the downward shift of the VBM is about twice larger than the upward shift the CBM, which is in contrast to the assumption of Palacios et al. We see from Fig. 3 that when Ti atoms occupy In 16h sites, an occupied singlet state occurs at 0.14 eV above the VBM and an unoccupied doublet state appears in the CB. On the Ti In 8c configuration; however, no in-gap state appears and only sulfur p states near the band edges move into the gap and reduce the band gap to 1.57 eV. In case of Nb, Fig. 3 shows that inclusion of Nb on In 16h and In 8c sites leads to occupied doublet in-gap states at 0.62 eV and 0.45 eV above the VBM, respectively. The electronic band structure of In 28 S 48 : Ti 16h 4 and In 28 S 48 : Nb 16h 4 are illustrated in Supplementary Fig. 2.
In the octahedral geometry of In 2 S 3 , the six neighboring ligands are oriented along the z-axis and the diagonal directions of the xy plane. Therefore, the d z 2 and d xy orbitals are subject to a stronger electrostatic repulsion compared to the d xz , d yz , and d x 2 Ày 2 orbitals. Consequently, transition metal substituents occupying these sites experience a corresponding crystal field that pushes the d xz , d yz , and d x 2 Ày 2 orbitals (also called t orbitals) into a lower energy compared to the d z 2 and d xy orbitals (also called e orbitals). In addition, due to the spin splitting, these levels further split into spin-up (t + and e + ) and spin down (t − and e − ) levels. As a result of crystal field and exchange-splitting, d levels of free Ti 3+ and Nb 3+ have respectively [t 1 þ e 0 þ t 0 À e 0 À ] and [t 2 þ e 0 þ t 0 À e 0 À ] electronic configurations.  When entering the octahedral In site, e + and e − levels of TM hybridize with S dandling-bond states (e(p)) and create two types of levels with e symmetry: one with a high contribution of six neighboring S atoms, referred to as dangling-bond hybrid (DBH) states and the other mostly localized on the TM atom, referred to as crystal field resonance (CFR) states 27 . Note that the e DBH and e CFR are bonding and antibonding levels, respectively. For the t + and t − levels, however, there is no host level with t symmetry in the similar energy range. Therefore, the t(d) states remain largely uncoupled and the amplitude of wavefunction is mainly localized on the TM, hence a CFR state 27 . Consequently, the order of electronic states of Ti and Nb on octahedral In sites are as: We note that realizing an ideal IB alloy depends on the position and occupancy of the t þ CFR state. The simplified electron-filling schemes of Ti and Nb on vacant In site (Vac 3À In ) are shown in Fig. 4. We see that for both Ti and Nb, the partially filled t CFR þ state splits into a filled lower-energy and an empty upper-energy branch, i.e., t CFR þ ! a 1 þ þ e 0 þ for Ti and t CFR þ ! e 2 þ þ a 0 þ for Nb. The a and e notations refer to the singlet and doublet natures of these states, respectively. Note that the partially filled degenerate states are Jahn-Teller active centers, which necessiates removing the degeneracy of the intermediate state through splitting it into filled and empty states. Therefore, the onset of a partially filled IB which can trigger simultaneous VB → IB and IB → CB transitions is extremely unlikely. However, if well delocalized, the a 1 þ and e 2 þ ingap states can increase the sub-band gap photon absorption, mostly through assisting in the IB → CB transition.
In 2011, Ho 28 used chemical vapor transport method to grow Nb-doped In 2 S 3 films and detected below gap transitions of 1.52 and 1.42 eV. Since these transitions were not detected in undoped In 2 S 3 , he assigned them to the creation of an IB. However, it was not clear whether the observed transitions stemmed from a deep defect level or a delocalized band. Thus, the question is whether the observed intermediate levels have an unwanted localized or a beneficial delocalized character. According to our calculations, In 2 S 3 :Nb allows IB to CB transitions of 1.55 and 1.66 eV, as 5% Nb ions occupy In 16h and In 8c sites, respectively. We note that although the spontaneous VB to IB and IB to CB transitions (which was the primary goal of constructing IB solar cells) can not be achieved, but if the wavefunctions of doped TM become delocalized throughout the lattice, we can have IB to CB subband gap transitions in addition to the VB to CB transitions.
Note that as we move from left to the right along the 3d row of the periodic table, we approach the heavier 3d metals, whose atomic t and e orbitals are deeper in energy. Due to this simple rule, by increasing the atomic number along the same row, the interactions between the inserted TM and anion dangling bonds become stronger. Therefore, the t orbital of the 3d transition metals after Ti should reside below the VBM. As discussed in our previous work, the t + state of the trivalent V splits into e + and a + sublevels, which reside at the top of VB and bottom of the CB, respectively. Hence, no in-gap state forms 20 . Going down the group 5B of periodic table from V to dubnium (Db) (or in general all d-block groups), the energy of d orbitals increases. Consequently, the t and e energies of the heavy Db should reside far above the energy levels of the light V.
Effect of concentration on the intermediate states formed by the incorporation of Ti and Nb As mentioned before, at low concentrations, transition metal impurities are likely to introduce undesired deep defect levels inside the band gap of the host semiconductor. However, at sufficiently high concentrations, their wavefunctions become delocalized and a beneficial IB forms. In the case of In 2 S 3 , we know from the previous subsection that Ti and Nb impurities induce respectively filled singlet and doublet states inside the band gap of In 2 S 3 . However, up to this point it is not possible to predict whether such states are localized defect levels or beneficial delocalized IBs. In order to gain insights on this issue, we present in Figs. 5 and 6 the concentration-dependent    6 Concentration-dependent DOSs of Nb in In 2 S 3 and their corresponding charge density isosurfaces. Spin-polarized densities of states computed for Nb-substituted In 2 S 3 along with the charge density isosurfaces associated with the features induced in the band gap for a, b one, c, d two, and e, f four impurities inside the supercell. The Nb atoms are located on the available 16h In sites in the cell and the pink box in density of states plots points out the defect states for which the correspondent charge density isosurfaces are produced using VESTA 41 (displayed at 3% of their maximum value). The Nb impurities are shown in red and marked with an arrow of the same color. In and S are shown as green and blue spheres, respectively. E. Ghorbani et al. spin-polarized DOS of the In 2 S 3 doped with Ti and Nb ions that are located on the In 16h sites. Furthermore, we show the partial charge density isosurfaces associated to the defect-induced states as the concentration of the doped TM increases. Since the charge density is associated to the defect wavefunction, its degree of (non) localization is associated to creation of localized or delocalized defect-induced states.
As we can see in Fig. 5a, c, e, the Ti-induced features on the electronic structures of In 2 S 3 , highlighted by the dashed box in each figure, do not show variations on their width, shape, or position as the Ti concentration increases, and only their height changes. Furthermore, their narrow shape points out that even at higher concentrations, the Ti-induced defect states would remain localized. In order to confirm this conclusion, we rely on the charge density isosurfaces of the Ti-induced defect states shown in Fig. 5b, d, f. As shown, even for a concentration of 5%, there is a rather negligible overlap between the states induced by Ti ions and host ions. This means that Ti-doped In 2 S 3 , at least for doping concentrations up to 5%, would not favor the formation of a beneficial IB.
On the other hand, as shown in Fig. 6a, c, e, Nb-induced states in In 2 S 3 are concentration dependent. When changing the doping concentration from 1.25 to 2.5%, the situation is analogous to the case of Ti and only the height of the defect states varies. However, when the Nb concentration reaches a value of 5%, the shape of the defect-induced states changes, too. Carrying out the same charge density analysis done for the case of Ti, we obtained the results presented in Fig. 6b, d, f. As illustrated, at the concentration of 2.5 and 5%, we see that the overlap of the defect wavefunction is achieved in all three dimensions and the introduction of a delocalized IB becomes probable. This result agrees with the wider intermediate band observed in the DOS (Fig. 6) and band structure ( Supplementary Fig. 2) of the Nb-doped In 2 S 3 , which amounts to 0.39 eV, compared to 0.15 eV obtained for the Ti case. Therefore, the potentially relevant Nb-doped In 2 S 3 with its well positioned IB can exhibit beneficial features for below gap IB → CB transitions.
In 2017, Wägele et al. 29 used X-ray diffraction, Raman spectroscopy, and scanning electron microscopy to analyze structural properties of the V-doped In 2 S 3 for various concentrations of V. They found that even with excessive amount of V (as high as 5.8 at %) in In 2 S 3 , no additional peak appeared in the X-ray spectra and only the lattice constant decreased, which is due to difference in atomic radii of V (1.35 Å) and In (1.55 Å). V and Nb are chemically similar elements and the atomic radii of Nb (1.45 Å) is closer to that of In. In this regard, we expect that In 2 S 3 is also stable against Nb incorporation.
The reader must be warned that the onset of a filled IB shifts the quasi-Fermi level to the top of the IB, which consequently limits the open circuit voltage (V OC ) of the device. Therefore, a key area of future research would be to identify possible mechanisms of depopulating the filled IB. Photo-filling or -emptying of an IB proposed by Strandberg and Reenass 30 could be helpful here as well. We note that the filled IB can be used on a device with spatially varying occupation of intermediate-band states. In such a device, the absorber has a filled IB in half of the device and an empty IB on the other half 31 . The theoretical efficiency of a device with spatially decoupled VB-to-IB and IB-to-CB is 52.83%, which is lower than the efficiency of an ideal partially filled IB but still significantly higher than the efficiency of the best single band gap solar cell. Based upon current state-of-the-art in the field, below gap absorption of Nb-doped In 2 S 3 can be also further improved through an effective non-compensated pn-doping approach, which has been successfully applied in TiO 2 32 .
Heat of solution of Ti and Nb in In 2 S 3 The heat of solution (H s ) of TMs in In 2 S 3 were calculated according to Eq. (7) and are shown in Fig. 7. Our results indicate that there is a difference between H s of TMs in In-rich and S-rich samples. For both Ti and Nb, substituting an In site is energetically more favorable under In-rich condition and their inclusion on S-rich samples becomes energetically more expensive. In particular, we see that the heat of solution of Ti atoms substituting the In sites are lower than Nb. To understand the difference in solubilities of Ti and Nb in In 2 S 3 , we analyze the mechanical and chemical contributions to the calculated heat of solution 33 . The mechanical (or structural) effect refers to the energy stored due to the size misfit of the impurity in the host crystalline matrix. The chemical effect represents the energy stored or released due to the constitution of impurity-host bonds. Depending on the nature of the bonds, this term can be positive or negative. According to Lozovoi et al. 33 , we define configurations A, A v , B, and B v as follows:  Note that A v and B v are hypothetical structures that we construct through using the relaxed A and B structures, respectively.
The change in energy due to structural effect consists of two steps; Host Removal (HR) and Structural Substitution (SS), which we define as following: (1) The contribution of Chemical Effect (CE) to heat of solution is defined as: for Ti, and for Nb. Table 1 lists the mechanical and chemical contributions to the heat of solution of Ti and Nb in In 2 S 3 . The energy values correspond to the incorporation of 1.25% TM in In-rich samples. We see that the removal of an In atom from the host configuration (the H HR mechanism) accounts for a large increase in the energy of the system. The effect of structural substitution (H SS ) is modest for Nb (0.01 eV). In case of Ti, however, there is a considerably larger lattice distortion contributed to the SS mechanism, which is due to the smaller size of Ti compared to Nb. Another distinguishing feature can be observed when the Ti (or Nb) are back into the vacant site (H CE ). Since both Ti and Nb are more electropositive than In, they form stronger bonds with the six surrounding S atoms and the term H CE causes a large decrease in the energy of the system containing the alloying TM. For Ti, contribution of H CE in decreasing the heat of solution is −8.33 eV. As a consequence, the chemical decrease almost cancels the mechanical increase (sum of H HR and H SS mechanisms) in the total H s . This results in a very low, but still positive heat of solution for Ti (0.08 eV). Note that the mechanical increase still dominates the chemical decrease, resulting in a positive heat of solution for Ti. The bonding between Nb and S is weaker compared to Ti. Therefore, the H CE is less negative for Nb and the energy consuming host removal and structural substitution can not be cancelled out by forming the weaker Nb-S bonds. Hence, the heat of solution of Nb is larger compared to Ti.
Our results indicate that the heat of solution, H s , of Ti and Nb are positive in both In-rich and S-rich atmospheres, suggesting that their dissolution into In 2 S 3 is an endothermic reaction. Therefore, both In-rich and S-rich samples can accommodate Ti and Nb atoms on octahedral In sites without leading to the formation of secondary phases.

DISCUSSION
By using first principles calculations, we analyzed structural, energetic, and electronic characteristics of Ti and Nb-doped In 2 S 3 . We showed that in n-type In 2 S 3 both Ti and Nb exist as neutral impurities (i.e., assuming the +III oxidation number) on In substitutional sites. The results show a tendency of Ti and Nb to occupy 16h octahedral In site rather than 8e tetrahedral site. In the octahedral symmetry of 16h In site, the d z 2 and d xy orbitals (i.e., the e branch) of the TMs experience a stronger repulsion compared to the d xz , d yz , and d x 2 Ày 2 orbitals (i.e., the t branch). Consequently, the d z 2 and d xy orbitals are shifted to higher energies, forming doubly degenerate e(d) states inside the CB of the In 2 S 3 for the two spin channels. The spin-up t branch (t + ), however, has lower energy and can be found inside the band gap. For Ti on 16h In site, the t + state consist of one occupied and two unoccupied levels. For Nb, however, the t + state accommodates two electrons. Hence, consisting of two occupied and one unoccupied levels. Due to the influence of strong electron correlation effects, the filled and empty states of the TM t + branch can not exist at the same energy level. Therefore, the threefold degeneracy of such level is removed due to an electron Jahn-Teller effect, which splits it into occupied and unoccupied sublevels. Prediction of such splittings depends on a correct description of the exchangecorrelation functional, which can be tested through fulfillment of the generalized Koopmann's theorem. Our choice of HSE functional (i.e., α = 0.25 and w = 0.13 Å −1 ) demonstrates an almost complete cancellation of the self-interaction error for Tiand Nb-doped In 2 S 3 . The results clearly indicate that Ti t + state splits into an occupied singlet at 0.14 eV above the VBM and an empty doublet inside the CB. The Nb t + state, however, splits into an occupied doublet at 0.62 eV above the VBM and an unoccupied singlet in the CB. We have demonstrated that the occupied singlet in-gap state of Ti 16h remains highly localized even at high concentration of 5%. The occupied doublet of the Nb t + , however, shows a considerable overlap with the interacting orbitals of the neighbouring atoms at concentrations above 2.5%. This finding suggest that Nb-doped In 2 S 3 can show below gap IB to CB transition. However, the VB to IB transition is prevented due to the absence of an empty state inside the gap. A comparison between the solubility of Ti and Nb in In 2 S 3 shows that under similar growth conditions, Ti can exist in In 2 S 3 at higher concentrations. Finally, we note that current study can be continued by analyzing the structural and electronic properties of other transition metals such as zirconium, hafnium, and tantalum in In 2 S 3 .

Computational details
Calculations in the framework of DFT were carried out using the Vienna ab initio Simulation Package 34,35 . The results were obtained using screened hybrid functional of Heyd, Scuseria, and Ernzerhof (HSE) 36 . The contribution of Hartree-Fock exchange, α, was set to the standard value of 0.25 and the screening range of the electron interaction was treated by adjusting the Thomas-Fermi screening parameter, w, to 0.13 Å −1 , leading to a band gap of 2.03 eV 13 , which is in excellent agreement with experimental reports 11,37 . The interactions between the ionic cores and the valence electrons were described with the projector augmented-wave (PAW) approach 38,39 . The semi-core 3p, 4s, and 3d states of Ti and 4p, 5s, and 4d states of Nb were treated as part of the valence. The wavefunctions were expanded up to a cutoff energy of 400 eV and the Brillouin zone was sampled using a 4 × 4 × 1 k-point mesh. In all calculations the spinpolarization effects were taken into account. The convergence thresholds for electronic self-consistency and residual Hellmann-Feynman force component on each atom were set to 10 −5 eV and 0.05 eV/Å, respectively. Calculations were performed for a unit cell containing 80 atoms. Both the ionic positions and the cell metric were allowed to relax.
In Table 2, we present bulk properties of the pure and TM-incorporated In 2 S 3 . A comparison between the bulk properties of the TM-incorporated on the In 8c and 16h sites and the PBE and HSE functionals is given in Supplementary  45 Å). Therefore, as expected, substitution of In by these transition metals reduces the size of unit cell. Furthermore, the band gap of In 2 S 3 slightly increases upon inclusion of the TM atoms, which is expected by considering the contraction of the cells.

Fulfillment of the generalized Koopmans' condition
As mentioned in the introduction, to assess suitability of HSE functional, one needs to test the fulfillment of the gKC. For this purpose, we quantify the non-Koopmans' energy (E NK i ) using where ϵ i is the highest occupied state in the system of N electrons, and E N and E N−1 are the total energies of the N-and (N − 1)-electron systems, respectively. In case of E NK i ¼ 0, the gKC is fulfilled and the applied functional presents a self-interaction-free description of the orbitals. We find that using the HSE functional, the E NK i is negligibly small for both Ti and Nb donors (smaller than 0.05 eV). Therefore, the applied hybrid functional correctly describes electron correlation in the system.

Calculation of formation energy and heat of solution
The defect formation energies were calculated at the theoretical equilibrium volume as where q is the charge state, ΔE[d q ] is the difference of the total energy of the cell with and without impurity, μ TM and μ In are the chemical potentials of the TM and In atoms, respectively. E VBM is the energy of the valence band maximum (VBM) of the pure cell, and E F is position of the Fermi level with respect to the E VBM . Δv q is a term aligning the average electrostatic potential of the pure and defect cell containing a neutral defect.
For n dopants in the neutral charge state, the defect formation energy is equivalent to the heat of solution H s,TM : H s;TM ¼ E TMn In 32Àn S 48 À E In32 S 48 À nμ TM þ nμ In : E TMn In 32Àn S 48 is the total energy of the defective supercell of indium sulfide containing various TM concentrations. The term E In32 S 48 accounts for the total energy of a non-defective supercell of In 2 S 3 with the same size as the one used for the defective systems, while n is the number of atoms added or removed from the supercell.
In Eqs. (6) and (7), the chemical potential of the incorporated TM (substituted In) is given by is the total energy of the TM (In) atom in its reference phase. Δμ TM (Δμ In ) is the maximum possible change of the chemical potential of the TM (In). For staying in the stable phase field of In 2 S 3 , the following constraints must hold: Δμ In 0; Δμ S 0; (8) 2Δμ In þ 3Δμ S ¼ ΔH f ðIn 2 S 3 Þ; where ΔH f (In 2 S 3 ) is the formation enthalpy of the β−In 2 S 3 . Here, we consider two experimental growth conditions: (i) In-rich limit: In 2 S 3 is in equilibrium with the metallic tetragonal Inphase. This determines the upper limit of μ In , Δμ In = 0, and the lower limit of μ S , Δμ S ¼ 1 3 ΔH f ðIn 2 S 3 Þ. (ii) S-rich limit: In 2 S 3 is in equilibrium with the orthorhombic S-phase, whence Δμ S = 0 and Δμ In ¼ 1 2 ΔH f ðIn 2 S 3 Þ. Note that the chemical potentials of In and S can vary between these two limits. Furthermore, the concentrations of Ti and Nb in In 2 S 3 are restricted by the formation of TiS 2 and NbS 2 as competing phases. Hence, the chemical potentials of Ti and Nb must satisfy the relations Eqs. (10) and (11), respectively, Δμ Nb þ 2Δμ S ΔH f ðNbS 2 Þ: The formation enthalpies of In 2 S 3 , TiS 2 , and NbS 2 were calculated with the HSE functional and are listed in Table 3. The chemical potentials of Ti and Nb are referenced to the energies of their pure elemental phases, i.e., hexagonal-close-packed Ti and Nb.

DATA AVAILABILITY
Data are available upon reasonable request from the corresponding author. Lattice parameters a and c (in Å) and band gap of pure and TM-substituted β−In 2 S 3 calculated with HSE functional. Experimentally reported values are given in parentheses.