Atomistic models of In and Ga diffusion in Cu(In,Ga)Se₂

The need for renewable energy sources has increased interest in solar energy technologies. Chalcopyrite Cu(In,Ga)Se₂ (CIGS) is a promising candidate absorber material for thin-film photovoltaic devices, as CIGS solar cells have shown the highest light-to-electricity conversion ever reported of any thin-film photovoltaic technology.¹ CIGS displays a high tolerance for crystalline disorder, which allows for significant control of the electronic properties of the material.² Indeed, CIGS [being a solid solution of CuInSe₂ (CIS) and CuGaSe₂ (CGS)] has a bandgap that can be tuned from 1.0 eV (pure CIS) to 1.7 eV (pure CGS) by controlling the In/Ga content. Commercially available CIGS absorbers typically are designed with alloy gradients, with the $[ Ga ] / ( [ Ga ] + [ In ] )$ (GGI) ratio varying with depth in the absorber to optimize carrier collection. Understanding the diffusion mechanisms of In and Ga in both CIGS alloys is, therefore, critical to optimizing cell performance.^1,3

There are numerous mechanisms that could potentially mediate the diffusion of group III species In and Ga in CIGS. In a crystal, diffusion proceeds through three primary mechanisms: (1) interchange by rotation, (2) migration through interstitial sites, or (3) diffusion via vacant lattice sites.⁴ A number of experimental studies have found evidence that the group III species In and Ga diffuse primarily via vacancies on the metal (group I: Cu, and group III) sublattices of CIGS.^4–6 More recent computational studies have confirmed this belief, finding that interstitial In_int is very unfavorable and unlikely to form, reducing the possibility of an interstitial diffusion mechanism.^7,8 Going further, computational studies have also found that V_In has a high formation energy and will not exist in great abundance, while V_Cu defects are expected to be much more prevalent.^7,8 This suggests the possibility of group III diffusion via vacant V_Cu sites, but only if group III species form group I antisites in high enough numbers. Reference 7 predicts this indeed happens, with all extra In atoms eventually expected to form In_Cu defects. Additional theoretical studies have confirmed a similar abundance of Ga_Cu defects in CIS, along with an abundance of In/Ga_Cu in CGS.^9,10 Building on these findings, we propose a model of group III diffusion, whereby In_Cu and Ga_Cu antisites diffuse via V_Cu defects.

In this work, we use density functional theory (DFT) to determine binding energies of Cu-sublattice defect complexes and show that they are well-modeled as coulombic interactions of point charges. We then attempt to reconcile our results with discrepancies in calculations of these binding energies found in previous works. We determine the migration barriers for group III diffusion on the Cu-sublattice, and justify the marked differences observed between In and Ga. We present both kinetic lattice Monte Carlo (KLMC) models and analytical models for group III diffusion in CIGS and use them to elucidate the nature of the different mechanisms involved. Last, we compare the results of these models to experiment and find strong consistency between them.

We calculated binding energies and migration barriers using Vienna Ab Initio Simulation Package (VASP) based on the projector augmented-wave method (PAW) and using a plane-wave basis set.¹¹ The energy cutoff of the plane-wave set was 385 eV, and the generalized gradient approximation (GGA) type Perdew–Burke–Ernzerhof (PBE)¹² exchange-correlation functional was employed. For binding energy calculations, we also employed hybrid Heyd–Scuseria–Ernzerhof (HSE06)¹³ exchange-correlation functionals, with the portion of the Hartree–Fock exchange set to $α=0.25$⁠, and an inverse screening length of $ω= 0.2 Å − 1$ (the default values for HSE06, having been shown to be suitable for most systems¹³). Most calculations were performed with a 64-atom cell and $4×4×4$ Monkhorst–Pack $k$-point mesh, while some calculations required a 128-atom cell with $4×4×2$ (if $c$-direction was doubled) or $2×4×4$ (if $a$-direction was doubled) $k$-point mesh. The extrapolations using 16-, 128-, and 432-atom cells also used a $4×4×2$ mesh. Ionic positions were relaxed until intra-ionic Hellmann–Feynman forces fell below 0.01 eV/Å, while the cell volume and shape were fixed to their ideal CIS or CGS values. The calculated lattice constants for CIS and CGS are shown in Table I and compared with experimental values. Our structural parameters show reasonably good agreement with experiment, with the $c$ axis generally larger in our calculations. For the calculation of migration energy barriers in VASP, the climbing-image nudged elastic band (CINEB) method was used.¹⁴ 

