Functional off-stoichiometry in Cu(In,Ga)Se₂. Part I: topotactic continuum of ordered defect compounds

Abstract

There is hardly a more bewildering aspect of the chalcopyrite material Cu(In,Ga)Se₂ (CIGSe), used as an absorber in thin-film solar cells, than its relation with off-stoichiometry, which is often manifested in the formation of ordered defect compounds (ODCs). Despite huge efforts devoted to determining the structure and properties of the ODCs, no comprehensive model has been provided. The convention is that the off-stoichiometry is mediated by the charge-neutral [2V_Cu + III_Cu] complexes forming in high concentrations. In this work, we refute this simplistic picture by treating the off-stoichiometry as an integral part of the crystal structure. In part I, herein, we report a discovery of the topotactic transformation series of stable ODC-like structures with different compositions that enable large Cu deficiency in the material. Crucially, the underlying structures are shown practically invisible to X-ray and neutron diffraction, as the parent lattice is preserved and barely deformed within the relevant compositional range. The model we deduce reconciles several lasting conflicts between experimental findings and first-principles calculations. In part II, this model is further validated by computing the electronic properties for the structures in the series and comparing the results with prior experimental findings. The developed structural series is uniquely self-consistent and is thus well-positioned to bolster computational analysis of chalcopyrite compounds and alloys.

---

1. Introduction

High tolerance to off-stoichiometry is a defining feature of Cu(In,Ga)Se₂ (CIGSe) solar absorbers that enabled the continuous development of thin-film solar cells for decades. Recently, this development culminated in a new power conversion efficiency record of 23.6% achieved with group-I-poor Ag-alloyed CIGSe.¹ The group-I deficiency in the record device is not a surprise considering that making the absorber overall Cu-poor is the standard fabrication strategy aimed at mitigating the detrimental segregation of Cu_xSe and inducing benign Cu depletion at grain boundaries.² The grain boundaries thus account for a fraction of the deviation from the nominal [Cu]/[In] = 1, but their small volume fraction implies the dominant role of the grains. The phenomenological ability of the absorbers to accommodate Cu deficiency is reflected in the wide homogeneity region of α-CIGSe in the equilibrium phase diagrams (see Fig. 1a and b). The added off-stoichiometry entails certain changes of the bulk properties, which at present cannot be properly assessed because the available structural models are dubious. Specifically, despite the broad homogeneity region of α-CIGSe, it is typically simulated as a pristine stoichiometric chalcopyrite compound or, when the role of Cu-deficiency is in focus, as chalcopyrite with [2V_Cu + III_Cu] complexes or isolated copper vacancies.^3–7 Such a dilute limit description is inapplicable because the average distance between complexes (assuming equal spacing on a three-dimensional grid) at [Cu]/([Ga] + [In]) = 0.9 is under 1.5 nm – sufficiently close to trigger complex clustering of some sort. This recognition could have been behind the idea of charge-neutral off-stoichiometric superclusters, which can form to accommodate off-stoichiometry without drastically changing the carrier concentrations or soaring scattering in α-CIGSe.⁸ However, the existence and structure of such superclusters remain obscure.

Fig. 1 Experimental phase diagram of the pseudo-binary (a) Cu₂Se–In₂Se₃ system adapted from ref. 9 and (b) Cu₂Se–Ga₂Se₃ system adapted from ref. 10. (c) The most important crystal structures reported for the α- and β-phases. The 1 : 3 : 5 ODC structural models were resolved using the neutron and X-ray diffraction experiments in ref. 11. The 1 : 5 : 8 ODC structural models were proposed in ref. 12. A more extensive list of literature structures is given in Table S1.

Structural models of β-CIGSe, which is a phase that forms alongside α-CIGSe at large Cu deficiencies, are also limited. Early studies from the 1960s–80s revealed that the β-CIGSe phase consists of a series of ordered defect compounds (ODCs), including Cu₃In₅Se₉, Cu₂In₄Se₇, Cu₈In₁₈Se₃₂, Cu₇In₁₉Se₃₂, CuIn₃Se₅, and CuIn₅Se₈ (referred to as 3 : 5 : 9, 2 : 4 : 7, etc.), and therefore spans a wide composition range (see Fig. 1a and b). However, this range is often reduced to only two compounds, 1 : 3 : 5 and 1 : 5 : 8, in solar cell studies. Experimentally, several crystal structure models were derived based on Rietveld refinements, where different fractional occupancies for cationic sites were proposed, as illustrated for the 1 : 3 : 5 ODC in Fig. 1c.¹¹ Computationally, a structural motif based on agglomeration of defect complexes was first discussed by Zhang et al.^13,14 Currently, the most stable structure of the 1 : 5 : 8 ODC proposed therein is confirmed stable relative to a mixture of 1 : 1 : 2 and 0 : 2 : 3,^13–15 and it is shown effective in reproducing the characteristic diffraction ODC peaks, motivating widespread adoption.^15–20 Later, Chang et al.¹² presented several alternatives for the 1 : 5 : 8 ODC structure, labelling the initial model by Zhang et al. as “Type-D”. These alternatives are now known to be unstable,²¹ but they still appear in computational studies (e.g. “Type-A” in ref. 22 and 23 and “Type-E” in ref. 24–28). The three most investigated 1 : 5 : 8 ODC models from Chang et al.¹² are illustrated in Fig. 1c.

Besides the 1 : 5 : 8 composition, very few ODC models compatible with computational analysis (i.e. without fractional site occupancies) have been suggested. The only model of note is the zinc-blende-derived 1 : 3 : 5 ODC studied in ref. 22, 23 and 29. As such, some authors resort to less rigorous ways to reproduce the off-stoichiometry, i.e. (i) treating most ODCs as weighted averages of 1 : 1 : 2 and 1 : 5 : 8 structures,^14,20,30 (ii) saturating the lattice with defects or complexes,^{3,18,28,31,32} (iii) integrating 1 : 1 : 2 and 1 : 5 : 8 domains quasi-randomly into large supercells,²⁴ etc. However useful, these approximations disregard stability concerns and thus validity of the results is questionable.

In part I of this contribution, we report a computational discovery of a topotactic series of stable (or near-stable) ODC structures in pseudo-binary Cu₂Se–In₂Se₃ and Cu₂Se–Ga₂Se₃ systems. We show that these series are underpinned by energetically favorable extended vacancy arrangements forming in different quantities and constellations. We describe the common structural motifs within the series, simulate diffraction patterns for the individual ODCs, and review prior experimental observation through the lenses of these new findings. We explain the wide homogeneity regions of both α-CIGSe and β-CIGSe by the continuous character of structural changes with off-stoichiometry. We conclude by outlining the implication of the generated knowledge for controlling phase composition and defect engineering, both of which are indispensable for the further development of photovoltaic devices. In part II, we investigate electronic properties of (near-)stable structures in the series and further validate the structural models via comparison with available experimental evidence.³³

2. Methods

2.1 Computational parameters

All first-principles calculations presented in this contribution were performed using the Vienna Ab initio Simulation Package (VASP)³⁴ employing the projector augmented wave (PAW)³⁵ method within density functional theory (DFT). Depending on the property of interest, the exchange-correlation functional was toggled between Perdew–Burke–Ernzerhof (PBE),³⁶ revised PBE for solids (PBEsol),³⁷ PBE with Hubbard U correction of 5 eV applied on the Cu d orbitals according to Dudarev et al. (PBE+U),³⁸ modified Becke–Johnson (mBJ),³⁹ strongly constrained and appropriately normed meta-GGA (SCAN),⁴⁰ where GGA stands for generalized gradient approximation, and hybrid Heyd–Scuseria–Ernzerhof with a 25% fraction of exact exchange (HSE06).⁴¹ Pseudopotentials with the following valence electron configurations were selected: Cu 3d¹⁰4s¹, Ga 4s²4p¹, In 5s²5p¹, and Se 4s²4p⁴. The Brillouin zones were represented by automatically generated Γ-centered uniform k-point grids of different densities specified below. Gaussian smearing with a 10 meV distribution width was chosen to define partial occupancies and enable batch processing of structures with finite and zero band gaps predicted using the PBE, PBEsol, PBE+U, and/or SCAN functionals. All calculations were performed in the non-spin-polarized regime, reasoning that the materials of interest are typically non-magnetic.

2.2 Data processing

Data processing was facilitated by pymatgen (Python Materials Genomics)⁴² and the structures were visualized using the Visualization for Electronic and STructural Analysis (VESTA).⁴³

2.3 High-throughput screening

Our strategy to identify stable ODC structures is essentially a brute-force screening of randomly generated structures, with a few tricks accelerating the search. First, a small initial dataset with about 450 Cu–In–Se structures containing up to 32 sites was generated using a custom-made script relying on the octet rule and oxidation state balance conditions (see below). This initial dataset was complemented by known literature structures. Second, the total energies of all those structures were calculated within DFT, and respective convex hulls were constructed for the ground-state stability analysis. Third, all structures within 10 meV per atom from the convex hull were selected to produce cluster expansion decomposition for total energies utilizing the Alloy Theoretic Automated Toolkit (ATAT).⁴⁴ The cluster expansion enabled quick pre-assessment of total energies, which was incorporated into the structure generation script as an additional filter to remove highly unstable structures on the fly. Fourth, the updated script was executed to generate more structures, ground-state stability of which was subsequently analyzed at the DFT level, and the resulting energies were used to update the convex hull and the cluster expansion. The generation/calculation/update loop was repeated several times until all possible Cu–In–Se structures in supercells containing up to 48 atomic sites (including vacancies) were considered. Then, hundreds more structures at three special compositions (i.e. 1 : 1 : 2, 1 : 3 : 5, and 1 : 5 : 8) with unit cells containing up to 64 sites were generated and analyzed. Considering the large dataset size, all DFT calculations at the screening stage were done using relatively coarse accuracy parameters (see below for specifics). In total, over 100 000 non-equivalent Cu–In–Se structures were considered, of which over 3000 were investigated at the DFT level. No full-fledged screening was attempted for the Cu–Ga–Se system, where the analysis was limited to isomorphs of relatively stable Cu–In–Se structures. The main structures discussed in this work can be accessed using the following identifier.⁴⁵

