**1. Introduction**

Semiconductor nanocrystals, known as quantum dots (QDs), have received significant attention in recent decades due to their potential relative to a multitude of applications including bioimaging [1–3], light-emitting diodes [4,5], and photovoltaic devices [6,7]. The essential property of these materials is the quantum size effect, whereby the band gap becomes strongly size-dependent when the radius of the nanocrystal is smaller than the excitonic Bohr radius. Early QD materials of interest included binary heavy metal chalcogenides such as CdS, CdSe, CdTe, PbS and PbSe. However, over the past decade, the improvement of solvothermal synthesis techniques has made it possible to produce high quality multinary chalcogenide nanocrystals which do not contain toxic elements [8]. By tuning the reactivity of the different precursors, a high degree of control over the stoichiometry can be achieved, enabling the QD properties to be controlled by varying the composition as well as by modifying the size. For example, CuInS2 QDs have been developed as promising candidates for various applications. Their band gap can be tuned between 1.5 eV (i.e., the band gap of bulk CuInS2) and approximately 2 eV, corresponding to photoluminescence emission in the near infrared and yellow regions of the spectrum, respectively [9]. The ability to emit within the so-called biological window makes CuInS2 QDs ideal fluorescent probes for bio-imaging [10–12], and their strong, broad absorption within the solar spectrum makes them a viable photoactive component in photovoltaic devices [13]. Although CuInS2 nanocrystals are less toxic than their heavy metal chalcogenide counterparts, there are some considerable toxicity concerns regarding indium, and furthermore, indium is a rare and expensive metal. Therefore, alternatives containing relatively earth-abundant constituent elements such as Zn, Sn and Fe have been investigated [14,15]. In this study we focus on the replacement of indium with two earth-abundant and inexpensive elements, namely Al and Fe.

**Citation:** Dickens, C.; Kinsella, A.O.J.; Watkins, M.; Booth, M. The Presence of Charge Transfer Defect Complexes in Intermediate Band CuAl1−*p*Fe*p*S2. *Crystals* **2022**, *12*, 1823. https:// doi.org/10.3390/cryst12121823

Academic Editor: Jacek Cwik ´

Received: 1 November 2022 Accepted: 10 December 2022 Published: 14 December 2022

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

CuAlS2 is similar to CuInS2 in that it is a [A<sup>I</sup> ][BIII][CVI]2 compound that is stable in the chalcopyrite crystal structure but has a much wider band gap (approximately 3.5 eV). Early computational studies by Zunger et al. found that CuAlS2 is further from the ideal zincblende structure than CuInS2, and that the interactions between the Cu 3*d* and S 3*p* states are stronger in CuAlS2, leading to a narrower upper valence band (VB) [16]. Relatively little work has been conducted on CuAlS2 nanocrystalline materials in comparison to CuInS2. Bhattacharyya et al. demonstrated that the type II band offset between CuAlS2 and ZnS or CdS in CuAlS2/ZnS [17,18] and CuAlS2/CdS [19] core/shell nanocrystals can be engineered to enable the effective band gap of the heterostructures to span the visible spectrum.

Chalcopyrite, CuFeS2, is perhaps the most earth-abundant copper ore and, despite having been known for thousands of years, remains to be fully understood, since transition metals with partially filled *d*-orbitals tend to introduce complex properties when included as a constituent in compounds such as chalcogenides. Neutron diffraction investigations of CuFeS2 have, for example, indicated that the magnetic moment of the Fe ion is substantially lower than that expected for a high spin Fe3<sup>+</sup> ion [20,21], an observation that is consistent with a model of CuFeS2 being a mixture of two distinct ionic states [Cu]+[Fe]3+[S]2<sup>−</sup> <sup>2</sup> and [Cu]2+[Fe]2+[S]2<sup>−</sup> <sup>2</sup> [22]. Experimental results regarding the oxidation states of Cu and Fe in CuFeS2 are inconclusive [23–27] and likely depend on the sample preparation. CuFeS2 is suspected by some to be a charge-transport-type insulator, with very small [28,29] or even possibly negative [30] charge-transfer energy, due to the particularly strong *pd* hybridisation between the Fe 3*d* and S 3*p* orbitals. Various authors have observed an optical resonance at approximately 500 nm in the absorbance spectrum of CuFeS2 nanomaterials and attributed it to the presence of an "intermediate band" (IB) between the VB and conduction band (CB) in CuFeS2 [31–37]. Several authors have demonstrated that this IB (also referred to as an "upper valence band" or "additional conduction band") leads to excellent photothermal conversion efficiency [31,32] and the potential for improved photovoltaic efficiency [38] or spintronic applications [39]. Gabka et al. demonstrated that the spectral position of the resonance peak in CuFeS2 could be red-shifted by increasing the [Cu]/[Fe] ratio (i.e., by decreasing Fe fraction) [40], and, more recently, Lee et al. demonstrated a stoichiometry-dependent transition between the dielectric resonance associated with the Fe 3*d* IB and a copper-vacancy-induced localised surface plasmon resonance (LSPR) in Fe-poor CuFeS2[36]. Gaspari et al. used a linear optical model to demonstrate that when an IB with fixed width is shifted towards the VB (i.e., reducing the VB–IB gap but increasing the IB–CB gap) the resonant absorption associated with the so-called Fröhlich condition, which for all-dielectric materials typically occurs in the deep UV, can be red-shifted into the visible region of the spectrum [32]. Yao et al. investigated the influence of the Fe and Cu oxidation states on the resonance feature and found that the incorporation of Fe2+ in CuFeS2 reduces the VB–IB gap, leading to a more prominent resonance feature [37]. In their work, they used the same model as Gaspari et al. but allowed the IB width to vary such that the VB–IB gap could be varied independently of the IB–CB gap. It was found that reducing the VB–IB gap in this manner blue-shifted the resonance feature in their simulated absorption spectra. These results illustrate that both the VB–IB gap and the IB width have a significant influence on the spectral position of the optical resonance.