We took a different approach than in these previous works, by demonstrating that the binding of In/Ga_Cu and V_Cu defects in CIGS can be modeled as coulombic interactions of point charges. As a result, the energy of an arbitrary configuration of these defects can be predicted from the electrostatic potential energy of an equivalently charged and spatially distributed array of point charges in a medium with static dielectric constant equal to that of CIGS. Such a model is computationally much more efficient than DFT for computing the energy of a given atomic configuration. The energy can be calculated almost instantaneously on a computing system that would ordinarily take several hours for a similar calculation using DFT. The efficiency of such a model is particularly useful for large KLMC simulations, as employed later in this work.

In order to demonstrate that the interactions can be modeled as coulombic, we calculated the energies of 128-atom cells using DFT with defect pairs at varying degrees of separation. If coulombic interactions can account for the binding (or repulsion) between defects, the change in DFT energy upon separating a pair of defects should correspond to the change in the Madelung energy of an empty cell containing only the isolated charged defects undergoing the same separation. To calculate the Madelung energy, we modeled In/Ga_Cu as $+2$ and V_Cu as $−$1 point charges in an empty cell, with the static dielectric constant of bulk CIS/CGS. The Ewald summation solver provided by Pymatgen was used for the calculations.¹⁸ Note that for charged systems, both VASP and Pymatgen impose a neutralizing background charge. The choice of the static dielectric constant proved to be a challenging issue. In particular, PBE-GGA functionals are known to overestimate the static dielectric constant relative to experimental values.^7,19 Thus, we cannot use experimental values directly, as they would not reflect the internal behavior of the DFT system to which we are comparing. Instead, the dielectric constant must be calculated from DFT, so that it can adequately capture the screening of coulombic interactions that occur within the DFT framework being studied. However, attempts to calculate the static dielectric constant for CIS/CGS with PBE-GGA functionals failed to produce stable converged results, and so we simply used best-fit values of 17 for CIS and 15 for CGS. These can be compared with reported experimental values for CIS, which range from 11.3 to 13.6,^20,21 and for CGS, which range from 11 to 11.5.^22–25 In contrast, Ref. 8 calculated a value for the dielectric constant of $≈10.65$ in CIS using HSE06 hybrid functionals, which is much closer to the experimental range. Therefore, we also performed DFT calculations using HSE06 functionals and scaled the corresponding Madelung energies by the dielectric constant 10.65. Due to computational constraints, we did not relax the ionic degrees of freedom in the HSE06 runs, but only relaxed the electronic degrees of freedom starting with fully relaxed PBE-GGA runs.