2.4 Structure generation script

The structure generation script was written in Python, relying heavily on the pymatgen library.⁴² Its basic function was to fill all cationic sites in the ideal zinc-blende supercells of various sizes with either vacancies or metal elements in different ratios, keeping the anionic sublattice filled by Se. The cationic ratios were defined such that the obtained structures fall on the pseudo-binary Cu₂Se–In₂Se₃ line or, in other words, assuming overall charge neutrality with the Cu, In, and Se ions having 2+, 4+, and 2− charged states (i.e. the expression 2[Cu] + 4[In] − 2[Se] = 0 was enforced). This condition is referred to as the oxidation state balance filter herein. Another filter, based on the minimal deviation from the octet rule (described in the next section), was incorporated to perform on-the-fly exclusion of all highly unstable structures. The third screening filter based on the cluster expansion was added at the later stages for more accurate energy pre-assessed during structure generation, as described above.

2.5 Minimal deviation from the octet rule

The octet rule is an empirical prerequisite of stability for chalcogenide structures with zinc-blende-derived lattices.^13,46,47 It states that the low-energy configurations have the sum of oxidation states of the four cations surrounding the anions equal to eight. At the 1 : 1 : 2 composition in the Cu–In–Se system, this condition translates into the requirement for all Se atoms to be surrounded by two Cu atoms and two In atoms, comprising a κ = 8 tetrahedron (where κ is the sum of cation oxidation states, see Fig. 2). Both the ground-state chalcopyrite and low-energy CuAu-type 1 : 1 : 2 structures in the I–III–VI systems are compliant with this rule. For the Cu-poor ODC structures, however, κ = 8 tetrahedra are insufficient and other building blocks must be introduced. The empirical stability condition for the ODCs was therefore updated to the minimal deviation from the octet rule, achieved by introducing tetrahedra with κ = 7 (two In, one Cu, and one cation vacancy) and κ = 9 (three In and one cation vacancy). The set of three tetrahedral types suffices to construct zinc-blende-derived structures of any composition between 1 : 1 : 2 and 1 : 5 : 8 inclusively, but further reduction in the Cu content requires κ = 6 (two In and two cation vacancies). In this work, the only low-energy structures with κ = 6 are the terminal In₂Se₃ and Ga₂Se₃ (denoted 0 : 2 : 3), but other stable ODCs with intermediate compositions and larger supercells may be identified in the future.

Fig. 2 Fraction of different building blocks (tetrahedra around Se ions) in the Cu–III–Se structures corresponding to the minimal deviation from the octet rule. The “vacancy content” is the fraction of lattice sites filled with vacancies. The tetrahedra are denoted by the sum of valence electrons for the four cations surrounding Se ions (κ values). The lines show linear dependencies and the markers represent the values for the topotactic series. The secondary on-top x-axes give alternative scales: “[Cu]” refers to the atomic fraction of copper (in at%) and “fraction of III₂Se₃” is defined for a mixture with Cu₂Se (equal to [III]/([Cu] + [III]) expressed in percent).

2.6 Coarse stability calculations

The initial identification of stable structures was performed by high-throughput screening, which required a reasonable compromise between computational speed and accuracy. Such a compromise was achieved by choosing the PBE functional, a k-point grid density of about 1500 points per reciprocal atom, a cut-off energy of 400 eV, and a force threshold of 10 meV Å⁻¹ for geometry optimization. All cell parameters were allowed to relax during geometry optimization.

2.7 Fine stability calculations

For eleven (near-)stable ODCs from within 0.1 meV per atom above the convex hull constructed based on our screening results, as well as for several special ODCs and literature structures, careful calculations of the formation energies with five different functionals were performed. For the calculations using the PBE, PBEsol, PBE+U, and SCAN functionals, the k-point density, the cut-off energy, and the force threshold for ionic relaxations were increased to about 3000 points per reciprocal atom, 550 eV, and 5 meV Å⁻¹, respectively. For the HSE06 calculations, the corresponding parameters were set to about 1500 points per reciprocal atom, 350 eV, and 50 meV Å⁻¹ because this functional is more computationally demanding.

2.8 Literature structures

A smaller set of 40 Cu–In–Se structures, not necessarily with a zinc-blende-derived lattice, was compiled from previous publications,^{12–28,48–58} the Materials Project repository,²⁹ the Open Quantum Materials Database (OQMD),⁵⁹ and the Inorganic Crystal Structure Database (ICSD),⁶⁰ by locating the structures reported for the I–III–VI systems (I = Cu and Ag; III = In, Ga, and Al; VI = S, Se, and Te) and performing isovalent substitutions. We intended to include all known I–III–VI structures reported anywhere, but in practice, some highly unstable structures or large quasi-random supercells were deemed irrelevant and thus ignored. The collected structures were added to the convex hull at the early stages of high-throughput screening and structure generation. The most notable of these are illustrated in Tables S1 and S2, SI.

2.9 Defect complex stability

To evaluate the relative stability of the [2V_Cu + III_Cu] complexes, additional PBE calculations were performed for the complexes in a tetragonal supercell containing 216 atoms. The computational setup from the coarse stability calculations was adapted for all parameters except the k-point grid, which was changed to 2 × 2 × 2 in this case. Two possible complex geometries were considered, but only the most stable one is referenced below as the computed energies were found to differ by only 20 meV per complex in both the Cu–In–Se and Cu–Ga–Se systems.

2.10 Diffraction patterns and geometries

Unless specified, all geometries described herein correspond to the PBEsol results because they are found to most accurately reproduce the experimental lattice (see below). The PBEsol-optimized structures are, in particular, used for the diffraction pattern simulations, which were done using a diffraction module in pymatgen⁴² for an X-ray wavelength of 1.542 Å (Cu Kα line) or a neutron wavelength of 1.797 Å. To emulate the spectral line broadening, Lorentzian smearing with a full width at half maximum (FWHM) of 0.3° was applied. For the analysis of different polytypes in the SI, as a shortcut, X-ray diffraction (XRD) patterns were simulated for PBE-optimized geometries with the lattice volume rescaled using a certain coefficient such that the volume of the most stable polytype is equal to the PBEsol-optimized value.

2.11 CA-type fraction

To describe the average local atomic neighborhood of the cations, a parameter entitled “CA-type fraction” was defined as the fraction of Cu atoms with local coordination resembling that in the CuAu-type CuInSe₂ and CuGaSe₂ structures. The resemblance was evaluated by matching the atomic clusters within the second coordination sphere of Cu ions (to the nearest group-III cations) with several standard clusters derived from the CuAu-type structure through modifications illustrated in Fig. S1. A similar parameter was already used in our previous study of CuInS₂ polytypes.⁶¹ The bottom line here was that the clusters counted as CA-type are closer to those in the CuAu-type than those in the chalcopyrite-type structures. A complementary parameter that quantifies the chalcopyrite-type coordination, entitled “CH-type fraction”, is defined analogously.

2.12 Vacancy arrangement types

To facilitate the analysis of vacancy arrangements in all generated ODCs, similar to the CA-type fraction, a set of clusters was derived from several stable ODC structures (see Fig. S2) and matched with the vacancy distributions in studied structures. With this procedure, a distinction was made between vacancies arranged into (001)_CH planes (where the subscript “CH” refers to the chalcopyrite conventional unit cell) and those comprising arrays with straight or zigzag connectivity. An automatic sampling was subsequently performed to group all ODC structures into distinct categories, as discussed below.

3. Results

3.1 Convex hull analysis

To determine the ground-state stability of all generated structures and thereby establish their thermodynamic feasibility, we first construct pseudo-binary convex hulls for both Cu₂Se–In₂Se₃ and Cu₂Se–Ga₂Se₃ systems using PBE-computed enthalpies. The convex hulls produced in this way are shown in Fig. 3a and b. For convenience, we herein present all enthalpies relative to the mixture of 1 : 1 : 2 (i.e. CuInSe₂ or CuGaSe₂) and 0 : 2 : 3 (i.e. In₂Se₃ or Ga₂Se₃). The red markers in convex hulls represent literature structures and the green/blue markers depict stable/unstable structures found by the screening. Since convex hulls of both systems are qualitatively similar, the behavior of the quaternary CIGSe alloy is expected to be analogous to that of the ternary systems.

Fig. 3 Stability analysis of zinc-blende-derived CIGSe structures. (a and b) Convex hulls based on the formation enthalpies (ΔH_f) and (c and d) decomposition enthalpies (ΔH_d) for the structures in the (a and c) Cu–In–Se and (b and d) Cu–Ga–Se systems. The shaded regions outline the homogeneity ranges in the experimental phase diagrams (see Fig. 1) for α (red), β (blue), and γ (green) phases at 700 °C for the Cu–In–Se system⁹ and at 800 °C for the Cu–Ga–Se system.¹⁰ All datapoints shown in this figure are from the PBE calculations, while the convex hulls obtained with other functionals are given in Fig. S3.