Some work has been performed on transition-metal-doped CuAlS2 [29,41–43]. Two absorption bands, located at 1.3 eV and 2.0 eV, were identified by Teranishi et al. in the optical absorption spectrum of Fe:CuAlS2 (i.e., Fe-doped CuAlS2) [41] that increased in intensity with increasing Fe content. Sato et al. investigated CuAl0.9Fe0.1S2 and noted a small divalent component in the Cu 2*p* XPS peak in CuFeS2 (again, consistent with the model of Pauling) but not in CuAl0.9Fe0.1S2, suggesting hybridisation between the unoccupied Fe 3*d* orbitals with the valence Cu 3*d* orbitals mediated by the S 3*p* valence orbitals [29]. Wang et al. doped CuAlSe2 with various transition metals and concluded that Ti:CuAlSe2 is more promising for intermediate-band-based photovoltaics since it possesses a partially filled IB as opposed to the completely unoccupied IB in Fe:CuAlSe2 [44].

Recently, Yadav et al. demonstrated the ability to synthesise quaternary earth-abundant CuAl*x*Fe1−*<sup>x</sup>*S2 nanocrystals with a band gap that is tunable across the entire UV/vis/NIR spectrum, whilst maintaining the chalcopyrite crystal structure [45]. As part of this work, they performed a limited DFT investigation into the electronic structure of the CuAl*x*Fe1−*<sup>x</sup>*S2 system showing the projected density of states (PDOS) calculated for the CuAl0.75Fe0.25S2 system, the composition of which was obtained by substituting two Al3<sup>+</sup> ions with Fe3<sup>+</sup> ions, i.e., by introducing two [FeAl] <sup>0</sup> antisite defects. However, it is plausible that the site preference for Fe substitutions in CuAlS2 is on Cu sites rather than on Al sites. First, since the effective ionic radius of Fe3+ in tetrahedral coordination is larger than that of Al3+ and smaller than that of Cu<sup>+</sup> (see Table 1), a [FeAl] <sup>0</sup> substitution is likely to introduce more local distortion in the lattice than a [FeCu] <sup>2</sup><sup>+</sup> substitution. Second, since the Cu-S and Fe-S bonds are significantly more covalent than the Al-S bond, which is mostly ionic in nature [46], a [FeAl] <sup>0</sup> substitution would result in a more covalent local environment. Navrátil et al. proposed that intrinsic point defects such as [FeCu] <sup>2</sup><sup>+</sup> antisites play an integral role alongside the charge transport phenomenon in explaining the low effective magnetic moment and weak ferromagnetism [47].

**Table 1.** Effective ionic radius of Cu+, Fe3+ and Al3+ in tetrahedral coordination [48] (the Fe3+ ion is assumed to be in the high spin state).