The results, shown in Fig. 1, show a close match between DFT and Madelung energies, demonstrating that the binding of In/Ga_Cu and V_Cu defects can be modeled by coulombic interactions alone. Even though for the PBE-GGA calculations we adjusted the dielectric constant to fit the Madelung energies to the DFT energies, this does not necessarily invalidate the results, as they still demonstrate that the energy variation with separation distance is well-modeled by considering only coulombic interactions. The only major deviation is in (b) for PBE-GGA calculations at first-nearest-neighbor (1NN) separation, but as this deviation does not appear in the equivalent and more accurate HSE06 calculation, it may be a unique feature of PBE-GGA functionals and not indicative of the real physical behavior. Note that in this work, we use “1NN” and related terms such as “2NN” to refer to nearest neighbors on the Cu-sublattice specifically, where “1NN” Cu sites are actually 2NN sites separated by Se atoms. We find that In and Ga have nearly identical behavior, despite differences in ionic radii. The CGS binding energies were also calculated but are not shown as they are very similar to the results in Fig. 1. The CGS binding energies tend to be slightly stronger than their CIS counterparts, explained by CGS having a slightly lower dielectric constant and slightly smaller lattice constants, both of which contribute to increased coulombic interaction. The differences between the PBE-GGA and HSE06 calculations in Fig. 1 are due to the difference in dielectric constant. An important consequence of our findings is that, due to the slow decay of the $r − 1$ coulombic energy dependence, we should expect long-range defect interactions, leading to loosely bound defect clusters, rather than tightly bound and well-defined defect complexes. This will impact not only the equilibrium configuration of defects but also the In/Ga diffusion mechanisms, as discussed later in Sec. III.

The results for PBE-GGA calculations are shown in Fig. 2, along with the binding energies found by using only the Madelung energies of equivalent supercells containing only the charged defects. For this, we took the Madelung energy of a supercell containing the complex composed of $n$ defects and subtracted the Madelung energies of $n$ supercells each containing one individual defect, ignoring perfect supercells as they contain no charged defects. Note that a value of $κ=17$ was used to scale the Madelung energies, which we found was appropriate for PBE-GGA, as previously explained. The results give converged values of $−$0.43 eV for In_Cu–V_Cu and $−$0.74 eV for In_Cu–2V_Cu. From the results it is clear that as supercell size increases, the values found with DFT converge to the values found with just the Madelung energy [which at infinite cell size is simply given by Eq. (2)]. The large deviations between the results found with DFT vs those found with the Madelung energy at small cell size can be explained by considering the different scaling of the two approaches. The Madelung energy scales like $L − 1$ as observed. However, Eq. (3) includes a term that scales like $L − 3$⁠, which accounts for the interaction of the localized charge distribution with the uniform background²⁸ and explains the large increase at small cell sizes. At infinite cell size, the two approaches converge to the same values and, thus, provide additional evidence for a model of binding of Cu-sublattice defects in CIGS that is coulombic.

Table II shows values for In/Ga_Cu–V_Cu and In/Ga_Cu–2V_Cu complexes using Eq. (2) (note that the use of the lower experimental value for $κ$ gives stronger binding energies in Table II than reported in Fig. 2). Our values for CIS are much stronger than those found in Ref. 10, but our values for CGS are very similar to theirs. Some difference may arise from their sampling of only a single cell size. Although they applied the correction scheme of Lany and Zunger,²⁹ the corrected energy may still vary significantly with cell size.²⁸ Still, it is not clear why their values match ours for CGS but not for CIS, or why they find such a large variation between CIS and CGS, where we find little.

In Ref. 7, they also report weaker antisite/vacancy binding in CIS than we find. However, the value for the In_Cu–2V_Cu binding energy in CIS is reported as the additional change in energy due to bringing a second V_Cu defect near a bound In_Cu–V_Cu complex. For consistency with our results, the value should be reported as $−$0.46 eV: the sum of the initial binding energy of the In_Cu–V_Cu pair and the energy change from adding a second V_Cu defect. This brings their value more in line with ours. Reference 7 used the finite-size-scaling extrapolation scheme of Ref. 27 to obtain converged binding energies. However, as Ref. 27 explains, the extrapolation can only be done with cells of the same symmetry (such as cubic cells of increasing size). In contrast to our application of this method, Ref. 7 used cells of widely varying symmetry (in terms of 16-atom unit cells: $2×1×1$⁠, $2×2×1$⁠, $3×3×1$⁠, $3×3×2$⁠, and $4×4×2$⁠), and consequently their extrapolated values are not necessarily reliable. Despite this, their results for separating In_Cu–V_Cu pairs using PBE-GGA functionals produce energy differences of $≈0.3$ eV, in line with values we find for a very similar procedure (see Fig. 1).