Fig. 3 reveals the existence of eleven ODCs with enthalpies on the convex hull (stable structures) or within 0.1 meV per atom above it (near-stable structures) in both ternary systems, along with the terminal phases from the literature (i.e. 1 : 1 : 2 and 0 : 2 : 3). This means that all ODCs in question have lower enthalpies compared with all other polytypes and all possible mixtures of other ODCs that yield the same composition. The small 0.1 meV per atom above-hull threshold represents the estimated tolerance for computational error. The stability of the structural series was confirmed using five different functionals (see Fig. S3), ensuring that the conclusions are reliable. All these structures are new, with the majority of compositions never before associated with a distinct ODC structure. For the 1 : 5 : 8 and 1 : 3 : 5 compositions, the discovered structures are simply more stable than the literature polytypes, proving the power of the employed high-throughput screening. For instance, at the PBE level of theory, the discovered ground states of CuIn₅Se₈ and CuGa₅Se₈ are about 2.2 and 4.1 meV per atom lower in energy than the most stable literature analog (“Type-D” in Fig. 1c) proposed by Zhang et al.^13,14 The computed decomposition enthalpies (as defined in ref. 62 and implemented in pymatgen,⁴² where negative values correspond to stable compounds) are close to zero (within 0.5 meV per atom) for all except the 1 : 5 : 8 ODC in the series, as shown in Fig. 3c and d. This situation means that the energy stimulus for the formation of other ODCs from a mixture of 1 : 1 : 2 and 1 : 5 : 8 is minimal and, at the same time, there is no substantial hindrance for the reverse decomposition reaction. The compositions at phase equilibrium must thus be sensitive to growth conditions and impurities, as discussed below.

3.2 Structural polymorphism of terminal phases

In the computed convex hull, a series of low-energy 1 : 1 : 2 and 0 : 2 : 3 polytypes were found within 2.5 meV per atom from the ground state. As expected, the chalcopyrite structure is the ground state at [Cu]/[In] = 1 and all other 1 : 1 : 2 polytypes are formed by integration of antiphase domains (often referred to as antisite domain boundaries), yielding mixtures of chalcopyrite-like (CH-type) and CuAu-like (CA-type) layers stacked to form (001)_CH/(100)_CA interfaces, as illustrated in Fig. 4. All of them obey the octet rule and their formation enthalpies depend linearly on the CuAu-type fraction, as we noted earlier for CuInS₂.⁶¹ This structural motif applies to CuGaSe₂ too, although the computed enthalpies there are almost four times higher, implying that CA-type ordering is suppressed with the increasing [Ga]/([Ga] + [In]) (GGI) ratio. These results echo the earlier finding of increasing formation energy of the antisite domain boundary upon Ga alloying⁶³ and justify the absence of the high-temperature sphalerite (disordered zinc-blende) phase of CuGaSe₂,¹⁰ in contrast to the cation-disordered sphalerite δ-CuInSe₂ occurring at above 810 °C.⁹

Fig. 4 Analysis of zinc-blende-derived 1 : 1 : 2 polytypes. (a) Formation enthalpies of 269 polytypes relative to the ground state (ΔH_f) versus the CA-type fraction (defined in the Methods section). The hexagons represent the most stable structures for the given CA-type fraction; all other structures are represented by the crosses. Examples of the hexagon-marked CuInSe₂ polytypes with CA-type fractions of (b) 14%, (c) 25%, (d) 50%, and (e) 75%.

The situation for the In₂Se₃ and Ga₂Se₃ endpoints is less conclusive. Experimentally, In₂Se₃ was observed in several polymorphic forms,^49–53,64 all of which belong to the hexagonal crystal family that could not be captured by our screening. Computationally, the ground state structure for In₂Se₃ is found to depend on the employed functional (see Table S2), but it is always of non-ODC type – the most stable identified ODC-like In₂Se₃ has the above-hull enthalpy of 8.8 and 17.6 meV per atom in the PBE and PBEsol calculations, respectively. In contrast, Ga₂Se₃ crystallizes in two phases, both derived from the zinc-blende lattice via the incorporation of vacancies at every third cationic site. In the high-temperature α-Ga₂Se₃ phase, the vacancy distribution is random.^54,65 In the low-temperature β-Ga₂Se₃ phase, the vacancies are ordered in the form of an ODC-like polytype, either so-called “monoclinic” (space group: Cc)^54,55 or “orthorhombic” (space group: Imm2).⁵⁶ The monoclinic structure is the most stable in PBEsol calculations, but the PBE functional surprisingly predicts a non-ODC ground state with Pm space group symmetry (see Table S2), in agreement with Liu et al.²² While the accuracy of these functionals can be debated, it is certain that a large set (82 in our PBE calculations) of ODC-like polytypes exists within 2.5 meV per atom from the most stable monoclinic β-Ga₂Se₃. This multiplicity of microstates is clearly a factor of disorder, but it might also contribute to stabilizing structures with vacancy ordering over the non-ODC alternatives. On the other hand, this multiplicity is also noted for In₂Se₃, which has 152 structures within 2.5 meV per atom of the most stable ODC-like polytype, yet no studies report synthesis of the zinc-blende-derived In₂Se₃ phase. Further comparative analyses are needed to determine the underlying causes of all highlighted discrepancies.

3.3 Structural trends for ODCs

Having considered the energetic argument for the existence of the topotactic series and polymorphism of the terminal phases, we now turn to the structural relation between ODCs and the chalcopyrite lattice. Fig. 5 shows a chalcopyrite unit cell projected along the [111]_CH axis and seven (near-)stable ODCs. The 0 : 2 : 3 structure is excluded due to the abovementioned uncertainty about the ground states. Slightly less stable structures of 1 : 5 : 8 and 4 : 6 : 11, which have only 0.02 and 0.01 meV per atom above-hull enthalpies in the Cu–In–Se system (PBE results), and some other notable unstable structures, are shown in Fig. S4–S8. The first feature of note is the visually apparent progression of vacancy arrangements with Cu deficiency, coupled with an increasing fraction of Se atoms having κ ≠ 8. These structural changes occur in a way that all depicted ODCs have strictly minimal deviations from the octet rule (see Fig. 2), confirming it as a valid prerequisite of stability. In fact, the most stable instance of a greater deviation from the octet rule (see Fig. S9 and S10, first row) is the 1 : 3 : 5 structure illustrated in Fig. S6 (with 3.3% of κ = 6, 33.3% of κ = 7, 23.3% of κ = 8, and 40% of κ = 9) at 1.3 meV per atom above the convex hull in both systems. Another common feature of low-energy ODCs is the CH-type coordination (see Fig. S9 and S10, second row). Indeed, akin to the 1 : 1 : 2 polytypes, the computed enthalpies of ODCs correlate with the CA-type fraction, although the dependence type (linear or otherwise) was not determined in this case.

Fig. 5 Crystal structures of eight selected (near-)stable structures from the topotactic series. The black lines connect vacancies separated by a (i.e. lattice parameter of the zinc-blende conventional cell) and the shaded grey regions are included to help visualize the vacancy arrangements.

Now, we look closely at specific ODC stoichiometries, starting with 1 : 5 : 8. The most stable 1 : 5 : 8 structure contains cationic vacancies arranged in (001)_CH planes, in accordance with the majority of literature models (see Fig. 1c and Table S1). The projections along the [010]_CH axis for the ground state, a slightly less stable, and the most stable literature 1 : 5 : 8 polytype (“Type-D”) are presented in Fig. S5. Within each vacancy plane, the vacancies are arranged into a square grid with the period of a (i.e. lattice parameter of the zinc-blende conventional unit cell, which is identical to that of the chalcopyrite cell). The planes are separated by a in the normal direction. This structural motif is prevalent for the 1 : 5 : 8 stoichiometry – all 88 CuIn₅Se₈ polytypes and all 17 CuGa₅Se₈ polytypes within 2.5 meV per atom from the convex hulls have vacancies exclusively in such planar arrangements (see Fig. S9 and S10, third and fourth rows). Considering the large number of low-energy polytypes, at least partial cationic disorder should be expected in 1 : 5 : 8 at elevated temperatures, akin to the 1 : 1 : 2 and 0 : 2 : 3 cases discussed above. Moreover, it means that different synthesis routes can yield slightly different ODC structures, explaining the variety of structural models proposed in the literature, as pointed out by Chang et al.¹² It can also be speculated that the extent of disorder is greater for CuIn₅Se₈ than for CuGa₅Se₈ at any temperature as the former has considerably more low-energy microstates. At the same time, the partial disorder should not compromise the overarching planar vacancy arrangement.