If it were indeed the case that [FeCu] <sup>2</sup><sup>+</sup> antisites were more abundant than [FeAl] 0 antisites, the presence of interstitial Cui or [CuAl] <sup>2</sup><sup>−</sup> antisite defects would be necessary to maintain the stoichiometric [Cu]/([Al]+[Fe]) ratio. Following the same reasoning regarding the site preference for Fe outlined above, the formation of [CuAl] <sup>2</sup><sup>−</sup> antisite defects would likely introduce significant strain in the lattice. However, it may be that neutral [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect pairs have relatively low formation energies. It is well-known that the formation energies of intrinsic point defects, specifically cation vacancies and cation antisites, are comparatively low in ternary chalcogenides, such as CuInS2 and CuInSe2, with the chalcopyrite crystal structure, due to the lattice strain generated by the different cation–anion bond lengths. Furthermore, these charged defects can cluster to form neutral defect complexes such as [InCu] <sup>2</sup><sup>+</sup> + 2[VCu] <sup>−</sup> or [InCu] <sup>2</sup><sup>+</sup> + [CuIn] <sup>2</sup><sup>−</sup> [49–52]. Harvie et al. used electron energy-loss spectroscopy mapping to image the elemental distribution of copper and indium in single chalcopyrite CuInS2 QDs and found that copper-rich and indium-rich domains can exist within individual nanocrystals [53]. Another consideration then is that the additional configuration entropy associated with such antisite defect complexes when compared to the same number of isolated [FeAl] <sup>0</sup> antisite defects may be a significant stabilising factor, especially in CuAl1−*p*Fe*p*S2 nanocrystals such as those synthesised at elevated temperature by Yadav et al. This possibility has not been explored in detail here but may be investigated using the methods developed by Grau-Crespo et al. [54]. Overall, then, it is not obvious whether [FeAl] <sup>0</sup> antisite defects or [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect complexes would be the more prominent Fe substitution mode in CuAlS2.

We expand on the work undertaken by Yadav et al. by examining these two distinct modes of Fe substitution in CuAlS2. We first show that the site preference for Fe substitution is indeed on Cu sites as opposed to Al sites and that a [FeCu] <sup>2</sup><sup>+</sup> defect produces a local relaxation that reduces the formation energy of nearby [CuAl] <sup>2</sup><sup>−</sup> defects. We then show that the lowest unoccupied Fe 3*d* orbitals exist inside the band gap of CuAlS2 for both substitution modes, but that the [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> antisite complexes induce a charge transfer from the highest occupied Cu 3*d* orbital to the lowest unoccupied Fe 3*d* orbital, resulting in a [FeCu] <sup>+</sup> + [CuAl] − defect complex and a significantly reduced band gap. Finally, we calculate the per-defect formation energies and the PDOS of the CuAl1−*p*Fe*p*S2 system for each Fe substitution mode, for several values of *p* between zero and one obtained by introducing uniformly distributed defects. We find that the per-defect formation energy for [FeAl] <sup>0</sup> antisites increases linearly with the Fe fraction across the whole range of compositions studied (from CuAlS2 to CuFeS2), whereas the per-defect formation energy of the [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect complexes is reduced as the Fe fraction increases. We comment on this result in the context of some of the experimental results presented in the literature.

#### **2. Methods**

DFT calculations for this work were performed within the open-source CP2K framework [55]. The majority of data displayed in this paper was obtained using a generalised gradient approximation (GGA) specifically the PBE functional [56]. The Gaussian and Plane Wave method was used with Gaussian basis sets [57] of triple-*ζ* quality, a plane wave cutoff of 600 Ry for geometry relaxations (750 Ry for cell optimizations) and corresponding GTH pseudopotentials [58,59] which include 11, 16, 3 and 6 valence electrons for Cu, Fe, Al and S, respectively. Self Consistent Field (SCF) convergence criteria was used with the orbital transform algorithm [60] with the convergence set to EPS\_SCF = <sup>1</sup> × <sup>10</sup><sup>−</sup>7, giving energetic convergence of below a micro Hartree in total energy. All systems were fully relaxed with respect to both internal atomic coordinates and cell parameters. The defective systems were constructed from a 64-atom supercell, 2 × 2 × 1 of the default 16 atom chalcopyrite structure, thus giving 16 cation A sites (Cu for CuAlS2), 16 cation B sites (Al for CuAlS2) and 32 anion sites (S) allowing both flexibility in configurations and the possibility of reasonable sampling. Brillouin zone sampling was restricted to the Γ point.

For convenience, all systems were assumed to be fully ferro-magnetic: Whilst CuFeS2 should be antiferromagnetic, and, doubtless, the heavily defective systems will have a complex spin ordering, the effect on overall energetics and positioning of single particle levels will be small and of limited importance for this study. We use the minority (beta) spin state channel when reporting band gaps.

When calculating the relative energies of defective systems, the chemical potentials of individual species were calculated using the now standard ab initio thermodynamics approach and derived from total energy calculations of the unary, binary and ternary compounds. In the few cases where charged defects were considered, Lany–Zunger [61] charge corrections were applied.

We refer to Fe substitution via [FeAl] <sup>0</sup> antisites as "mode *α*" and to Fe substitution via [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect complexes as "mode *β*" (Figure 1). The configuration containing a single Fe substitution via mode *α* is referred to as '*α*1', and the configuration containing a single Fe substitution via mode *β* is referred to as '*β*1'. In the case of two Fe substitutions, for which we have studied several configurations, each configuration is labelled, for example, '*α*2.1', '*α*2.2', '*α*2.3', etc. (see Figure 1). When referring to a system with more than two Fe substitutions, in general we will use the *p* notation (as in CuAl1−*p*Fe*p*S2) or refer to specific configurations containing *n* defects per supercell as *α<sup>n</sup>* or *βn*. Since the 2×2×1 supercell contains 16 Al sites, the value of *p* for a configuration containing *n* defects is *p* = *n*/16. We used the *p* = 0.125 (*n* = 2) and *p* = 0.25 (*n* = 4) configurations to generate what we call the "inverted" *p* = 0.875 (*n* = 14) and *p* = 0.75 (*n* = 12) systems, respectively, where the "inverted" configuration contains *n* defects located at the non-defective sites of the (16 − *n*) configuration.

**Figure 1.** Illustrations showing the configurations *α*2.1, *α*2.2, *α*2.3, *β*2.1, *β*2.2, *β*2.3 and *β*2.4.

In order to determine the most uniformly distributed "*α*" configurations for a particular value of *p*, we utilised the ASE package [62] to generate all possible locations of *n* [FeAl] <sup>0</sup> defects. The average Fe-Fe distance in a given configuration was calculated using a minimum image convention (e.g., Fe*<sup>i</sup>* to Fe*<sup>j</sup>* would use the shortest of the distances between Fe*<sup>i</sup>* and all periodic images of Fe*j*). We find that the possible configurations split into multiple groups containing configurations with very similar mean Fe-Fe separation. From the group with the largest mean Fe-Fe separation, we selected the configuration with the lowest standard deviation about the mean as the most uniformly distributed configuration. To generate the *β* systems, the [FeAl] <sup>0</sup> were replaced with [CuAl] <sup>2</sup>−, and the nearest Cu site in the relaxed geometry was replaced by an Fe to form a [FeCu] <sup>2</sup><sup>+</sup> defect.

Owing to the apparent importance of on-site electronic interactions and correlations for understanding the properties of CuFeS2, DFT investigations based on LDA and GGA have struggled to provide reliable predictions of electronic structures of it and related materials. As a pragmatic approach, we have checked that the results are not qualitatively sensitive to the inclusion of exact exchange. After completion of PBE calculations, additional SCF calculations were performed using a generalised PBE0 type functional with 8, 16, or 24% exact exchange on the optimized PBE geometry.

## **3. Results and Discussion**

#### *3.1. Single (p* = 0.0625*) Fe Substitution*

We calculate the band gap of CuAlS2 to be approximately 1.65 eV (Figure 2), which is significantly smaller than the experimental value of 3.49 eV [63], as is typical for PBE [64]. The VB of CuAlS2 consists of mostly Cu 3*d* and S 3*p* character, implying pd hybridisation, consistent with the strong covalent character of the Cu-S bonds.

**Figure 2.** PDOS of CuAlS2 (copper = red, sulfur = yellow, aluminium = blue).

For the case of a single Fe substitution, we compared the formation energy of a [FeAl] <sup>0</sup> defect (configuration *α*1) with that of a single neutral [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect pair (configuration *β*1). Because of the symmetry of the CuAlS2 lattice, the site selection has a negligible effect on the calculated defect formation energy of the *α* substitution mode. In the case of the *β* substitution mode, the Cu site closest to the [CuAl] <sup>2</sup><sup>−</sup> defect was chosen for the Fe substitution. We also calculated the formation energy of single isolated [FeCu] 2+ and [CuAl] <sup>2</sup><sup>−</sup> defects. The calculated values for all defect types are shown in Table 2.

**Table 2.** Calculated defect formation energies in CuAlS2.


We find that the formation energy of a [FeCu] <sup>2</sup><sup>+</sup> antisite (0.471 eV) is significantly lower than that of a neutral [FeAl] <sup>0</sup> antisite (2.789 eV), as expected. To reiterate, we primarily attribute this to the difference in effective ionic radii and cation–anion bond lengths of the relevant ions in tetrahedral coordination. We also find that the combined formation energies of the isolated [FeCu] <sup>2</sup><sup>+</sup> and [CuAl] <sup>2</sup><sup>−</sup> defects (4.583 eV) is almost 1 eV larger than the formation energy of a neutral [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect pair (3.599 eV). We suggest, therefore, that when a relatively weak Cu-S (BDECu-S = 2.86 eV [65]) bond is replaced with a relatively strong Fe-S bond (BDEFe-S = 3.34 eV [65]), the surrounding Cu-S and Al-S bonds weaken in response and that it is this local weakening of the Al-S bonds that reduces the formation energy of [CuAl] <sup>2</sup><sup>−</sup> antisites close to the [FeCu] <sup>2</sup><sup>+</sup> defect. Indeed, when analysing the distribution of Al-S bond lengths in the system containing a single [FeCu] <sup>2</sup><sup>+</sup> antisite, we find that for Al sites far away from (i.e., more than 4 Å) the Fe substitution, the Al-S bond lengths at each site are uniform (between 2.26 Å and 2.27 Å with a standard deviation of less than 0.005 Å). However, close to (i.e., within 4 Å of) the Fe substitution, the four Al-S bond lengths at each Al site vary to a much greater degree. The mean values remain consistent between 2.26 Å and 2.27 Å, but the Al-S bond involving the sulfur that is also bonded to the substituted Fe weakens significantly (this Al-S bond length increases to more than 2.29 Å), and the other three strengthen slightly (see Figures A1–A3).

When a single Fe substitution is made, the atomic PDOS indicates that the fully occupied Fe 3*d* orbitals are quite broad and exist throughout the VB suggesting hybridisation with the S 3*p* orbitals (Figure 3), consistent with the strong covalency of the Fe-S bond in chalcopyrite [26,29]. The unoccupied Fe 3*d* orbitals exist inside the band gap of the CuAlS2 host and display some S character, suggesting weak pd hybridisation. In the case of configuration *α*1, the five unoccupied 3*d* orbitals resemble the canonical localised d-electron states Fe <sup>3</sup>*dz*<sup>2</sup> , Fe <sup>3</sup>*dx*<sup>2</sup>−*y*<sup>2</sup> , Fe <sup>3</sup>*dxy*, Fe <sup>3</sup>*dxz* and Fe <sup>3</sup>*dyz* (see Figure A4) and appear in two clearly separated regions within the band gap (Figure 3a), corresponding to the minority spin *e* and *t*<sup>2</sup> orbitals separated by a crystal field splitting of approximately 0.5 eV, which is smaller than that observed in Fe:ZnS [66], for example. We observe some internal structure in each of these 3*d* bands, which implies that the degeneracy of the two *e* orbitals (*dz*<sup>2</sup> and *dx*<sup>2</sup>−*y*<sup>2</sup> ) and three *<sup>t</sup>*<sup>2</sup> orbitals (*dxy*, *dxz* and *dyz*) has been lifted, presumably because the tetrahedral coordination is not ideal due to the varying cation–anion bond lengths.

In the case of configuration *β*1, the asymmetric crystal field is even more evident, with the five 3*d* orbitals appearing to be well-separated (Figure 3b)). This asymmetry is reflected in the Fe-S bond lengths: the four Fe-S bond lengths in the *α*<sup>1</sup> configuration are all approximately 2.297 Å (with a standard deviation of <0.001 Å), whereas the four Fe-S bond lengths in the *β*<sup>1</sup> configuration are 2.380 Å, 2.367 Å, 2.354 Å and 2.268 Å (corresponding to a standard deviation of 0.050 Å). By comparing the density of states of configuration *β*<sup>1</sup> (Figure 3b)) with that of configuration *<sup>α</sup>*1, we observe that the Fe 3*dx*<sup>2</sup>−*y*<sup>2</sup> orbital, which is the lowest unoccupied state in configuration *α*1, is the highest occupied state in configuration

*β*1, and the Cu 3*dxy* orbital, which is occupied in configuration *α*1, is shifted up in energy to about 100 meV above the VB edge. This is observed in the system containing only the isolated [FeCu] <sup>2</sup><sup>+</sup> antisite (see Figure A5), but the occupied Fe 3*d* and unoccupied Cu 3*d* are essentially degenerate (25 meV separation). We therefore attribute it to a charge transfer from Cu 3*dxy* to Fe 3*dx*<sup>2</sup>−*y*<sup>2</sup> , that is significantly enhanced by the presence of the [CuAl] <sup>2</sup><sup>−</sup> antisite and which results in a local [Cu]2+[Fe]2+[S]2<sup>−</sup> <sup>2</sup> ionic state. If we include the oxidation state explicitly in our defect notation, this corresponds to

$$[\mathrm{Fe}\_{\mathrm{Cu}^{\mathrm{I}}}^{\mathrm{III}}]^{2+} + [\mathrm{Cu}\_{\mathrm{Al}^{\mathrm{III}}}^{\mathrm{I}}]^{2-} \rightarrow [\mathrm{Fe}\_{\mathrm{Cu}^{\mathrm{I}}}^{\mathrm{II}}]^{+} + [\mathrm{Cu}\_{\mathrm{Al}^{\mathrm{III}}}^{\mathrm{II}}]^{-}. \tag{1}$$

Various studies have highlighted the importance of the Fe2+ content on the properties of CuFeS2 [37,47]. For example, Yao et al. have shown that the intensity of the plasmon-like feature in CuFeS2 grows with increasing Fe2+ content [37]. It was suggested by the authors that the Fe2+ state corresponds to an additional electron in the lowest unoccupied Fe 3*d* orbital, which effectively narrows the VB–IB gap, as we also suggest here. We therefore expect that this possible route to forming Fe2+ in Fe-substituted CuAlS2 may be relevant to understanding the nature of the optical resonance observed by Yadav et al. in the absorption spectrum of CuAl1−*p*Fe*p*S2 [45].

**Figure 3.** PDOS for configurations (**a**) *α*<sup>1</sup> and (**b**) *β*<sup>1</sup> (copper = red, sulfur = yellow, aluminium = blue and iron = green).

#### *3.2. Double (p* = 0.125*) Fe Substitution*

In the case of two Fe substitutions, we have investigated the effect of the precise configuration on the defect formation energies and PDOS. We find that the per-defect formation energy for the three configurations containing two [FeAl] <sup>0</sup> defects, labelled *α*2.1, *α*2.2 and *α*2.3 in Figure 1, is similar in all three cases (2.811 eV, 2.813 eV and 2.820 eV, respectively) and only slightly higher than the single defect case (2.789 eV). However, we do see some small variation in the PDOS (Figure 4), specifically that the variations in the local distortion already are sufficient for defect levels to span the range of energies associated with an IB such as that observed in CuFeS2.

**Figure 4.** PDOS for configurations (**a**) *α*2.1, (**b**) *α*2.2 and (**c**) *α*2.3 (copper = red, sulfur = yellow, aluminium = blue and iron = green).

In the case of [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect complexes we examined four configurations (*β*2.1, *β*2.2, *β*2.3 and *β*2.4 shown in Figure 1). We find that the per-defect formation energies in the *β*<sup>2</sup> configurations are more variable than for the *α*<sup>2</sup> configurations: for configurations *β*2.1, *β*2.2, *β*2.3 and *β*2.4, we calculated it to be 3.455 eV, 3.789 eV, 3.732 eV and 3.541 eV, respectively. We note that the two configurations with lower formation energies (*β*2.1 and *β*2.4) are those in which the [FeCu] <sup>2</sup><sup>+</sup> defects exist 'in between' the [CuAl] <sup>2</sup>−, which will become increasingly likely as the Fe fraction is increased via this Fe substitution mode. The PDOS for configurations *β*2.1, *β*2.2, *β*2.3 and *β*2.4 are shown in Figure 5. It is evident that the presence of the charge transfer is insensitive to small variations in the configuration.

#### *3.3. Multiple (p* > 0.125) *Fe Substitutions*

The calculated per-defect formation energies for both substitution modes across the entire composition space between CuAlS2 and CuFeS2 are shown in Figure 6. The per-defect formation energy of the [FeAl] <sup>0</sup> antisite increases linearly with the Fe fraction. In contrast, the formation energy per defect of the [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> complex decreases relatively quickly as *p* increases from 0 to 0.5, such that the two modes have almost indistinguishable defect formation energies for the ideal mixed system CuAl0.5Fe0.5S2.

**Figure 5.** PDOS for configurations (**a**) *β*2.1, (**b**) *β*2.2, (**c**) *β*2.3 and (**d**) *β*2.4 (copper = red, sulfur = yellow, aluminium = blue and iron = green).

In general, to generate the *β<sup>n</sup>* configurations from the *α<sup>n</sup>* configurations, we replaced the *n* [FeAl] <sup>0</sup> defects with [CuAl] <sup>2</sup><sup>−</sup> defects and chose to substitute the closest Cu to each [CuAl] <sup>2</sup><sup>−</sup> site for Fe to form *n* [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect pairs. We have investigated how this choice affects the formation energies at either very low (*p* < 0.2) or very high (*p* > 0.8) Fe content, where the configuration space is reasonably accessible. The four different values of the per-defect formation energy at *p* = 0.125 for the *β* configurations (red curve in Figure 6) correspond to the *β*2.1, *β*2.2, *β*2.3 and *β*2.4 systems (see Figure 1). It is clear that slight variations in the *β*<sup>2</sup> configurations can lead to significantly different defect formation energies. The reason for this fluctuation is likely the additional degree of freedom associated with the [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect pair relative to the isolated [FeAl] <sup>0</sup> defect. The order of the *β*<sup>2</sup> configurations in terms of defect formation energy is as follows: *β*2.2 > *β*2.3 > *β*2.4 > *β*2.1. It appears, therefore, that when the two [FeCu] <sup>2</sup><sup>+</sup> defects are located in between the two [CuAl] <sup>2</sup><sup>−</sup> defects (see Figure 1), the formation energy of the cluster as a whole is decreased. This is consistent with our reasoning in Section 3.1 regarding the lower formation energy of the [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect pairs relative to the sum of the formation energies of the isolated [FeCu] <sup>2</sup><sup>+</sup> and [CuAl] <sup>2</sup><sup>−</sup> defects and is also consistent with the observation of compositional domains in similar systems [53]. What is not clearly visible in Figure 6 is that we have also performed similar calculations for the three *α*<sup>2</sup> configurations *α*2.1, *α*2.2 and *α*2.3 (blue curve at *p* = 0.125), but the variation in per-defect formation energies is almost imperceptible. As the Fe fraction increases, we would expect this additional configurational freedom in the *βn*><sup>1</sup> systems and, therefore, the variation in the calculated defect formation energies, to also decrease. Indeed, we also calculated the variation in the per-defect formation energy of the four *β*<sup>14</sup> systems (at *p* = 0.875) generated by 'inverting' the *β*<sup>2</sup> systems and found that the fluctuation is significantly smaller than at *p* = 0.125. However, the full exploration of this hypothesis is beyond the scope of this study.

**Figure 6.** Defect formation energies (per-defect) for[FeAl] <sup>0</sup> antisites (blue circles) and [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] 2− defect complexes (red crosses). The multiple values at *p* = 0.125 and *p* = 0.875 correspond to the *α*2.1, *α*2.2, *α*2.3, *β*2.1, *β*2.2, *β*2.3 and *β*2.4 configurations and their inverted systems, respectively.

To explain the observed trends in the per-defect formation energies of the two Fe substitution modes, we examine the composition-induced structural changes in both systems. [AI ][BIII][CVI]2 compounds in the chalcopyrite crystal structure are derived from the zinc-blende structure by replacing the group II element of a binary chalcogenide compound with group I and group III elements in an ordered manner such that the average number of valence electrons per atom is four (i.e., according to the extension of the Grimm–Sommerfeld rules to ternary compounds by Goryunova [67]). In such ternary systems, where the cationic sub-lattice contains two different cation species, the two cation–anion bond lengths can differ (*dA*−*<sup>C</sup>* = *dB*−*C*), potentially causing significant strain in the crystal lattice. The tetragonal distortion parameter *η* = *c*/2*a*, where *c* is the lattice parameter in the *z* direction and *a* is the lattice parameter in the *x* direction, characterises the compression (*η* < 1) or tension (*η* > 1) along the *z* direction. Figure 7a) shows the tetragonal distortion parameter *η* for the range of quaternary compositions between CuAlS2 and CuFeS2. It can be seen that *η* approaches a value of approximately 0.997 when *p* ≥ 0.75, which is significantly closer to the ideal chalcopyrite structure (*η* = 1) than the value of 0.983 observed when *p* = 0. This indicates an expansion of the primitive cell along the *z* direction, which is consistent with the observed increase in primitive cell volume of about 1% (Figure 7b)).

**Figure 7.** (**a**) Tetragonal distortion parameter *η* for both modes of Fe-substitution over the range 0 < *p* < 1; and (**b**) Unit cell volume for both modes of Fe-substitution over the range 0 < *p* < 1 (blue circles are *α* configurations and red crosses are *β* configurations. The multiple values at *p* = 0.125 and *p* = 0.875 correspond to the *α*2.1, *α*2.2, *α*2.3, *β*2.1, *β*2.2, *β*2.3 and *β*2.4 configurations and their inverted systems, respectively.

In CuAl1−*p*Fe*p*S2, the variation in the cation–anion bond lengths increases substantially as *p* increases. As can be seen in Figure A6, the variation about the mean Fe-S, Al-S and Cu-S bond lengths is greater for systems containing the defect complexes than it is for systems containing the isolated antisite defects, reflecting the additional degree of freedom associated with the precise [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> arrangement. As we mentioned in the introduction, the large configuration entropy in materials with multiple components in a disordered configuration can be a significant stabilising factor [54,68]. We therefore anticipate that the [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> complexes become even more competitive with increasing Fe fraction than is suggested by the data presented in Figure 7.

When the Fe substitutions occur in the form of a [FeAl] <sup>0</sup> antisite defect, *d<sup>p</sup>* Fe-S increases gradually with the Fe fraction, only broadening slightly for intermediate values of *p* (Figure A6c)). When analysing the distributions of Cu-S and Al-S bond lengths in the *α* configurations, we see that as soon as a single [FeAl] <sup>0</sup> defect is introduced into the supercell (*p* = 0.0625), the *d<sup>p</sup>* Cu-S values appear to split into two clusters. We observe a similar splitting, albeit less prominent, in the *d<sup>p</sup>* Al-S values for 0 < *p* < 1 (Figure A6b)), which may be rationalised as follows. The secondary cluster of shorter bond lengths is likely to be made up of those Al-S bonds close to the [FeAl] <sup>0</sup> antisite, which are strengthened due to the replacement of the Al-S bond, which has a bond dissociation energy BDEAl-S = 3.88 eV [65], with a slightly weaker bond (BDEFe-S = 3.34 eV [65]). This local strengthening of the Al-S bonds close to the existing [FeAl] <sup>0</sup> antisites hinders the formation of additional [FeAl] <sup>0</sup> antisites.

In the case of Fe substitutions via the introduction of [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect complexes, we again see that as soon as a single defect is introduced, the *d<sup>p</sup>* Cu-S values appear to split into two distinct clusters (Figure A6d)). Since the bond dissociation energy of Cu-S is significantly lower than that of Fe-S, the [FeCu] <sup>2</sup><sup>+</sup> defects substitute a relatively weak bond with a relatively strong bond, and the surrounding Cu-S and Al-S bonds weaken. It is likely that this local weakening of the Al-S bonds reduces the formation energy of

[CuAl] <sup>2</sup><sup>−</sup> antisites close to existing [FeCu] <sup>2</sup><sup>+</sup> defects. This is consistent with our earlier observation (Section 3.1) that the combined formation energies of the isolated [FeCu] <sup>2</sup><sup>+</sup> and [CuAl] <sup>2</sup><sup>−</sup> defects is larger than that of the neutral [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect complex. It is also possible that the local weakening of Cu-S bonds reduces the formation energy for additional [FeCu] <sup>2</sup><sup>+</sup> antistites. This would be consistent with the observation (Figure 6) that the per-defect formation energy of a [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect complex initially decreases with increasing defect density.

Figure 8a,b show the PDOS of CuAl0.75Fe0.25S2 obtained via uniformly distributed [FeAl] <sup>0</sup> antisite defects (configuration *α*4) and [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect complexes (configuration *β*4), respectively. We still see quite a large difference between the *α* and *β* configurations. In the *β*<sup>4</sup> configuration, the highest occupied minority spin states have significantly more Fe 3*d* character than in the *α*<sup>4</sup> configuration, and the lowest unoccupied minority spin states have significantly more Cu 3*d* character.

**Figure 8.** PDOS for (**a**) *α*<sup>4</sup> and (**b**) *β*<sup>4</sup> systems (i.e., *p* = 0.25) (copper = red, sulfur = yellow, aluminium = blue and iron = green).

Figure 9a,b show the PDOS of CuAl0.25Fe0.75S2 obtained via uniformly distributed [FeAl] <sup>0</sup> antisite defects (configuration *α*12) and [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect complexes (configuration *β*12), respectively. It is clear that once the Fe content is this high, both systems begin to resemble CuFeS2. However, there is still a major difference between the *α*<sup>12</sup> and *β*<sup>12</sup> configurations, in that the latter has a significantly narrower VB–IB gap.

Figure 10 shows the calculated band gap of CuAl1−*p*Fe*p*S2 obtained via uniformly distributed [FeAl] <sup>0</sup> antisite defects and [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect complexes. In their 2021 article, Yadav et al. calculated a band gap of 1.78 eV for the CuAlS2 system, a band gap of 0.21 eV for the CuFeS2 system and a band gap of 0.38 eV for the CuAl0.75Fe0.25S2 system. It can be seen from Figure 10 that our calculations are in approximate agreement. We see that a single Fe substitution in the form of a [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect pair reduces the band gap from 1.6 eV to approximately 0.1 eV, whereas for a single [FeAl] <sup>0</sup> defect the band gap is reduced to 0.6 eV. By the time *p* = 0.5, the band gap of both the *α* and*β* configurations is approximately 0.2 eV. When *p* > 0.5, we see that the band gap of the *α* configuration increases again to about 0.4 eV for *p* = 0.75 before dropping to about 0.2 eV for *p* = 1.0, whereas the band gap of the *β* configuration remains at about 0.2 eV.

**Figure 9.** PDOS (**a**) *α*<sup>12</sup> and (**b**) *β*<sup>12</sup> systems (i.e., *p* = 0.75) (copper = red, sulfur = yellow, aluminium = blue and iron = green).

**Figure 10.** Band gap of the CuAl1−*p*Fe*p*S2 obtained via uniformly distributed [FeAl] <sup>0</sup> antisite defects (blue circles) and [FeCu] <sup>2</sup><sup>+</sup> + [FeCu] <sup>2</sup><sup>+</sup> defect complexes (red crosses). The multiple values at *p* = 0.125 and *p* = 0.875 correspond to the *α*2.1, *α*2.2, *α*2.3, *β*2.1, *β*2.2, *β*2.3 and *β*2.4 configurations and their inverted systems, respectively.

As was mentioned previously, Yao et al. observed a plasmon-like resonance at 500 nm in the absorbance spectrum of CuFeS2 and attributed it to the presence of Fe2<sup>+</sup> [37]. We notice that in the data presented by Yadav et al. a similar resonance, also centred at 500 nm, appears in the absorbance spectrum of CuAl*x*Fe1−*<sup>x</sup>*S2 nanocrystals only when *x* < 0.5 (in our notation this corresponds to *p* > 0.5 in CuAl1−*p*Fe*p*S2) rather than increasing gradually in amplitude as soon as Fe is introduced [45]. We therefore speculate that for systems with *p* > 0.5 in which the [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> complex becomes relatively abundant, the charge transfer that we observe for the *β* configurations may be responsible for the experimentally observed plasmon-like absorption feature. This also is consistent with the fact that Sato et al. identified a small divalent component in the Cu 2*p* XPS peak in CuFeS2 but not in CuAl0.9Fe0.1S2 [29]. We therefore expect that a comprehensive XPS study of CuAl1−*p*Fe*p*S2 nanocrystals will enable the determination of the precise Fe fraction at which the Cu2<sup>+</sup> component emerges. We also expect that a detailed XRD study, similar to the one performed by Yao et al., who demonstrated the site preference for Mn on the Cu site in CuInSe2 by performing Rietveld refinements of X-ray diffraction data [69], will shed light on the prevalence of the [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect complexes we propose here.

#### **4. Conclusions**

In summary, we have used density functional theory to investigate Fe substitution in CuAlS2 via two modes, namely [FeAl] <sup>0</sup> antisites and [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect complexes. We find that the formation energy of [FeCu] <sup>2</sup><sup>+</sup> antisites is significantly lower than that of [FeAl] <sup>0</sup> antisites. Furthermore, we find that the presence of [FeCu] <sup>2</sup><sup>+</sup> promotes the formation of nearby [CuAl] <sup>2</sup><sup>−</sup> antisites and that neutral [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> defect complexes are formed competitively when *p* ≥ 0.5. Analysis of electron density and density of states reveals that charge transfer within these defect complexes leads to the formation of local Cu2+/Fe2+ ionic states, whereas the dominant ionic state in the [FeAl] <sup>0</sup> substituted systems is Cu+/Fe3+. We speculate that charge transfer processes such as the one discussed in this article can lead to the formation of broad IBs with narrow VB–IB gaps in similar transition metal-containing quaternary and penternary chalcogenides. A clearer understanding of these processes and how they can be utilised to tune the electronic structure of materials in order to generate desirable optoelectronic properties may pave the way for their application in photovoltaics and spintronics.

**Author Contributions:** Conceptualization, M.B. and C.D.; data curation, C.D.; formal analysis, C.D. and A.O.J.K.; funding acquisition, M.W. and M.B.; investigation, C.D. and A.O.J.K.; methodology, M.W. and C.D.; project administration, M.B.; resources, M.W.; software, M.W.; supervision, M.B. and M.W.; validation, M.W.; visualization, C.D. and A.O.J.K.; writing—original draft, C.D. and M.B.; writing—review and editing, M.B., C.D. and M.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported via our membership of the UK's HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/R029431), this work used the ARCHER2 UK National Supercomputing Service and the UK Materials and Molecular Modelling Hub for computational resources, MMM Hub, which is partially funded by EPSRC (EP/T022213) for YOUNG.

**Data Availability Statement:** Data available on request due to storage size of data set.

**Acknowledgments:** M.B. would like to acknowledge Richard Robinson for clarifying the band structure used to generate their simulated absorption spectra in reference [37]. A.K. would like to thank the University of Lincoln Undergraduate Research Opportunities Scheme for funding his summer research project.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


**Figure A1.** Al-S bond lengths for Al close to [FeCu] <sup>2</sup><sup>+</sup> antisite defect (only half of the nearest neighbour Al sites are shown for purposes of clarity).


**Figure A2.** Al-Fe separation and Al-S bond lengths for all Al sites in CuAlS2 containing an isolated [FeCu] <sup>2</sup><sup>+</sup> antisite defect (blue cells indicate a larger value relative to the mean and red cells indicate a smaller value relative to the mean).

**Figure A3.** Al-S bond lengths plotted against Al-Fe separation for CuAlS2 containing an isolated [FeCu] <sup>2</sup><sup>+</sup> antisite defect.

**Figure A4.** 3*d* orbitals in configuration *α*<sup>1</sup> (top row) and configuration *β*<sup>1</sup> (bottom row). Yellow circles = sulfur, red circles = copper, blue circles = aluminium and green circles = iron.

**Figure A5.** PDOS of CuAlS2 containing an isolated [FeCu] <sup>2</sup><sup>+</sup> antisite defect (copper = red, sulfur = yellow, aluminium = blue and iron = green).

**Figure A6.** (**a**) Cu-S bond lengths in [FeAl] <sup>0</sup> substituted CuAl1−*p*Fe*p*S2; (**b**) Al-S bond lengths in [FeAl] 0substituted CuAl1−*p*Fe*p*S2; (**c**) Fe-S bond lengths in [FeAl] <sup>0</sup> substituted CuAl1−*p*Fe*p*S2; (**d**) Cu-S bond lengths in [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> substituted CuAl1−*p*Fe*p*S2; (**e**) Al-S bond lengths in [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> substituted CuAl1−*p*Fe*p*S2; (**f**) Fe-S bond lengths in [FeCu] <sup>2</sup><sup>+</sup> + [CuAl] <sup>2</sup><sup>−</sup> substituted CuAl1−*p*Fe*p*S2.

**Figure A7.** Band gaps for the *α* (blue circles) and *β* (red crosses) configurations in (**a**) major; and (**b**) minor spin channels.