Reference 17 reports binding energies of the In/Ga_Cu–2V_Cu complexes in CIS much closer to ours, although still weaker. Because the focus of their work is comparing the binding energies of In vs Ga complexes, rather than identifying the absolute magnitude of the energies, they did not attempt to extrapolate their results even though they are still changing markedly with increasing cell size (despite using the correction scheme of Freysoldt et al.³⁰). This suggests that after extrapolating to infinite cell size, they may find stronger binding energies than currently reported, bringing them even more in line with our values. Additionally, they concur with our results that In and Ga have similar binding energies.

As discussed in Sec. I, previous experimental and computational work leads us to focus on a diffusion mechanism for In and Ga in CIGS dominated by 1NN vacancy-mediated hops of In_Cu and Ga_Cu antisites on the Cu-sublattice. This mechanism involves an In_Cu or Ga_Cu hopping onto a neighboring vacant V_Cu site. To calculate the migration barriers of these processes in CIS and CGS, DFT CINEB calculations were performed for vacancy-antisite exchanges. Figure 3 shows the energy along the diffusion path taken by the moving antisite, and Table III shows the resulting migration barriers for the different processes. Our In_Cu and Ga_Cu barriers agree very well with those of Refs. 3 and 17, and our V_Cu barriers agree with previous literature as well.^2,7,31 Although In_Cu and V_Cu have similar diffusion barriers, it is still expected that Cu diffusion occurs much faster than In diffusion. This arises due to most In in the crystal being immobile on group III sites, so the overall diffusivity of In is significantly lower than the In_Cu antisite diffusivity, a point discussed further in Sec. III. Two important observations can readily be made. First, the Ga process in CIS has a much higher barrier than the In process in CIS and similarly the Ga process in CGS has a much higher barrier than the In process in CGS. Second, the In process in CIS has a higher barrier than the same In process in CGS, which is likewise true for the Ga process in CIS compared to the Ga process in CGS, and for the Cu process in CIS compared to the Cu process in CGS. Both observations can be explained by quantifying the “average distortion” that the lattice undergoes during the migration. To determine the average distortion of a given state, we calculated the distance every atom in a simulation cell is displaced from its ideal relaxed position. From this, we determined the average displacement distance of an atom in the cell by averaging the displacements at all sites (excluding the sites involved in the exchange). Last, we found the change in average displacement between the ground and saddle states of each process, which we will refer to as the “average distortion” caused by the process. The results are presented in Table III. Positive values of average distortion indicate a greater degree of average displacement in the saddle state compared to the ground state. We expect that a larger value of average distortion will correspond to a larger barrier, as the system must undergo a greater change from ideality for the defect to migrate, which requires greater energy. In fact, this is exactly the trend that is observed. In both CIS and CGS, the In_Cu and Cu_Cu exchanges have similar barriers to each other and correspondingly have similar values of average distortion. Meanwhile, in both CIS and CGS, the Ga_Cu exchange has a significantly higher energy barrier than the equivalent In_Cu and Cu_Cu processes, and likewise has a significantly higher value of average distortion. Furthermore, all CIS processes have higher barriers than their equivalent CGS processes and likewise the CIS processes have significantly greater average distortions than their CGS counterparts.