The most stable 1 : 3 : 5 structure, in contrast, is conceptually different from the literature one – the cationic vacancies arrange into … → [100]_CH → [010]_CH → [100]_CH → … zigzag arrays (see Fig. S6), not straight channels along the [100]_CH axis as in the structure from ref. 22, 23 and 29. The PBE-computed enthalpies of the literature 1 : 3 : 5 structure are 6.5 and 5.4 meV per atom above the identified ground states in the Cu–In–Se and Cu–Ga–Se systems, respectively. This is not surprising considering that the literature polytype contains κ = 10 tetrahedra (i.e. 60% of κ = 7, 20% of κ = 9, and 20% of κ = 10), whereas our ground state has the minimal deviation from the octet rule (i.e. 40% of κ = 7, 20% of κ = 8, and 40% of κ = 9). Interestingly, the ground state 1 : 3 : 5 structure is devoid of the (001)_CH vacancy planes, despite experimental XRD patterns showing such reflections (see the next section). A dedicated effort to find 1 : 3 : 5 structures with the vacancy planes was made, leading to the discovery of a structural series of 1 : 1 : 2/1 : 5 : 8 multilayers stacked along the [11̄0]_CH direction such that the overall composition matches that of 1 : 3 : 5, as illustrated in Fig. S7. Enthalpies of such structures increase with the multilayer thickness, but crucially, the 1 : 5 : 8 layer can be thick enough to accommodate the vacancy planes. The PBE-computed enthalpy of the most stable vacancy-bearing CuIn₃Se₅ (CuGa₃Se₅) structure is 0.6 (1.1) meV per atom, which is sufficiently low to make these structures accessible at the temperatures of solid-state synthesis and solar cell fabrication.

All intermediate (near-)stable ODCs with 0.5 ≤ [Cu]/[III] < 1.0 have the structural motif best exemplified by Cu₄In₆Se₁₁ in Fig. 5. It consists of large chalcopyrite domains separated by a string of [100]_CH vacancy arrays, which constitute a two-dimensional Cu-poor layer parallel to the (011̄)_CH crystallographic planes. In the Cu-poor regions, all Se ions neighboring vacancies have either κ = 7 or κ = 9, whereas all Se ions inside the chalcopyrite region have the usual κ = 8. Increasing the thickness of the chalcopyrite region produces stable ODCs with higher [Cu]/[III] (e.g. Cu₆In₈Se₁₅ and Cu₈In₁₀Se₁₉), whereas decreasing it yields stable Cu₂In₄Se₇. As such, by stacking chalcopyrite layers of different thicknesses, the 0.5 ≤ [Cu]/[III] < 1.0 compositional range can be spanned with the ODCs.

A categorically distinct motif is observed for the 2 : 4 : 7 and 3 : 5 : 9 ODCs studied by Yarema et al.⁵⁷ These structures were generated under the condition of equidistant vacancies, yielding κ = 10, κ = 6, and/or κ = 5 tetrahedra. These tetrahedra signify excessive deviation from the octet rule, which is reflected in all formation enthalpies relative to the ground states exceeding 15 meV per atom, as listed in Table S1. The instability is slightly less severe for the alternative 3 : 5 : 9 structure, uniquely derived from the wurtzite lattice in ref. 58. Still, the computed formation enthalpies for Cu₃In₅Se₉ and Cu₃Ga₅Se₉ in this structure are over 10 meV per atom above the ground state. Hence, a decision was made to disregard all three structures from ref. 57 and 58 in favor of the more stable polytypes.

Special attention is deserved by the near-stable 10 : 12 : 23 because it corresponds to [Cu]/[III] ≈ 0.83. This composition is at the edge of the single-phase region of α-CIGSe (see Fig. 1) and typical for devices, meaning that the actual absorbers are reminiscent of this structure. The [100]_CH projection of the near-stable 10 : 12 : 23 structure is depicted in Fig. 6a. Just like for the 4 : 6 : 11 ODC described above, this structure is composed of Cu-poor regions sandwiched between chalcopyrite domains. These Cu-poor regions therefore accommodate the majority of the off-stoichiometry in Cu-poor CIGSe, not isolated point defects or neutral [2V_Cu + III_Cu] defect complexes, as assumed by the traditional models.^14,66–68 To estimate the thermodynamic favorability of these arrangements, we compared the PBE-computed enthalpy of the isolated defect complex with that for the discovered 10 : 12 : 23 structure (see the Methods section), finding that the former is 0.3 eV per complex higher in both the Cu–In–Se and Cu–Ga–Se systems. This is not surprising considering that [2V_Cu + III_Cu] complexes increase the deviation from the octet rule, with two κ = 10 tetrahedra being added per complex. In other words, the ODC-like arrangement is energetically favorable at the absolute zero temperature. To probe the thermal evolutions of these arrangements, we also computed the energies of several 10 : 12 : 23 polytypes that were created purposefully in supercells with doubled length along the out-of-plane direction relative to the Cu-poor layers and different separations between them. The PBEsol results in Fig. 6b reveal that the interaction between the Cu-poor layers is weak, with the per-vacancy enthalpy difference being much lower than thermal energy (k_BT ≈ 26 meV at room temperature). As such, the loss of periodicity in the layer stacking is practically inevitable. The actual Cu-poor α-CIGSe must thus occur as a mixture of irregular, partly intertwined Cu-poor regions separating chalcopyrite domains that still contain equilibrium concentrations of [2V_Cu + III_Cu] complexes and other point defects as predicted for the stoichiometric chalcopyrite material.

Fig. 6 (a) The identified ground-state structure of the 10 : 12 : 23 ODC in CIGSe (where the Cu-poor planar defects are evenly spaced by about 30 Å). (b) Relative enthalpy for the 10 : 12 : 23 structures with different separations between the Cu-poor planar defects. These calculations were performed using the PBEsol functional and the fine stability calculation setup (see the Methods section).

Finally, two stable ODCs were found in the 0.2 < [Cu]/[III] < 0.5 region – 4 : 10 : 17 and 3 : 7 : 12 depicted in Fig. 5. These structures are also composed of intertwined chalcopyrite and Cu-poor regions, but it might be more convenient to view them as multilayers of 2 : 4 : 7 and 1 : 3 : 5 ODCs with different thicknesses, as illustrated in Fig. S8. By stacking the layers in different sequences, the structural series could be extended to span the entire 0.2 ≤ [Cu]/[III] ≤ 0.5 range. This exercise, however, needs additional structure generation runs for larger supercells to supplement the series in future studies.

3.4 Simulated diffraction patterns

The reason why different ODC structures are rarely recognized may be that their X-ray/neutron diffraction patterns are similar by virtue of the shared parent lattice. To test this hypothesis, we simulated and compared diffraction patterns for all stable PBEsol-optimized ODCs. The key findings are presented in Fig. 7 and complemented by Fig. S11–S19.

Fig. 7 (a and b) Simulated XRD patterns and (c and d) effective lattice constants of (near-)stable (a and c) Cu–In–Se and (b and d) Cu–Ga–Se ODC structures. The lattice geometries were optimized with the PBEsol functional. The literature data points represent experimental values from ref. 9, 25 and 69–75. The estimated effective lattice parameters for Cu–(In_0.5,Ga_0.5)–Se are presented in Fig. S11. The lattice parameters optimized with different functionals are compared in Fig. S12. The respective simulated neutron diffraction patterns are given in Fig. S13.

The most intense XRD peaks for ODCs are indeed similar to those of the chalcopyrite 1 : 1 : 2, with only a slight shift being visible at higher magnification (see insets of Fig. 7a and b). These peaks originate from the parent zinc-blende lattice and shift due to the known lattice contraction upon copper depletion.⁶⁹ The lattice contraction can be expressed via the effective lattice constant, where Ω is the volume per 8 lattice sites or, equivalently, for the conventional chalcopyrite unit cell. The PBEsol-computed trends and experimental data from ref. 9, 25 and 69–75 are given in Fig. 7c, d and S11. Clearly, the lattice constants decrease linearly with respect to the III₂Se₃ molar fraction, in good agreement with experimental observations. The results delivered by other functionals are less fulfilling (see Fig. S12), which is expected considering that PBEsol was designed to reproduce lattice geometries.³⁷ For the Cu–Ga–Se system, the error in the effective lattice constant is about 2% for PBE and PBE+U calculations, which is significant but not critical considering that the dependance on off-stoichiometry is reproduced. For the Cu–In–Se system though, all functionals except PBEsol fail to capture the lattice contraction, in part because its magnitude is roughly half of that in the Cu–Ga–Se system.

Besides the main parent lattice peaks, the XRD patterns of the ODCs contain a series of weak features at various 2θ angles, attributed to vacancy superstructures. These features are more intense in the Cu–In–Se than in the Cu–Ga–Se system and they are also visible in the simulated neutron diffraction patterns, although with lower intensities (see Fig. S13). The superstructure peaks are the strongest for the 1 : 5 : 8 structures, where the patterns are in agreement with measurements on highly crystalline single-phase β-CIGSe.^11,76,77 Some of these peaks appear as intensified chalcopyrite features (e.g. for CuInSe₂: (011) at 17.1°, (013) at 27.7°, (022) at 34.6°, (121) at 35.6°, and (015)/(123) at 41.9°), while others correspond to otherwise extinct reflections (e.g. for CuInSe₂: (002) at 15.2°, (110) at 21.7°, and (114) at 38.0°).⁷⁸ The appearance of peaks extinct for chalcopyrites is usually explained by the lattice symmetry change to either stannite (space group: I4̄2m) or CuAu-type (space group: P4̄2m),⁷⁵ caused by the long-range ordering of vacancies into the (001)_CH planes.^13,14 Fig. S17 presents simulated XRD patterns for the 10 most stable 1 : 5 : 8 polytypes, all of which have vacancies exclusively arranged in (001)_CH planes (see the previous section). These patterns are indistinguishable and, as a result, averaging them for the total of 89 structures does not smear out the peak intensities. This observation confirms that the superstructure peaks in 1 : 5 : 8 originate from the (001)_CH vacancy planes that withstand partial disordering. Noteworthily, the simulated diffraction patterns are similar for some literature structures as well, especially for those closer to the convex hull (see Fig. S17 and Table S1), suggesting that some of them can still serve as decent ODC models for structural studies.

For comparison, multiple superstructure peaks are seen for the 1 : 1 : 2 polytypes in Fig. S14. In this case, the superstructure peaks are the reflections that appear for polytypes but absent in the “disordered” pattern (simulated for the cation-disordered zinc-blende sphalerite phase). The superstructure peak positions depend on the CA-type fraction, and hence they smear out in the polytype-averaged XRD patterns. Only a couple of superstructure peaks persist after averaging (e.g. at 34.8° and 47.3° for CuInSe₂) as they nearly overlap for different polytypes. Whether these signatures are observable is up for debate, but their low intensity would certainly make for challenging detection.

For the 0 : 2 : 3 ODC-like structures, even more superstructure peaks emerge for polytypes (see Fig. S18), but only a few of them overlap and persist after averaging (e.g. at 14.5°, 29.3°, and 39.0° for fictional In₂Se₃ and at 15.4°, 31.0°, and 41.3° for Ga₂Se₃). Compared with the literature, disordered zinc-blende and the lowest-energy ODC yield XRD patterns similar to the majority of measurements. However, there are also experimental studies with patterns better matched by our polytype-averaged simulations.^79,80 A possible explanation is that those publications produced Ga₂Se₃ with interlaced crystal morphology, akin to that of CuInS₂ nanoparticles and thin films grown by sputtering.^61,81

The above examples illustrate that the impact of structural multiplicity on diffraction patterns is not universal. Fig. 7a and b reveals a series of superstructure XRD peaks, albeit less intense, for (near-)stable ODCs of other (not 1 : 5 : 8) compositions. These peaks appear to evolve from weak chalcopyrite reflections that split into several smaller peaks spanning a progressively wider 2θ range at lower [Cu]/[III] ratios. Surprisingly, their positions are independent of the distance between the Cu-poor regions composing the structures, as shown for 4 : 6 : 11 in Fig. S15 and for 10 : 12 : 23 in Fig. S19, which means that the average patterns retain these reflections. They should therefore be visible in partially disordered Cu-poor α-CIGSe. However, the superstructure peaks might still vanish if the cationic disorder is coupled with an increased CA-type fraction, which is the case for the stoichiometric 1 : 1 : 2 chalcopyrites. Detailed simulations of disordered structures using, for instance, Monte Carlo methods are needed to refine the XRD pattern predictions.

Surprisingly, the XRD reflections from the vacancy planes do not show up in the simulated patterns for the ground state 1 : 3 : 5 structure, nor are they visible for the literature 1 : 3 : 5 structure, as shown in Fig. S16. This finding can be explained by the obvious lack of vacancy planes in those structures (see Fig. 5), but it contradicts a well-established experimental fact – these peaks are present for both the 1 : 3 : 5 and 1 : 5 : 8 compositions.^11,76,77,82 Looking at the XRD patterns for other 1 : 3 : 5 structures, we see that a special multilayer-type low-energy polytype depicted in Fig. S7b exhibits such reflections (see Fig. S16) because it contains vacancy planes within the constituent 1 : 5 : 8-like domains. Since the above-hull enthalpy of this structure is only 0.6 and 1.1 meV per atom in the Cu–In–Se and Cu–Ga–Se systems, respectively, it might still form locally within β-CIGSe grains. It is also possible that this multilayer structure is stabilized by some overlooked factor, like vibrational entropy. However, when analyzing diffraction patterns for the rest of this series (see Fig. S7), we note that the vacancy plane reflections appear only for specific multilayer thicknesses, likely due to the alternating constructive/destructive interference from the adjacent 1 : 5 : 8 domains. While it is still possible that only those polytypes that produce constructive interference form in actuality, we acknowledge that this assumption lacks evidence and thus requires further analysis.

4. Discussion

4.1 Structural trends

Our calculations reveal that CIGSe occurs in many closely related structural forms with different arrangements of cationic vacancies. These structures underpin the gradual changes in lattice parameters with respect to the [Cu]/[III] ratio.^25,83–85 The existence of the ODC series, where each component is stable and has nearly zero decomposition enthalpy, explains why both α-CIGSe and β-CIGSe have wide homogeneity regions, which make the material tolerant to off-stoichiometry in a practical sense⁸⁶ but sensitive to processing conditions when it comes to the crystal structure.¹¹ It is even more important, however, to recognize these structures as representatives of a continuum that can be spanned by stacking ODC layers in various constellations. This mechanism is reminiscent of the idea of charge-neutral off-stoichiometric superclusters proposed by Rockett,⁸ except that the Cu-poor regions discovered herein cannot be produced by aggregation of defect complexes. The structural continuum allows the absorber to undergo topotactic transitions (varying the [I]/[III] ratio without changing the lattice type) during the deposition process, which typically includes periods of Cu-poor growth, whether during three-stage co-evaporation or metal precursor selenization. In the completed devices, β-CIGSe is often (if not always) present as either a thin layer at the heterojunction^85,87–90 or embedded grains,^86,91–94 forming low-energy coherent interfaces with α-CIGSe due to small lattice mismatch, especially at low GGI.^95,96 The resulting morphologies inevitably impact the absorber properties while being almost invisible to X-ray and neutron diffraction, as revealed by our simulated patterns. The only sign of ODCs in XRD are the weak superstructure peaks, normally noticeable for well-ordered crystals of β-CIGSe, but not recognizable in typical polycrystalline thin films. However, even in the best circumstances, the superstructure peak assignment is ambiguous, as evident from the variety of structural β-CIGSe models.^{11,21,25,76,97,98} In this regard, the discovered stable structures can serve as more accurate base models for future Rietveld refinements and short-range order analyses.^12,70

4.2 Thermodynamic stability

Despite the overall consistency of the topotactic ODC continuum with experimental observations, one serious discrepancy stands out. Experimentally, both Cu₂Se–In₂Se₃ and Cu₂Se–Ga₂Se₃ phase diagrams feature a broad two-phase (α + β) miscibility gap (see Fig. 1). This gap means that ODCs with [I]/[III] ≈ 0.45–0.75 are unstable, in contradiction to the predicted stability of the entire continuum. Our explanation for this (laid out below) is based on the entropic contribution that is not explicitly calculated but can be qualitatively understood based on structural trends. In the future, a full study of thermal effects using Monte Carlo simulations would help complete the picture.

The main sources of entropy in non-magnetic gapped materials are vibrations and configurational disorder. For the comparison of ODC structures with the common parent lattice, similar to disordered alloys,^99,100 vibrational entropy is expected to have relatively small effect on phase boundaries. As the aim here is to capture the general topology of phase diagrams, it is the configuration part that we focus on. Fig. 8 shows our tentative estimation for the configuration entropy contribution (−T·ΔS, where T is the absolute temperature) to the free energy (ΔG = ΔH_f − T·ΔS, where ΔH_f is the computed formation enthalpy), based on the Boltzmann definition of entropy (ΔS = k_B·ln W, where W is the number of microstates). The sketched entropy curve has two energy minima dictated by the structural trends as follows:

Fig. 8 Schematics of the free energy curve for off-stoichiometric CIGSe. The enthalpy and entropy contributions are illustrated separately, with the latter proposed as the source of the empirical miscibility gap.

• At [Cu]/[III] = 1, neglecting the contribution from stoichiometric defects (such as the antisite domain boundaries), the configurational entropy approaches zero because the entire bulk has the same ordered chalcopyrite structure (only one microstate is present, W = 1).

• Near [Cu]/[III] ≈ 0.8, where the off-stoichiometry is accommodated by Cu-poor planar defects, the configuration entropy is higher because the separation between the Cu-poor regions is random, owing to the weak interaction between them. In addition, the chalcopyrite domains are large enough to accommodate [2V_Cu + III_Cu] complexes at random, which raises the configurational entropy further.

• Near [Cu]/[III] ≈ 0.5, which corresponds to the 2 : 4 : 7 composition, the configurational entropy is low again because the Cu-poor planar defects are tightly packed and no variation in separation is possible. Moreover, the small volume of chalcopyrite domains constrains the formation of [2V_Cu + III_Cu] complexes.

• Near [Cu]/[III] ≈ 0.2, the configurational entropy is higher again due to the large multiplicity of low-energy structures, evident from the large number of low-energy 1 : 5 : 8 structures in both Cu₂Se–In₂Se₃ and Cu₂Se–Ga₂Se₃ systems. This logic applies all the way down to [Cu]/[III] = 0, which also exhibits structural multiplicity (although the shape of the entropy curve there is less important for the present discussion).

The energy minima in the sketched free energy curve correspond to α-CIGSe and β-CIGSe, with the equilibrium compositions defined using the common tangent. The region between them is the empirical miscibility gap. Note that outside of the miscibility gap, all ODC structures comprising α-CIGSe and β-CIGSe belong to the same topotactic ODC series described in the previous sections.

This explanation is not without contradictions, however. It logically follows that the miscibility (α + β) gap should narrow at lower temperatures as the entropy contribution dwindles. In contrast, Herberholz et al. concluded a wider miscibility gap for the Cu–In–Se ingots produced at 100 °C compared with those obtained at 500 °C (i.e. both single-phase regions are narrower at 100 °C).^101,102 Similar assumptions have been made when sketching low-temperature sections of phase diagrams elsewhere.^9,10 While the empirical evidence is limited, we recognize the possible contradiction and emphasize the need for further analysis.