Figure 4 shows the saddle and ground state configurations of Ga_Cu–V_Cu and In_Cu–V_Cu swaps in CIS. All six processes reported in Table III have qualitatively the same saddle point configuration. The migrating atom moves through the tetrahedral interstitial site, but at the saddle point, it is distorted away from the true tetrahedral site, toward the two vacancies. The distance between the migrating atom and Se1 remains relatively constant throughout the migration so that the migrating atom is effectively “pivoting” around this selenium. We find that for all six processes, the saddle point is symmetric between the initial and final states. This contrasts with Ref. 3, who find an asymmetric saddle point configuration for the Ga_Cu swap in CIS. This discrepancy may arise from differences in their calculation method from ours that are not reported, such as the lattice constants used and the degrees-of-freedom allowed to relax. They hypothesize that the large difference between Ga and In barriers in CIS may have some connection to this asymmetry, but since we find similar energy barrier differences with wholly symmetric configurations, the asymmetry may not be a key factor. As previously explained, we quantified the amount each atom was displaced from its ideal relaxed location. For both In and Ga migration processes, Se1 undergoes by far the greatest increase in displacement between the ground and saddle states, in line with its role as the “pivot” for the migrating atom. However, the average displacement of atoms in the In process ground state is much greater than that of the Ga process ground state, while the average displacement of atoms in the two saddles is more comparable. Because the In process begins already further from ideality, moving to the saddle does not require as large of an average distortion as for the Ga process, which begins closer to ideality in a more stable configuration (hence, the higher migration barriers found for Ga compared to In processes). This trend is also observed in CGS, despite Ga being the native lattice species. Of the nine nearest neighbors labeled in Fig. 4, only Se2, Se3, and Se4 behave significantly differently between the In and Ga processes. As these three are all neighbors to the ground state (and in the case of Se2, not a neighbor to the saddle state), this confirms that the difference between In and Ga processes originates primarily with the ground state.

In addition to the 1NN hops, previous studies have found that second-nearest-neighbor (2NN) hops of V_Cu have a comparable migration barrier to 1NN V_Cu hops and, therefore, should be considered. We use “2NN” to refer specifically to the second-shortest distance between Cu sites, not to lattice connectedness. Of special interest is the situation where the initial and final sites of a 2NN V_Cu hop are both 1NN to the same antisite. Although there is no overall change in binding energy, the expected path approaches very closely to the antisite, and so the possibility of a distinctly different barrier from the isolated 2NN hop exists. Table IV shows the barriers for 2NN hops of V_Cu in the presence of different antisites and native copper. Our value of 1.18 eV for a free 2NN V_Cu hop in CIS agrees well with Ref. 2. In CIS, the presence of an antisite at 1NN increases the barrier only modestly, but it has a more substantial effect in CGS. This is likely due to the smaller lattice constants of CGS compared to CIS, which brings the antisite closer to the migrating atom and raises the energy. In both CIS and CGS, when an antisite is present, the migrating atom is distorted away from the antisite at the saddle point. In both cases, the saddle point configuration is asymmetric when an antisite is present, and the potential energy landscape around the saddle point is very flat and broad. Because 2NN hops only exist in $ab$ planes (see Fig. 5), the inclusion of 2NN hops in our models may impact anisotropy, as discussed in Sec. III 2NN hops of In_Cu and Ga_Cu were found to have such high barriers that they do not significantly impact group III diffusion.

To account for alloyed CIGS, we considered how the formation energies and migration barriers of the relevant defects vary with GGI. Our DFT calculations show no energy change when exchanging $G a Cu + I n III → I n Cu + G a III$ or $( G a Cu + V Cu ) 1 NN+ I n III → ( I n Cu + V Cu ) 1 NN + G a III$⁠, suggesting In_Cu and Ga_Cu antisites form with equal favorability. References 10 and 31 have confirmed that In_Cu and Ga_Cu have formation energies similar to each other in both CIS and CGS, and that the formation energy of V_Cu is the same in both CIS and CGS. Likewise, Refs. 35 and 36 have found that in bulk group III rich (thus Cu poor) CIS or CGS, the formation energy of V_Cu is unchanged over a wide range of Fermi energies in both CIS and CGS, and likewise the formation energy of In_Cu in CIS is similar to that of Ga_Cu in CGS. Taken together, these findings suggest that defect concentrations are unaffected by changes in GGI outside of the change in antisite ratio reflected by Eq. (8). Additionally, previous Monte Carlo studies suggest that for temperatures above the range of interest considered here (⁠ $> 400 K$⁠), In and Ga mix homogeneously on the group III sublattice. Therefore, we are justified in only considering system-level GGI and ignoring the possibility of local GGI variation.^37–39 However, because there is a large variation in migration barriers between CIS and CGS, for the given GGI, we used a linear interpolation of the barriers between their CIS and CGS end points. This assumption was confirmed by running several NEB calculations with randomly arranged In and Ga atoms (corresponding to a GGI of 0.375) on group III sites. The average barrier for the In_Cu–V_Cu swap was nearly identical to the value expected from a linear interpolation, likewise for the Ga_Cu–V_Cu swap.