4.3 Sensitivity of phase boundaries

Since the small configuration entropy determines the phase diagrams, as described in the previous section, any extra energy source can have a strong impact on the location of phase boundaries. One possible example is the incorporation of alkali elements, which tend to pile up in the β-phase of pristine and Ag-alloyed CIGSe absorbers.^86,103,104 The preferential accumulation can lower the energy of one CIGSe phase more than that of the other, shifting the delicate energy balance and thereby altering the phase equilibrium. Indeed, 0.1 at% Na addition was shown to expand the single-phase region of α-CISe,^101,102 and changes in the β-CIGSe phase amount were noted upon alkali incorporation.^105,106 Bearing in mind that the energy contributions brought about by different alkalis are not identical, it should also come as no surprise that more β-CIGSe patches form in the chalcopyrite films grown on potassium-rich high-strain-point glass than on the sodium-rich soda-lime glass substrate.⁹¹ Besides, because the thermodynamic stimulus for decomposition in the (α + β) range is just the small entropy gain, the phase composition may be extra sensitive to non-equilibrium processing. It has indeed been shown that the single-phase chalcopyrite region of α-CIGSe is in effect extended by quenching.^107,108 Moreover, colloidal synthesis of Cu–In–Se nanoparticles was reported to produce compounds with compositions in the miscibility gap,^57,109 which could be due to the stabilization effects of free surfaces, disorder, impurities, and/or relatively low growth temperature (about 280–310 °C). These discrete pieces of evidence support the governing role of configuration entropy, although a proper quantitative examination of the outlined theory is certainly needed.

4.4 Defect chemistry outlook

The ODC continuum discovered here invokes fundamental change to the established view on the defect chemistry of chalcopyrite materials. As a start, the term “Cu vacancies” should no longer be used to describe off-stoichiometry because stable ODCs cannot be derived by incorporation of [2V_Cu + III_Cu] complexes. This is evident from the misfit between the chalcopyrite domains in Cu₁₀In₁₂Se₂₃ depicted in Fig. 6a – multiple manipulations are needed to produce this structure (see Fig. S20). Hence, the association of cationic vacancies with any element is unreasonable – they must be treated as an integral part of the crystal structure, as lattice sites in their own right. The near-stoichiometric ODC-like structure of 10 : 12 : 23 is thus the best description of Cu-poor CIGSe. All atoms filling vacancies in this structure would constitute interstitial defects, not substitutional. This claim applies not only to external dopants but also to Cu ions, which are expected to yield electrically active donors when occupying the neutral cation vacancy site, most likely in the “+1” charged state (judging by the nominal valences). The bulk defect chemistry of CIGSe^15,110,111 in this case still applies, but only locally, within the chalcopyrite domains, where uncompensated Cu vacancies can still be the dominant acceptors as in the stoichiometric material. These bulk defects would interplay electronically with defects in the Cu-poor regions, defining the macroscopic Fermi level. This mechanism might hold the keys to understanding metastabilities and degradation in CIGSe absorbers.^112,113 Mobility of point defects through the vacancy arrays may also be enhanced because the ODC-like structure provides channels for fast ion transport. This process in particular could be behind the faster mass transport observed in Cu-poor CIGSe absorbers.¹¹⁴

Alkali metal ions added during the pre-, peri-, or post-deposition treatments, as well as those in-diffusing from the glass substrates, could also incorporate preferentially into the Cu-poor planar defects, filling the neutral vacancy sites. This would not be surprising considering the propensity for Na segregation in other off-stoichiometric imperfections, such as grain boundaries, dislocations, and nanoclusters in bulk CIGSe.^115,116 Similar to alkali interstitials in bulk chalcopyrites,^117,118 the resulting alkali-on-vacancy defects are expected to act as donors in the “+1” charged state. This prediction is in stark contrast with stoichiometric CIGSe, which strongly favors charge-neutral substitutional alkali-on-copper defects.^117–119 The “+1” defect charge is also likely for alkali ions occupying vacancy sites in [2V_Cu + III_Cu] complexes but the resulting [Na_Cu + III_Cu + V_Cu] complexes are projected to have higher formation enthalpy and hence occur in lower concentrations, due to the extra 0.3 eV needed to form the [2V_Cu + III_Cu] complex in the first place. Either way, the charged alkali-on-vacancy species would be able to change the defect chemistry makeup, especially for the alkali-loaded thin β-CIGSe layers at heterojunctions.⁸⁹ This hypothesis has been crudely assessed for unstable (“Type-D”) CuIn₅Se₈ before,¹⁹ suggesting that alkali elements prefer to form stable donors. On a similar note, other extrinsic dopants could favor formation in proximity of the Cu-poor planar defects. For example, focusing on the same unstable (“Type-D”) CuIn₅Se₈ polytype, the solubility limit for Cd inside β-CIGSe due to in-diffusion from the CdS buffer was estimated to be in the order of 1% at low temperature (under the idealized assumption of a metallic Cd source).^16,17 With the knowledge of the stable ODCs, an extended study of point defects in Cu-deficient CIGSe can be conducted in a systematic manner. The presence of electrically active alkali-related dopants in ODCs, if confirmed, would help reconcile the observations of varying doping densities¹²⁰ and enhanced optoelectronic characteristics¹²¹ with the charge-neutral nature of the dominant alkali-on-copper defects in stoichiometric CIGSe.^117–119

In addition, while alkali ions in stoichiometric CIGSe have low migration barriers,^4,6,118 the barriers along the vacancy arrays in the Cu-poor regions are likely even lower. Ion movement in this scenario would not only contribute to fast Na diffusion,^117,122 but also explain reversible in-diffusion under prolonged voltage bias (i.e. potential-induced degradation)^123,124 and the suspected electric polarization yielding hysteresis behavior due to ionic conductivity in Ag-alloyed CIGSe.^86,112

5. Summary

We performed a high-throughput computational investigation of off-stoichiometric structures in Cu(In,Ga)Se₂ (CIGSe) and proposed an entirely new picture that resolves several lasting discrepancies and offers a new direction for materials engineering. By screening over 100 000 zinc-blende-derived Cu–In–Se and Cu–Ga–Se structures, we discovered a topotactic series of stable and near-stable ordered defect compounds (ODCs) forming a continuum with respect to [Cu]/[III] in both systems. These ODCs are responsible for the wide single-phase regions of α-CIGSe and β-CIGSe in the phase diagram, making CIGSe tolerant to Cu deficiency. The stable ODC structures are stabilized by different constellations of vacancies in the parent lattice, which are enabled by re-arrangement of cations in a way that minimizes deviation from the octet rule. The structural kinship of ODCs makes them challenging to detect by diffraction methods, with the exception of Cu(In,Ga)₅Se₈ which yields several unique peaks. Together with the small decomposition enthalpies, the challenging detection is the likely reason why the ODC series has remained obscured until now. Still, the experimental lattice contraction in Cu-poor absorbers is closely reproduced by the topotactic ODC continuum, which overall provides a better description of the structural evolution in Cu-poor CIGSe than the conventional model based on [2V_Cu + III_Cu] defect complexes. Nevertheless, the defect complexes can still play an important role in explaining the emergence of the miscibility gap between α-CIGSe and β-CIGSe and in enabling a proper description of the defect chemistry of the former. The discovery of the ODC continuum is thus imperative for more accurate calculations of the electronic properties (see part II of this study), phase transitions, and defect chemistries of chalcopyrite absorbers in the practically relevant, Cu-deficient composition ranges. It also inspires a similar search for ODCs in other functional materials, with kesterites being the most obvious candidates. The generated knowledge is expected to stimulate further improvements in thin-film solar cells via better control over the deposition and defect engineering in absorber films.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

The most stable structures discovered here have been deposited at the Materials Cloud Archive repository and can be accessed via the following identifier: https://doi.org/10.24435/materialscloud:j9-15.⁴⁵ Supplementary information (SI): complementary figures and tables. See DOI: https://doi.org/10.1039/d5ta07043k. Other raw data, processing scripts, and more detailed descriptions can be provided by the corresponding author upon reasonable request.

Acknowledgements

This work was partially supported by the Wallenberg Initiative Materials Science for Sustainability (WISE) funded by the Knut and Alice Wallenberg Foundation. The authors are also thankful for the financial support from the Swedish Foundation for Strategic Research (grant no. RMA15-0030) and the strategic research area STandUP for Energy. The computations and data handling were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS), partially funded by the Swedish Research Council through grant agreement no. 2022-06725. The authors also thank Dr Jes K. Larsen and Dr Oleksandr I. Malyi for the insightful scientific discussions.

References

1. J. Keller, K. Kiselman, O. Donzel-Gargand, N. M. Martin, M. Babucci, O. Lundberg, E. Wallin, L. Stolt and M. Edoff, Nat. Energy, 2024, 9, 467–478.
2. S. Siebentritt, L. Gütay, D. Regesch, Y. Aida and V. Deprédurand, Sol. Energy Mater. Sol. Cells, 2013, 119, 18–25.
3. H. Zhao, M. Kumar and C. Persson, Phys. Status Solidi C, 2012, 9, 1600–1603.
4. L. E. Oikkonen, M. G. Ganchenkova, A. P. Seitsonen and R. M. Nieminen, J. Appl. Phys., 2013, 114, 083503.
5. L. E. Oikkonen, M. G. Ganchenkova, A. P. Seitsonen and R. M. Nieminen, J. Appl. Phys., 2013, 113, 133510.
6. T. Maeda, A. Kawabata and T. Wada, Jpn. J. Appl. Phys., 2015, 54, 08KC20.
7. S. Nakamura, T. Maeda and T. Wada, Jpn. J. Appl. Phys., 2013, 52, 04CR01.
8. A. Rockett, Thin Solid Films, 2000, 361–362, 330–337.
9. U.-C. Boehnke and G. Kühn, J. Mater. Sci., 1987, 22, 1635–1641.
10. J. C. Mikkelsen, J. Electron. Mater., 1981, 10, 541–558.
11. S. Lehmann, D. F. Marrón, M. León, R. Feyerherm, E. Dudzik, E. J. Friedrich, M. Tovar, Y. Tomm, C. Wolf, S. Schorr, Th. Schedel-Niedrig, M. Ch. Lux-Steiner and J. M. Merino, J. Appl. Phys., 2011, 109, 013518.
12. C.-H. Chang, S.-H. Wei, J. W. Johnson, S. B. Zhang, N. Leyarovska, G. Bunker and T. J. Anderson, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 68, 054108.
13. S. B. Zhang, S.-H. Wei and A. Zunger, Phys. Rev. Lett., 1997, 78, 4059–4062.
14. S. B. Zhang, S.-H. Wei, A. Zunger and H. Katayama-Yoshida, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 9642–9656.
15. J. Pohl and K. Albe, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 245203.
16. J. Kiss, T. Gruhn, G. Roma and C. Felser, J. Phys. Chem. C, 2013, 117, 10892–10900.
17. J. Kiss, T. Gruhn, G. Roma and C. Felser, J. Phys. Chem. C, 2013, 117, 25933–25938.
18. H. Xiao and W. A. Goddard, J. Chem. Phys., 2014, 141, 094701.
19. E. Ghorbani, J. Kiss, H. Mirhosseini, G. Roma, M. Schmidt, J. Windeln, T. D. Kühne and C. Felser, J. Phys. Chem. C, 2015, 119, 25197–25203.
20. M. Malitckaya, T. Kunze, H.-P. Komsa, V. Havu, E. Handick, R. G. Wilks, M. Bär and M. J. Puska, ACS Appl. Mater. Interfaces, 2019, 11, 3024–3033.
21. C. D. R. Ludwig, T. Gruhn, C. Felser and J. Windeln, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 174112.
22. W. Liu, H. Liang, Y. Duan and Z. Wu, Phys. Rev. Mater., 2019, 3, 125405.
23. X. Chen, W. Liu and Y. Duan, J. Phys.: Condens. Matter, 2020, 33, 075401.
24. A. Sharan, F. P. Sabino, A. Janotti, N. Gaillard, T. Ogitsu and J. B. Varley, J. Appl. Phys., 2020, 127, 065303.
25. T. Maeda, W. Gong and T. Wada, Jpn. J. Appl. Phys., 2016, 55, 04ES15.
26. F. Jiang and J. Feng, Appl. Phys. Lett., 2006, 89, 221920.
27. S. Kumar, S. Joshi and S. Auluck, Mater. Chem. Phys., 2015, 162, 372–379.
28. L.-H. Tu, N. T. T. Tran, S.-K. Lin and C.-H. Lai, Adv. Energy Mater., 2023, 13, 2301227.
29. A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder and K. A. Persson, APL Mater., 2013, 1, 011002.
30. S. Wei and A. Zunger, J. Appl. Phys., 1995, 78, 3846–3856.
31. S.-H. Han, C. Persson, F. S. Hasoon, H. A. Al-Thani, A. M. Hermann and D. H. Levi, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 74, 085212.
32. D. W. Houck, E. I. Assaf, H. Shin, R. M. Greene, D. R. Pernik and B. A. Korgel, J. Phys. Chem. C, 2019, 123, 9544–9551.
33. K. V. Sopiha, J. Keller, C. Persson, J. J. S. Scragg, C. Platzer-Björkman and M. Edoff, J. Mater. Chem. A, 2025 10.1039/D5TA07043K.
34. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
35. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
36. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
37. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
38. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 1505–1509.
39. F. Tran and P. Blaha, Phys. Rev. Lett., 2009, 102, 226401.
40. J. Sun, A. Ruzsinszky and J. P. Perdew, Phys. Rev. Lett., 2015, 115, 036402.
41. A. V. Krukau, O. A. Vydrov, A. F. Izmaylov and G. E. Scuseria, J. Chem. Phys., 2006, 125, 224106.
42. S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson and G. Ceder, Comput. Mater. Sci., 2013, 68, 314–319.
43. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.
44. A. van de Walle, Calphad, 2009, 33, 266–278.
45. K. Sopiha, J. Keller, C. Persson, J. Scragg, C. Platzer-Björkman and M. Edoff, Materials Cloud Archive, 2025, 2025, 134.
46. S.-H. Wei, L. G. Ferreira and A. Zunger, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 45, 2533–2536.
47. S.-H. Wei, S. B. Zhang and A. Zunger, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, R2478–R2481.
48. B. Krebs, A. Schiemann and M. Läge, Z. Anorg. Allg. Chem., 1993, 619, 983–988.
49. K. Osamura, Y. Murakami and Y. Tomiie, J. Phys. Soc. Jpn., 1966, 21, 1848.
50. A. Likforman and M. Guittard, C. R. Acad. Sci., 1974, 279, 33–35.
51. A. Likforman, D. Carré and R. Hillel, Acta Crystallogr., Sect. B: Struct. Sci., 1978, 34, 1–5.
52. A. Pfitzner and H. D. Lutz, J. Solid State Chem., 1996, 124, 305–308.
53. M. Küpers, P. M. Konze, A. Meledin, J. Mayer, U. Englert, M. Wuttig and R. Dronskowski, Inorg. Chem., 2018, 57, 11775–11781.
54. D. Lübbers and V. Leute, J. Solid State Chem., 1982, 43, 339–345.
55. G. Ghémard, S. Jaulmes, J. Etienne and J. Flahaut, Acta Crystallogr., Sect. C: Cryst. Struct. Commun., 1983, 39, 968–971.
56. G.-Y. Huang, N. M. Abdul-Jabbar and B. D. Wirth, J. Phys.: Condens. Matter, 2013, 25, 225503.
57. M. Yarema, N. Yazdani, O. Yarema, N. Đorđević, W. M. M. Lin, D. Bozyigit, S. Volk, A. Moser, A. Turrini, P. A. Khomyakov, M. Nachtegaal, M. Luisier and V. Wood, Adv. Mater., 2024, 36, 2406351.
58. A. Moser, O. Yarema, N. Rusch, N. Đorđević, W. M. M. Lin, D. Bozyigit, N. Yazdani, M. Yarema, M. Luisier and V. Wood, ACS Nanosci. Au, 2024, 5, 21–28.
59. S. Kirklin, J. E. Saal, B. Meredig, A. Thompson, J. W. Doak, M. Aykol, S. Rühl and C. Wolverton, npj Comput. Mater., 2015, 1, 1–15.
60. M. Hellenbrandt, Crystallogr. Rev., 2004, 10, 17–22.
61. J. K. Larsen, K. V. Sopiha, C. Persson, C. Platzer-Björkman and M. Edoff, Adv. Sci., 2022, 9, 2200848.
62. C. J. Bartel, A. W. Weimer, S. Lany, C. B. Musgrave and A. M. Holder, npj Comput. Mater., 2019, 5, 4.
63. K. Park, B.-H. Jeong, H. Y. Lim and J.-S. Park, J. Appl. Phys., 2021, 129, 025703.
64. T. Gödecke, T. Haalboom and F. Sommer, J. Phase Equilib., 1998, 19, 572.
65. H. Hahn, W. Klingler and Z. Anorg, Chem, 1949, 259, 135–142.
66. S. Chen, A. Walsh, X.-G. Gong and S.-H. Wei, Adv. Mater., 2013, 25, 1522–1539.
67. C. Stephan, S. Schorr, M. Tovar and H.-W. Schock, Appl. Phys. Lett., 2011, 98, 091906.
68. C. Stephan, T. Scherb, C. A. Kaufmann, S. Schorr and H.-W. Schock, Appl. Phys. Lett., 2012, 101, 101907.
69. C. Rincón, S. M. Wasim, G. Marín, J. M. Delgado and J. Contreras, Appl. Phys. Lett., 2003, 83, 1328–1330.
70. E. Haubold, P. Schöppe, S. Eckner, S. Lehmann, I. Colantoni, F. d'Acapito, F. di Benedetto, S. Schorr and C. S. Schnohr, J. Alloys Compd., 2019, 774, 803–812.
71. J. L. M. de Vidales, H. Charbonnier, R. Diaz, J. M. Merino, F. Rueda and M. Leon, Jpn. J. Appl. Phys., 2000, 39, 336.
72. J. M. Merino, J. L. M. de Vidales, S. Mahanty, R. Diaz, F. Rueda and M. Leon, J. Appl. Phys., 1996, 80, 5610–5616.
73. J. M. Merino, S. Mahanty, M. Leon, R. Diaz, F. Rueda and J. L. M. de Vidales, Thin Solid Films, 2000, 361–362, 70–73.
74. C. Stephan, S. Schorr and H. W. Schock, MRS Proc., 2009, 1165, M09.
75. C. Stephan, Structural trends in off stoichiometric chalcopyrite type compound semiconductors, PhD thesis, Helmholtz-Zentrum Berlin, 2011.
76. W. Paszkowicz, R. Lewandowska and R. Bacewicz, J. Alloys Compd., 2004, 362, 241–247.
77. S. Lehmann, D. F. Marrón, M. Tovar, Y. Tomm, C. Wolf, S. Schorr, T. Schedel-Niedrig, E. Arushanov and M. C. Lux-Steiner, Phys. Status Solidi A, 2009, 206, 1009–1012.
78. M. Souilah, A. Lafond, N. Barreau, C. Guillot-Deudon and J. Kessler, Appl. Phys. Lett., 2008, 92, 241923.
79. J. K. Larsen, O. Donzel-Gargand, K. V. Sopiha, J. Keller, K. Lindgren, C. Platzer-Björkman and M. Edoff, ACS Appl. Energy Mater., 2021, 4, 1805–1814.
80. K. George, C. H. (Kees) de Groot, C. Gurnani, A. L. Hector, R. Huang, M. Jura, W. Levason and G. Reid, Chem. Mater., 2013, 25, 1829–1836.
81. X. Shen, E. A. Hernández-Pagan, W. Zhou, Y. S. Puzyrev, J.-C. Idrobo, J. E. Macdonald, S. J. Pennycook and S. T. Pantelides, Nat. Commun., 2014, 5, 5431.
82. S. Yamazoe, H. Kou and T. Wada, J. Mater. Res., 2011, 26, 1504–1516.
83. M. Souilah, A. Lafond, C. Guillot-Deudon, S. Harel and M. Evain, J. Solid State Chem., 2010, 183, 2274–2280.
84. M. León, J. M. Merino and E. Arushanov, Mold. J. Phys. Sci., 2006, 5, 373–381.
85. A. Stokes, M. Al-Jassim, A. Norman, D. Diercks and B. Gorman, Prog. Photovolt.: Res. Appl., 2017, 25, 764–772.
86. J. Keller, L. Stolt, K. V. Sopiha, J. K. Larsen, L. Riekehr and M. Edoff, Sol. RRL, 2020, 4, 2000508.
87. O. Cojocaru-Mirédin, Y. Fu, A. Kostka, R. Sáez-Araoz, A. Beyer, N. Knaub, K. Volz, C.-H. Fischer and D. Raabe, Prog. Photovolt.: Res. Appl., 2015, 23, 705–716.
88. Y. Zhao, S. Yuan, Q. Chang, Z. Zhou, D. Kou, W. Zhou, Y. Qi and S. Wu, Adv. Funct. Mater., 2021, 31, 2007928.
89. S. Ishizuka and P. J. Fons, Phys. Rev. Appl., 2021, 15, 054005.
90. T. Nishimura, S. Toki, H. Sugiura, K. Nakada and A. Yamada, Prog. Photovolt.: Res. Appl., 2018, 26, 291–302.
91. K. V. Sopiha, J. K. Larsen, O. Donzel-Gargand, F. Khavari, J. Keller, M. Edoff, C. Platzer-Björkman, C. Persson and J. J. Scragg, J. Mater. Chem. A, 2020, 8, 8740–8751.
92. J. R. Tuttle, D. S. Albin and R. Noufi, Sol. Cells, 1991, 30, 21–38.
93. J. Keller, N. Martin, O. Donzel-Gargand, K. Kiselman, U. Zimmermann, L. Stolt, C. Platzer-Björkman and M. Edoff, Sol. RRL, 2024, 8, 2301018.
94. J. Keller, K. V. Sopiha, O. Stolt, L. Stolt, C. Persson, J. J. S. Scragg, T. Törndahl and M. Edoff, Prog. Photovolt.: Res. Appl., 2020, 28, 237–250.
95. D. Schmid, M. Ruckh, F. Grunwald and H. W. Schock, J. Appl. Phys., 1993, 73, 2902–2909.
96. M. A. Contreras, H. Wiesner, J. Tuttle, K. Ramanathan and R. Noufi, Sol. Energy Mater. Sol. Cells, 1997, 49, 239–247.
97. W. Hönle, G. Kühn and U.-C. Boehnke, Cryst. Res. Technol., 1988, 23, 1347–1354.
98. T. Hanada, A. Yamana, Y. Nakamura, O. Nittono and T. Wada, Jpn. J. Appl. Phys., 1997, 36, L1494.
99. A. M. Holder, S. Siol, P. F. Ndione, H. Peng, A. M. Deml, B. E. Matthews, L. T. Schelhas, M. F. Toney, R. G. Gordon, W. Tumas, J. D. Perkins, D. S. Ginley, B. P. Gorman, J. Tate, A. Zakutayev and S. Lany, Sci. Adv., 2017, 3, e1700270.
100. C. K. Gan, Y. P. Feng and D. J. Srolovitz, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 235214.
101. R. Herberholz, H. W. Schock, U. Rau, J. H. Werner, T. Haalboom, T. Godecke, F. Ernst, C. Beilharz, K. W. Benz and D. Cahen, in Conference Record of the Twenty Sixth IEEE Photovoltaic Specialists Conference - 1997, 1997, pp. 323–326.
102. R. Herberholz, U. Rau, H. W. Schock, T. Haalboom, T. Gödecke, F. Ernst, C. Beilharz, K. W. Benz and D. Cahen, EPJ Appl. Phys., 1999, 6, 131–139.
103. J. Keller, P. Pearson, N. Shariati Nilsson, O. Stolt, L. Stolt and M. Edoff, Sol. RRL, 2021, 5, 2100403.
104. J. Keller, H. Aboulfadl, L. Stolt, O. Donzel-Gargand and M. Edoff, Sol. RRL, 2022, 6, 2200044.
105. C. Stephan, D. Greiner, S. Schorr and C. A. Kaufmann, J. Phys. Chem. Solids, 2016, 98, 309–315.
106. T. Kodalle, L. Choubrac, L. Arzel, R. Schlatmann, N. Barreau and C. A. Kaufmann, Sol. Energy Mater. Sol. Cells, 2019, 200, 109997.
107. T. Gödecke, T. Haalboom and F. Ernst, Z. Metallkd., 2000, 91, 622–634.
108. D. Colombara, H. Elanzeery, N. Nicoara, D. Sharma, M. Claro, T. Schwarz, A. Koprek, M. H. Wolter, M. Melchiorre, M. Sood, N. Valle, O. Bondarchuk, F. Babbe, C. Spindler, O. Cojocaru-Miredin, D. Raabe, P. J. Dale, S. Sadewasser and S. Siebentritt, Nat. Commun., 2020, 11, 3634.
109. O. Yarema, M. Yarema and V. Wood, Chem. Mater., 2018, 30, 1446–1461.
110. M. Malitckaya, H.-P. Komsa, V. Havu and M. J. Puska, Adv. Electron. Mater., 2017, 3, 1600353.
111. S. Lany and A. Zunger, J. Appl. Phys., 2006, 100, 113725.
112. P. Pearson, J. Keller, L. Stolt and C. Platzer-Björkman, Phys. Status Solidi B, 2023, 260, 2300170.
113. P. Pearson, J. Keller, L. Stolt, M. Edoff and C. Platzer-Björkman, Phys. Status Solidi B, 2022, 259, 2200104.
114. D. J. Schroeder, G. D. Berry and A. A. Rockett, Appl. Phys. Lett., 1996, 69, 4068–4070.
115. P.-P. Choi, O. Cojocaru-Mirédin, R. Wuerz and D. Raabe, J. Appl. Phys., 2011, 110, 124513.
116. A. Karami, M. Morawski, H. Kempa, R. Scheer and O. Cojocaru-Mirédin, Sol. RRL, 2024, 8, 2300544.
117. Z.-K. Yuan, S. Chen, Y. Xie, J.-S. Park, H. Xiang, X.-G. Gong and S.-H. Wei, Adv. Energy Mater., 2016, 6, 1601191.
118. M. Malitckaya, H.-P. Komsa, V. Havu and M. J. Puska, J. Phys. Chem. C, 2017, 121, 15516–15528.
119. H. Aboulfadl, K. V. Sopiha, J. Keller, J. K. Larsen, J. J. S. Scragg, C. Persson, M. Thuvander and M. Edoff, ACS Appl. Mater. Interfaces, 2021, 13, 7188–7199.
120. C. P. Muzzillo, Sol. Energy Mater. Sol. Cells, 2017, 172, 18–24.
121. O. Ramírez, M. Bertrand, A. Debot, D. Siopa, N. Valle, J. Schmauch, M. Melchiorre and S. Siebentritt, Sol. RRL, 2021, 5, 2000727.
122. A. Laemmle, R. Wuerz, T. Schwarz, O. Cojocaru-Mirédin, P.-P. Choi and M. Powalla, J. Appl. Phys., 2014, 115, 154501.
123. V. Fjällström, P. Szaniawski, B. Vermang, P. M. P. Salomé, F. Rostvall, U. Zimmermann and M. Edoff, IEEE J. Photovolt., 2015, 5, 664–669.
124. V. Fjällström, P. M. P. Salomé, A. Hultqvist, M. Edoff, T. Jarmar, B. G. Aitken, K. Zhang, K. Fuller and C. K. Williams, IEEE J. Photovolt., 2013, 3, 1090–1094.

---