A difficulty that arose was developing a fast and efficient algorithm to calculate the binding energy of the system, given our coulombic binding model. The long-range $r − 1$ decay of coulombic interactions requires a sufficiently long cutoff to achieve well-converged results, but this can be computationally intensive. We implemented the approach outlined by Ref. 26, which requires only finding the electrostatic potentials at the two sites on which the process occurs, including only the stationary charges up to a certain cutoff. The cutoff is defined radially from the point halfway between both sites, so that the total charge used to evaluate both sites is unchanged. This has been shown to speed convergence.²⁶  Figure 6 shows convergence of the calculated diffusivity of In for increasing cutoff. As can be seen, above 20 angstroms the results show tight convergence. All the following calculations used 20 angstrom cutoffs to ensure the best results.

We calculated diffusivities for In and Ga in CIS over a range of compositions and temperatures using our KLMC model. We chose to report Ga diffusivity in CIS as opposed to in CGS in order to have a direct comparison between In and Ga diffusivity, and due to the lack of experimental data for diffusivity in very high Ga content CIGS. To be completely physically realistic, a very small number of Ga atoms should have been added to the simulations to reflect a dilute Ga content that does not alter the overall properties of the bulk CIS system. However, it was found to be difficult to obtain sufficient statistics to estimate the MSD with such low Ga content. We decided instead to artificially convert all antisites to Ga while still maintaining CIS parameters, in order to obtain more robust MSD data. Since our simulation has both In_Cu and Ga_Cu bind via the same electrostatic interactions, the equilibrium configuration of defects should only depend on the overall concentration of vacancies and antisites, not the ratio of Ga_Cu/In_Cu. In turn, as the hop rate of an individual diffusing Ga_Cu antisite only depends on temperature and the neighboring defect configuration, conversion of In_Cu to Ga_Cu should not affect Ga diffusivity. A run with half Ga_Cu and half In_Cu antisites confirmed that no change in Ga diffusivity occurs from a change in the ratio of Ga_Cu/In_Cu.

Figure 7 shows results from the simulations for different temperatures and compositions. Figure 7(a) depicts the diffusivity of antisites, which is simply the diffusivity of the indicated species extracted from the KLMC simulations, as given in Eq. (5). Figure 7(b) gives the overall diffusivity, which is the antisite diffusivity multiplied by the fraction of the diffusing species that is mobile. This takes into account the fact that most of the In and Ga in the system exists frozen on the group III sublattice and does not contribute to diffusion; thus, the overall diffusivity is lower than the antisite diffusivity. It can be readily observed that the diffusivity of antisites [Fig. 7(a)] is mostly independent of the cation ratio, suggesting that the diffusing behavior of an individual antisite is not strongly affected by the density of other antisites and vacancies nearby, at least within these composition ranges. Therefore, the compositional dependence of the overall diffusivity of In and Ga in CIGS [Fig. 7(b)] is primarily due to the availability of mobile antisites, which increases in more Cu-poor material. Reference 5 finds a severalfold decrease in the diffusivity of Ga in CIS at 1000 K when increasing the k_Cu/k_In ratio from 0.75 to 0.95, which is generally consistent with our findings. It is also evident from our results that Ga diffusion is predicted to be a couple orders of magnitude lower than In diffusion in CIGS having the same composition, which is consistent with experimental observations. We discuss this in greater detail below.

The activation energies for In and Ga diffusion were found from the average of the 3 Arrhenius slopes for different compositions. In diffusion was found to have an activation energy of 1.24 eV and Ga diffusion a value of 1.54 eV. These activation energies are quite reasonable. For In diffusion, the activation energy of 1.24 eV is noticeably higher than the migration barrier of 1.1 eV as listed in Table III, whereas for Ga the value of 1.54 eV is closer to the barrier of 1.47 eV. When an antisite swaps with a vacancy, if the exchange does not bring the antisite into contact with a new vacancy, it is likely to simply reverse the swap with the original vacancy, leading to no net diffusion. For long-range diffusion to occur, after an antisite-vacancy exchange, the vacancies must re-orient themselves to provide a pathway for a subsequent exchange that does not merely reverse the prior swap.⁴⁰ Because there are a variety of mechanisms for this to occur, it is difficult to determine a single barrier to this re-orientation process. These mechanisms include an antisite capturing a vacancy from far away, a vacancy disassociating from an antisite, and a vacancy orbiting around an antisite to a different 1NN location. Furthermore, the rates of all these mechanisms are influenced by the local defect landscape. A benefit of KLMC is its ability to capture all these complex behaviors with a relatively simple model. Despite this complexity, even using just NEB calculations of a few representative processes suggests that these mechanisms generally fall in the range of 1.2–1.3 eV. For In diffusion, the low migration barrier causes the In antisites to swap back and forth many times before a vacancy can re-orient itself on a new site. Therefore, the overall diffusion of In is limited by the effective “vacancy re-orientation” barrier, readily explaining why the activation energy is in the range of 1.2–1.3 eV. With Ga however, the migration barrier is so high that after a Ga_Cu–V_Cu exchange, the vacancies will re-orient many times before the Ga exchanges again, and subsequent Ga hops will be equally likely in any direction. Thus, the overall diffusion of Ga is limited by Ga_Cu–V_Cu exchange, and so the activation energy is nearly the exchange energy.

It was also observed that when considering $a$⁠, $b$⁠, and $c$ axis diffusion separately, the diffusivity in the $c$ direction was generally between 1.5 and 2 times higher than that in the $a$ or $b$ direction. As mentioned previously, because the Cu-sublattice of CIGS is identical to the diamond lattice elongated by a factor of about $2$ in the $c$ direction, each 1NN hop moves the defect proportionately more along the $c$ direction than along the $a$ or $b$ directions. Thus, it is expected that about a twofold higher diffusivity would be observed in the $c$ direction. For In, the doubling of the $c$ direction diffusivity is depressed to only a 1.5–1.9 factor increase. This is believed to be due to the effect of 2NN V_Cu hops, which provide additional pathways for vacancy reorientation in the $ab$ plane only, boosting $ab$ diffusivity by increasing the frequency of uncorrelated hops. For Ga, a simple twofold $c$ direction increase is observed, explained by Ga diffusion being primarily limited by Ga/V exchange rather than V reorientation. These results suggest the possibility of anisotropic diffusion in CIGS, which is particularly interesting as we have not identified any work that recognizes such behavior in CIGS.

We defined two free parameters for the effective activation energies of 1V and 2V+ complexes and fit the analytical model to the KLMC results for different temperatures and compositions by minimizing the total percent error. For In, the solution converged on identical activation energies for 1V and 2V+ of 1.26 eV, which is also nearly identical to the overall activation energy of 1.24 eV found directly from the KLMC simulations. The average percent error in the diffusivity from the analytic vs the KLMC model is only 2.9%. As explained previously, because the diffusion of In is limited by the effective “vacancy-reorientation” barrier, not the antisite-vacancy exchange, both 1V and 2V+ complexes have the same limiting barrier, so the results we found are expected.

For Ga, we found a 1V effective activation energy of 1.47 eV and a 2V+ energy of 1.75 eV. The average percent error in the diffusivity is 4.4%, similar to the In diffusivity error. The difference of effective activation energy here is also reasonable. The 1V activation energy is identical to the 1V antisite-vacancy exchange barrier because overall diffusion of Ga is limited by the antisite-vacancy exchange barrier. The 2V+ antisite-vacancy exchanges have higher barriers because they require changes in binding energy that increase the barrier by about 0.15 eV for a 2V complex hop and even higher for 3V+ hops, in line with the value of 1.75 eV found for 2V+ complex activation energies. An interesting consequence of this higher activation energy for 2V+ complexes is that diffusion of Ga is dominated by 1V complex hopping. Even at 300 K, when our model predicts 98.7% of antisites should be in 2V+ complexes, diffusion is surprisingly still dominated by 1V complexes, which account for 99.8% of the diffusivity at this temperature. Figure 7 plots the analytic models against the KLMC results, which show very good agreement.

To validate our models against experimental results, we performed simulations to match the conditions of the In radiotracer diffusion study by Ref. 43. In this work, they report temperatures for each run and a GGI value of 0.3. They do not give values for valence ratio or cation ratio, but the process cited for the fabrication of their samples⁴⁴ is known to produce CIGS with valence ratio of nearly 1 and cation ratio of 0.85.^45–48 Therefore, we found it a reasonable comparison to use these values of valence ratio, cation ratio, and GGI in our simulations. Figure 9 shows the comparison of our KLMC models with their experimental results. We find an excellent agreement between our model and the reported data. Because our KLMC model relies solely on values calculated from ab initio methods (with the exception of the experimental dielectric constant, which itself falls within a narrow range of values and is not observed to strongly influence diffusivity in the simulations), these results provide strong evidence for the validity of our model for In diffusion in CIGS. Although we could not find studies on tracer diffusion of Ga in CIGS, there exist several studies on various types of in- and inter-diffusion of Ga. While these studies show a wide range of reported values, most suggest a diffusivity of Ga much smaller than that of In, consistent with our findings shown in Fig. 7.^{5,6,45,49–51} Reference 45 finds an activation energy for Ga in-diffusion from GaAs into single-crystal CIS of 1.49 eV, very close to our value of 1.54 eV for Ga tracer diffusion. Future work will involve developing continuum models to reproduce the specific design of these in- and inter-diffusion experiments, allowing for a more direct comparison between model and experiment.

In conclusion, we have developed models for In and Ga diffusion in CIGS starting from first-principles. We have calculated the binding energies of defect complexes existing on the Cu-sublattice of CIGS and have determined that the binding can be modeled as electrostatic interactions of point charges. We have found the migration barriers for processes occurring on the Cu-sublattice as well, explaining the much higher Ga diffusion barrier by quantifying the change in lattice distortion during the migration. We have developed KLMC simulations for In and Ga diffusion in CIS, determining the activation energies for diffusion and finding minimal dependence on composition. We find that In migration is characterized by a high probability of reversals, while Ga diffusion is dominated by uncorrelated hops.

Based on analysis of the mechanisms revealed, we developed an analytic model for In and Ga diffusion suitable for use in a continuum simulation and have shown that it predicts complex formation and diffusivity as seen in the KLMC simulations with great accuracy. Last, we have compared our models with the experimental results of Ref. 43 for In diffusion and have found strong agreement. Likewise, our model has predicted diminished diffusivity of Ga relative to In by some orders of magnitude, a result echoed by experimental studies.

This work was supported in part by the U.S. Department of Energy’s Office of Energy Efficiency and Renewable Energy (EERE) under the Solar Energy Technology Office Award Number DE-EE0008556. This work was facilitated through the use of the Hyak supercomputer system at the University of Washington, funded by the UW Student Technology Fee. This material was based in part upon work supported by the state of Washington through the University of Washington Clean Energy Institute.

The authors have no conflicts to disclose.

Aaron Gehrke: Conceptualization (supporting); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Software (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). David Sommer: Conceptualization (supporting); Methodology (supporting); Supervision (supporting); Writing – review & editing (equal). Scott Dunham: Conceptualization (lead); Funding acquisition (lead); Methodology (equal); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.