Sequential quenching to predict semiconductor defect concentrations from formation & migration energies: the case of CdTe:As doping

Abstract

Defect concentrations in semiconductors are strongly influenced by thermal history during growth and cooldown, yet most defect calculations assume either instantaneous quenching from high temperature or that full-equilibrium is maintained – two limiting cases rarely approached in reality. Here, we introduce sequential quenching (SQ) as a 3rd type of defect calculation utilizing defect formation and migration energies to model defect concentrations subject to diffusion-limited kinetics in samples cooled at finite rates. In SQ, the concentration of each defect is frozen at a characteristic temperature determined by its diffusion rate, distance to sources/sinks, and cooling rate. Because different charge-states interact through charge neutrality but freeze at different temperatures, the sequence of freeze-in events is non-commuting. Critically, not all room-temperature SQ solutions can be predicted from full equilibrium (EQ) or full-quenching (FQ) calculations – erroneous predictions are likely without SQ. We illustrate SQ using the example of As-doped CdTe, for which experimental data show differences in doping with cooling rate and between polycrystalline thin-films for photovoltaics and bulk crystals. SQ calculations reveal that fast-diffusing defects such as Cd-interstitials remain mobile to lower temperatures and freeze-in at larger characteristic distances, leading to strong compensation and n-type behavior in rapidly cooled or bulk samples. Slower cooling and reduced characteristic distances suppress donor freeze-in and enhance p-type activation. These results establish SQ as a physically transparent and computationally efficient framework for connecting cooling conditions, sample geometry, and defect kinetics to dopant activation in CdTe and related materials.

---

Introduction

The properties of materials are strongly influenced by the concentrations of point defects and defect complexes, which govern carrier density, compensation, and transport behavior in semiconductors.^1–3 The sensitivity of semiconductors to impurities and native defects at ppm–ppb concentrations is unparalleled in material properties and explains why semiconductor processing is so demanding. Populations of defects are not static; they evolve in response to the thermal and chemical history during crystal growth, irradiation, post-growth thermal processing, and even device operation over long times. Their evolution is governed not only by short-range diffusion and complex formation, but also by long-range transport from/to sources/sinks like dislocations, grain boundaries, surfaces, interfaces, and other extended defects. Therefore, bulk crystals and thin films may exhibit different defect ensembles, as may homoepitaxial thin films vs. heteroepitaxial ones having higher dislocation counts. The prime example discussed herein is the widespread observation that efficiently doping CdTe and CdSeTe polycrystalline thin films p-type using group-V elements (grV = N, P, As, Sb) is more challenging than for melt-grown, larger crystals. Another example is that radiation damage in GaN pn diodes grown on freestanding GaN templates can self-anneal at room temperature, while damage was permanent for diodes grown on sapphire with copious threading dislocations. The interpretation is that the dislocations served as sinks for interstitials, preventing Frenkel pair annihilation.^4–10

During crystal growth, local thermal and chemical equilibrium – as for bulk crystallization from the melt – or at least steady-state dynamic equilibrium – as for vapor phase growth processes like molecular beam epitaxy (MBE) or organometallic vapor phase epitaxy (OMVPE) – can be assumed at the growth interface. Usually, the best assumption one can make is that the growth interface and its defect numbers reach full- or dynamic- equilibrium at the highest temperature experienced, since it offers the highest chance of overcoming all kinetic barriers. It is not guaranteed, however, that full chemical and kinetic equilibrium will be reached; for example, it is common to grow structures involving Mg acceptor doping of AlGaInN nitride semiconductors in OMVPE and then to subject them to an activation anneal at a higher temperature in inert gas to remove the passivating H atoms.^11–13 These materials may be thermodynamically unstable at these temperatures in the absence of extreme pressures or reactive nitrogen sources, but the process is optimized to kinetically limit decomposition.

Since diffusion and other kinetic processes become slower as temperature decreases, the next question is “How do details of the cooling process affect defect populations?” Point defects may diffuse to sample surfaces, annihilate each other (as in vacancy-interstitial Frenkel pairs),¹⁰ be generated or annihilated by edge dislocation climb,¹⁴ form complexes like Coulomb-bound donor-acceptor pairs,^15,16 condense into clusters (e.g., voids or oxygen-related precipitates in Si¹⁷), or cross phase boundaries. The extent of all of these processes depends on the time spent at each temperature – the T(t) trajectory. The rates of most kinetic processes are Arrhenian, thus the total extent of processes like diffusion or reaction varies inversely with cooling rate. In principle, each sample having different dimensions and microstructural sources/sinks would require a unique, multidimensional drift-diffusion-Poisson-reaction simulation to predict its final distributions of defects. Such in-depth studies are exactly what is needed to understand specific devices and processes,^18–20 but this level of specificity and expensive computation is anathema to making high-level, more generalizable predictions. Additionally, it is rare for all of the needed kinetic parameters and boundary conditions to be available, adding uncertainty to the results. Therefore, a defect calculation approach that approximates defect kinetics and the effects of a coarse-grained description of sample size and microstructure is desired in order to better predict outcomes for a variety of real-world samples.

Herein, we introduce sequential quenching (SQ) as a method for approximating defect concentrations taking into account defect thermodynamics, diffusion kinetics, thermochemical history, and the characteristic distance separating surfaces or other sinks/sources (e.g., crystal size, film thickness, grain size, or distance between dislocations). Conventional calculations of defect concentrations typically assume either full thermodynamic equilibrium (EQ) during infinitely slow cooling or complete freeze-in during instantaneous full quenching (FQ). EQ conceptually corresponds to cooling rate γ = 0 or all kinetic barriers being <k_BT for all temperatures, while FQ corresponds to γ = ∞ or the limiting kinetic barrier being >>k_BT for all temperatures. Interstitials with low migration barriers maintain EQ even to below room temperature (e.g., Li in Si), while vacancy-mediated diffusers like substitutional dopants frequently exhibit behavior close to FQ. These limiting cases fail to represent real-world growth and processing conditions, for which cooling rates may span several orders of magnitude and defect evolution is governed jointly by thermodynamics, kinetics, and sample geometry.^21,22 Kröger referred to these situations as various types of partial equilibrium – one may assume that a subset of defects has fixed concentrations below threshold temperatures during processing, while others continue to equilibrate based on the thermodynamic conditions. SQ was developed to allow estimation of defect concentrations for finite γ between the two typical limits and as a function of characteristic distance.

Fig. 1(a) illustrates how the concentrations (denoted by [ ] in the text) of generic neutral defects 1 and 2 would behave during cooling under FQ, SQ, and EQ assumptions. Despite all starting from the same initial conditions at the highest temperature, each scenario results in a unique ensemble of defects at room temperature. SQ is illustrated for both slow and rapid cooling; larger concentrations of each defect (larger supersaturation) are maintained for rapid cooling. Fig. 1(b) summarizes the key thermodynamic and kinetic assumptions underlying the three calculation approaches. Under EQ, the cooling behavior of neutral defects is governed only by their own thermodynamics. In EQ, FQ, and SQ cases, we assume that charge carriers equilibrate at all temperatures.

Fig. 1 (a) Concentration-temperature trajectories for generic, uncharged defects 1 and 2 having different formation energies under assumptions of equilibrium (EQ), full quenching (FQ), and sequential quenching (SQ) at slow and fast cooling rates. In SQ, each defect (or charge state) remains equilibrated down to some freezing-in temperature below which the kinetics are too slow to allow further changes. Arrhenius kinetics result in higher freezing temperatures for faster cooling rates (γ). Impurities at concentrations well below solubility limits can be modeled as fixed-concentration defects, which is equivalent to assuming FQ or SQ with T_freeze at T ≥ T_max. (b) Table outlining thermodynamic and kinetic aspects of each calculation procedure: SQ uses a unique combination of EQ and FQ for each defect or charge state. (c) Concentrations of donor-like D and acceptor-like A₁, A₂ defects, as well as band electrons and holes (n, p) for an EQ scenario taking charge balance into account, and in which all defects have equilibrium concentrations less than the intrinsic carrier density n_i at all temperatures. At room temperature, the sample would be intrinsic, meaning that n and p dominate charge balance. (d) As for (c), under an SQ scenario in which A₁ freezes in before A₂ and D remains able to equilibrate at room temperature (e.g., D might be a mobile interstitial). In (d), the sample is intrinsic at high temperatures, but A₁ and p dominate charge balance, and the sample is p-type as a result of A₁ freezing in at high temperatures, resulting in large supersaturation.

In Fig. 1(a), we illustrated neutral defects which behave independently according to their own thermodynamics and kinetics. Charged defects interact via the charge balance constraint, and thus the results for any charged defect will vary depending on the behaviors of all others and of band carriers. Fig. 1(c) illustrates the concentrations vs. temperature assuming EQ for a hypothetical semiconductor with one fast-diffusing donor-like (+/0) defect D, one acceptor-like (0/−) defect A₁ with high migration barrier, and a 2nd acceptor A₂ with low migration barriers. Their formation energies are such that, under EQ, intrinsic carriers dominate charge balance below ∼2/3 of the max T shown, while the low formation energy of D results in charge balance between D and electrons at high temperatures. However, when considering quenching from intermediate temperatures, we must examine the ordering of defect concentrations.

From Fig. 1(c), one could predict that, if the defect system froze in at high temperatures, n-type doping would result, while intrinsic doping would result if the defect system froze in at lower temperatures (since [D], [A1], and [A2] would all be below about 2/3 of the max T scale shown. In contrast, Fig. 1(d) depicts an SQ calculation for the same semiconductor and defects, which predicts p-type doping for finite cooling rate because acceptor A₁ with a high migration barrier freezes in at high temperature. Charge balance becomes dominated by the very supersaturated A₁ acceptors and holes for all lower temperatures, while A₂ and D keep decreasing in number to lower temperatures because of their lower kinetic barriers, and so might only be detectable using capacitive techniques like DLTS. For clarity: in this example, the formation and migration energies follow the ordering 0 < E_form,D < E_form,A1 < E_form,A2 and 0 < E_mig,D < E_mig,A2 < E_mig,A1.

The central result is that sequential freeze-in of different defects (and charge states) produces room-temperature defect ensembles that are not contained within the limiting EQ or FQ solutions, even when all scenarios share the same initial high-temperature state.²³ In other words, situations will exist in real-world semiconductor processing for which typical defect analysis is completely incapable of predicting the real outcomes. This is another dramatic illustration that examining defect formation enthalpies and phase equilibrium using density functional theory (DFT) effectively at 0 K is insufficient in terms of predicting real-world behaviors.^{3,18,19,24–31} Thus, we assert that SQ is not simply a more complicated analysis that will, in the end, yield marginal additional insight. Rather, combining both defect thermodynamics and kinetics, as in SQ, is critical for accurate predictions.

CdTe is technologically important in bulk crystal form for radiation detectors and as substrates for infrared, but in micron-scale polycrystalline form is also the most commercially-produced thin film photovoltaic technology.⁸ Especially in this context, the native defect system tends to compensate for the intentional dopants. To wit, group-V (gr-V = N, P, As, Sb) acceptors can reach activation ratios (p/[grV]) of 50% up to approximately 1–2 × 10¹⁷ cm⁻³ in melt-grown bulk crystals, yet it is difficult to achieve few-% activation in films grown by physical vapor deposition techniques suitable for high-throughput solar cell manufacturing.^4,8,32–34 Doping into the desired 10¹⁵–10¹⁷ cm⁻³ range using closed space sublimation (CSS) or vapor transport deposition (VTD) may require [As] at the 0.1 at% level.^5–7,35 While the phenomenology of low gr-V doping activation is well-established, the microscopic causes remain unproven. Because native defects or the dopant itself in configurations other than substitutional acceptors might compensate holes, as opposed to some detectable other impurity, hypothesis testing has focused on computation.

Around a decade ago, the best available density functional theory computations using hybrid functionals concluded that P_Te, As_Te, and Sb_Te acceptors were unstable with respect to a distorted AX donor-like configuration, leading to dopant self-compensation.^36–38 Additionally, various complexes were also proposed.^8,26,39,40 More recent computations taking into account CdTe's large spin-orbit coupling in the valence band predict gr-V on Te acceptor ionization levels in better agreement with experiments, especially for Sb.^4,30,41 Notably, these works predict that AX behavior should not occur, and instead propose donor-like Cd interstitials (Cd_i), defect complexes involving As (As_Te-Cd_i) as key contributors to compensation^{29,30,42–44} whose populations depend sensitively on thermal history and cooling conditions. As_Cd “antisite” donors may also appear in Te-rich cases as discussed herein.³⁰ In this paper, we explore the implications of combined defect thermodynamics and kinetics on As acceptor doping in CdTe.

In order to allow process- and sample-dependent modelling of defect concentrations without full-blown PDE simulations, we introduce sequential quenching (SQ) as part of KROGER.^3,45 KROGER is a computational framework for defect concentration calculations focusing on merging computed formation energetics with real-world thermochemistry and thermochemical constraints specific to processing methods and finite-temperature effects on defect free energies. In SQ, defect species freeze at characteristic temperatures determined by their migration barriers, the cooling rate, and the characteristic distance, while more mobile defects continue to equilibrate to lower temperatures. This approach captures essential non-equilibrium physics without requiring fully spatially resolved reaction-diffusion simulations, enabling efficient estimation of defect concentrations in finite-sized samples. Especially in the past decade, as the computational cost of computing defect migration barriers from DFT has come down and accuracy has increased, many other works have, in their postprocessing analysis, characteristic temperatures at which individual defects’ concentrations become frozen-in based on attempt frequency or diffusion length considerations. Our SQ method generalizes and automates these calculations and allows for cross-coupling of defect diffusivities and concentrations – e.g., diffusion of substitutional impurities will depend on the vacancy concentration and impurities will add “drag” to the vacancies’ diffusion.

Our SQ framework implements analysis based on the ratio of the characteristic diffusion length to sample dimensions or inter-sink/source distances. It can thus robustly estimate how defect behavior will change with radius in bulk boules and be different between large crystals, thin films, and polycrystals. The condition of charge neutrality, which arises because the self-energy of a charge density of ppm defects dwarfs defect formation energies for volumes larger than ∼1 µm³, mutually couples charged defects. If any defect's concentration freezes in, it will affect the resulting defect ensemble. In this work, we apply SQ to As-doped CdTe as a representative and technologically-relevant case study, demonstrating how cooling rate, characteristic distance, and defect migration kinetics jointly govern compensation and gr-V activation.

Remarkably, the SQ physically-transparent defect analysis framework explains the long-baffling, experimental observations that high p-type doping and activation are straightforward in bulk (>1 cm) crystals, but that polycrystalline thin films with thickness and grain size of a few µm have low activation and hole density. Optimized thermochemical treatments can achieve high activation in epitaxial thin films.^35,46 Especially if Cd_i is taken to be the dominant source of compensation, these observations are difficult to rationalize: wouldn’t it be easier for Cd_i, which (like many other cations: Cu, Au, Na, Li, etc.) remains mobile at room temperature in the large CdTe lattice, to escape to interfaces and grain boundaries in polycrystalline thin films? The ability of the relatively-computationally-light KROGER with SQ to estimate defect concentrations and predict such differences in defect behaviors as a function of sample geometry and microstructure makes it a unique and powerful tool for defect analysis in specific processing conditions.

Calculation methods

Sequential Quenching (SQ) is a calculation procedure enabled by the ability in KROGER⁴⁵ to dynamically declare individual charge states, defects, or elements as “frozen” – meaning that their concentrations are fixed for all lower temperatures – or “open” meaning that their numbers continue to change in response to temperature, Fermi level, and relevant element chemical potentials.³ SQ allows us to account for diffusion-limited kinetics to model defect concentrations resulting from finite cooling rates, rather than the limits of full equilibrium (no kinetic limitations) or full quenching (no diffusion or reaction). The concepts behind the SQ calculation are simple. A set of defects and charge states is defined with formation energies, migration energies, diffusion prefactors, charges, and other relevant attributes. Another attribute for charge states lists the native defect(s)/charge state(s) X_i that mediate its diffusion.¹⁵ At the highest temperature, T_max, the concentrations of all defects are specified as initial conditions. Herein, we explore both the cases of specified as chemical potential µ_As and specified, fixed total As concentration, but equilibrating all defects through charge balance. However, other initial conditions are possible: for example, concentrations of dopants and Frenkel pairs could be specified in order to predict defect evolution during annealing following ion implantation.

Starting at T_max, SQ proceeds down a list of temperatures usually ending at T_min = 300 K, and at each temperature, the diffusivity D^q_j of each charge state q of defect j is computed from the current concentration(s) of its mediating native defects and its Arrhenius law:

in which D^q_j,oo is the diffusion prefactor for the mechanism mediated by defect or charge state X_i (or X_l, X_m…) and E_m is the corresponding migration energy. If n diffusion mechanisms are possible simultaneously for a defect, they are simply added such that the total diffusion constant for defect X^q_j is D^j,tot = …+ D^q−1_j(X_l) + D^q−1_j(X_m) + D^q−1_j (X_l) + D^q−1_j(X_m) + …

A defect is still able to change its concentration As long as at least one diffusion mechanism is operating fast enough to allow it to diffuse from/to sources/sinks.

We define T_freeze for a defect/charge state as the temperature below which its remaining diffusion length during cooldown is insufficient to communicate with sources/sinks over distance L. All charged defects remain mutually coupled through charge neutrality at every temperature, so freezing one species constrains the entire ensemble. We note that our current treatment does not capture defect reactions such as complex formation, which will be addressed in later work. Similarly, grain boundaries and interfaces are not treated explicitly with distinct migration energetics in the present formulation, but could possibly be captured by users by elaborating this SQ approach by considering multi-step intragrain and grain boundary diffusion paths with their characteristic lengths L.

Substitutional defects typically have small D_o and large E_m, while interstitials have the opposite trends. Especially for continuous cooling, the exponential dependence on E_m dominates the freeze-in behavior, while D_o shifts freeze-in temperatures only weakly. Herein, we adopt D_oo = 0.25 cm² s⁻¹ for all defects, which allows the best agreement with the activation data of Nagaoka et al.⁶ This would correspond to assuming the vibrational frequency Γ ≈ 6 × 10¹³ Hz in the microscopic relation D_o = Γa_o², which is perhaps slightly high. Ideally, this prefactor would be determined individually for each diffusion mechanism, as would the global adjustable parameter f denoting the fraction of bandgap temperature dependence occurring in the conduction band. Because D_o shifts T_freeze logarithmically, our qualitative conclusions are robust to order-of-magnitude variation in D_o. We verified that using D_o/10 and 10D_o changes absolute freeze-in temperatures but preserves the key trends in activation versus γ and L. Future work can refine mechanism-specific prefactors from phonon/attempt-frequency analysis.

In our implementation of SQ, we assume that a defect freezes-in at the temperature T_freeze, at which its root-mean-square diffusion length for the remainder of the cooling process falls below the characteristic distance between sources/sinks. Freeze-in means that it could no longer equilibrate its numbers via diffusion to/from sinks for all lower temperatures. Interestingly, closed-form solutions are possible for the integrated diffusion length of Arrhenius diffusion during continuous cooling for constant cooling rate (linear) and Newtonian (exponential) T(t) trajectories.^47,48 For linear cooling T(t) = T_max −γt with t ∈ [0,(T_max − T_min)/γ] the solution is exact, while the integral for T(t) = T_min + (T_max − T_min)exp(−γt/(T_max − T_min)) does not converge on t ∈ [0, ∞] but the integral using T(t) ≈ T_max exp(−γt/T_max) provides a very good approximation (D(T(t)) and thus the integral x²(D(T(t))) are well approximated, not T(t) itself). For the linear case, carrying out the integration from t = 0 to t = (T_current − T_min)/γ gives the integrated diffusion length of a defect during cooling from the current temperature T_current to T_min as:

in which and . For the case of exponential free cooling, the approximation mentioned yields:The exponential integral Ei(z) herein is defined as, and Γ(a,z) denotes the upper incomplete Gamma function . Caution is advised in implementing these results in different programming languages because of differing definitions of special functions. In Mathematica Ei[z] = ExpintegralEi[z], in SciPy Ei[z] = expi(z), however, in Matlab we utilize Ei(z) = −real(expint(−z)), which is valid for z real and positive. Also, in Matlab the definitions of the incomplete Gamma functions differ from Mathematica, but luckily, Γ(a,z) is equal to -real(expint(−z)) in Matlab for z both real and positive.

Throughout, γ ≡ −dT/dt denotes the cooling rate, and L is the characteristic transport length (taken as half the smallest relevant dimension, e.g., half grain size or half film thickness). A defect (or charge state) is labeled open if its concentration continues to re-equilibrate at the current temperature, and frozen if its concentration is held fixed for all lower temperatures. When a defect becomes frozen, its total concentration remains fixed for all lower temperatures, however, the admixture of its charge states is self-consistently recalculated at each temperature according to a Gibbs distribution in the same way typical for FQ calculations and as described in Arnab et al.³ where the concentration for each charge state is

In the case of no coupling between the defect of interest and any others, these equations could simply be solved for T_freeze. This situation could be found, for example, for interstitials that are neutral or that never become one of the defects dominating charge balance. However, for a general system of defects coupled via charge balance and diffusion mediation, these represent the current best estimates for its maximum diffusion length, assuming that the concentration of any diffusion-mediating defects never increases at temperatures below T_current. In each SQ calculation, we sweep over multiple cooling rates and characteristic distances, thus, we step the calculation from T_max to T_min in small steps (e.g., 1 K for rigor) and at each temperature, check which, if any, defects freeze-in for permutations of these two parameters.

Herein, we use the example of CdTe and explore the case where interfaces and grain boundaries are assumed to be the dominant sources/sinks providing, in the absence of 2nd phase precipitation, the dominant pathway for elimination of defects such as Cd_i during cooling with Cd-rich conditions (or Te_i in Te-rich). The characteristic source/sink distances are taken as ½ of the grain size, film thickness, or bulk crystal's smallest dimension. The incorporation of quasichemical rates for reactions like V_Cd + Cd_i = ∅ or As_Te⁻ + Cd_i²⁺ = (As_Te-Cd_i)⁺ is beyond the scope of the present work, but we account for As_Te-Cd_i complexes by treating them as a separate defect with a large migration energy.

Each SQ calculation is begun at a specified high temperature, with the user specifying initial concentrations, whether each element should be treated as fixed-total number or fixed-chemical potential, and whether any defects or charge states should be frozen for the entire calculation. The case of the dopant concentration [As] being specified (as in the green defect labeled “fixed conc.” in Fig. 1(a)) is appropriate for liquid phase crystal growth in which the dopant concentration in the well-mixed liquid is below the solubility limit for the maximum temperature where As might diffuse and is locked-in at time of crystallization. The case in which µ_As is set by equilibrium with Cd₃As₂, amorphous As_Te, or As₂Te₃ is appropriate for cases of high high-doping in which As may precipitate out as 2nd phases.^4,6,49

For clarity, we separate the SQ inputs into (i) thermochemical boundary conditions and finite-temperature band-edge assumptions that determine µ_i(T) and the charge-neutral Fermi level, and (ii) a defect energetics database (formation energies, transition levels, and migration barriers) that controls equilibrium populations and kinetic freeze-in. We describe (i) first, followed by (ii). As in our prior work on defects in β-Ga₂O₃, we compute chemical potentials from the best-available thermochemical data optimized for the host Cd-Te binary phase system. We define scenarios of CdTe in equilibrium with elemental Cd or Te, corresponding to gedanken experiments with CdTe touching the element or in equilibrium with its equilibrium vapor pressure. We define the 0 of our electron energy scale as the valence band energy at 0 K, and use a free parameter f to parameterize the fraction of bandgap temperature dependence E_g(T) occurring in the conduction band. Unfortunately, we do not have access to the vibrational free energies of each charge state of each defect, nor are the temperature dependence of charge transition levels of the defects available experimentally at material processing temperatures.⁵⁰ Thus, we effectively assume that the charge transition levels of all defects remain fixed and that the band edges move versus temperature. In some sense then, the f parameter is used herein to represent how the temperature dependencies of the dominant defect(s)’ charge transition levels move relative to the band edges. In future calculations, incorporating vibrational free energies of all defects, f should be set strictly according to the absolute band edge energies. We also utilize the same approximation for the effects of vibrational entropy on all the charge states of each defect in terms of average mode counting. We define the representative mode frequency for the host from its Debye temperature, and then the −TS_vib term for 3 quantum harmonic oscillator modes is added for each atom added (e.g., for interstitials) and 3 subtracted for each atom removed (e.g., for vacancies). This results in no change for substitutionals and antisites, which may explain why we find that As_Cd donor-like defects may be significant, whereas other analyses have implicated Cd_i only. We do note that Nagaoka reported XFH data suggesting that As_Cd may exist for heavily As-doped samples, although it remains somewhat unclear whether that result might have been related to defect clustering or Cd₃As₂ precipitation.⁵¹ Additionally, in unpublished NMR studies, Pougue, Nagaoka, Scarpulla, and Rockett found solid-state NMR signatures in similar samples that suggested unique atomic configurations. For these reasons, we present our conclusions herein from the best SQ analysis possible at this time, but anticipate that the doping behavior of gr-V's may still reveal further surprises.

At each temperature, defect concentrations are obtained by solving for the Fermi level under global charge neutrality, including contributions from charged defects and free carriers. When no elemental constraints are imposed, a one-dimensional root-finding procedure is used. When fixed element concentrations are specified, a multi-variable optimization is used to simultaneously satisfy charge neutrality and elemental balance. For constrained calculations, a hybrid optimization approach is used, consisting of a particle-swarm search followed by simplex refinement. The particle swarm uses population sizes of approximately 300–15000 and objective limits of 10¹⁰ (initial search) and 10⁸ (refined search), while the final simplex step uses tolerances of 10⁻⁵ for both variable and function convergence. For full-quench calculations, the Fermi level is scanned with a step size of k_BT/5. These solver settings are user-defined and are selected such that residual charge and element-balance errors remain negligible compared to the defect and carrier concentrations of interest. Convergence is verified by ensuring that further tightening of these parameters does not change the Fermi level or dominant defect concentrations. Diffusion is modeled using an Arrhenius form (eqn (1)), where migration barriers E_m are taken from first-principles calculations and are assumed temperature-independent. In the present work, diffusion prefactors D₀ are assigned uniform values across defects to isolate the role of migration barriers and cooling conditions. When multiple diffusion mechanisms exist for a given defect, their contributions are summed. Further details of the numerical implementation and thermodynamic framework are provided in our previous paper on KROGER.³

First-principles defect energetics were obtained from density functional theory (DFT)^52,53 calculations performed with the projector augmented-wave (PAW) method^54,55 as implemented in VASP.^56,57 To accurately describe the CdTe band gap and defect energetics, we employed the screened hybrid functional of Heyd–Scuseria–Ernzerhof (HSE),^58,59 including spin–orbit coupling (SOC) in the total energies. Following our established CdTe defect workflow, the fraction of Hartree–Fock exchange was set to 33%, and the plane-wave energy cutoff was 300 eV. With this setup, the calculated equilibrium lattice constant of CdTe is a = 6.545 Å and the band gap is E_g = 1.504 eV, consistent with experimental values.⁶⁰

Point defects and defect complexes were modeled in 216-atom supercells. All structures were fully relaxed until residual forces were below 0.02 eV Å⁻¹. Brillouin-zone integrations were performed at the Γ point; tests with off-Γ meshes do not alter the qualitative conclusions for the defect trends considered here.

Defect formation energies were evaluated using the standard formalism.⁶¹ For a defect (or defect complex) D in charge state q,

where E_tot(D^q) and E_tot(bulk) are total energies of the defective and pristine supercells, n_i counts atoms of species i added (n_i > 0) or removed (n_i < 0), and µ_i are the corresponding chemical potentials constrained by CdTe phase stability. The Fermi level ε_F is referenced to the valence-band maximum (VBM) E_v of bulk CdTe and varies within the band gap. Finite-size corrections E^q_corr account for spurious electrostatic interactions and potential alignment in charged supercells;⁶¹ electrostatic potential alignment between charge states was performed by comparing the average electrostatic potential in a region far from the defect. Thermodynamic charge-transition levels ε(q/q′) were obtained from the crossings of E_f (D^q) and E_f (D^q′).

Migration barriers were determined using the climbing nudged elastic band (cNEB) method^62,63 with at least six intermediate images between the initial and final configurations. To reduce the computational cost while maintaining consistency with hybrid-functional energetics, minimum-energy paths were first relaxed at the GGA level using PBEsol, followed by single-point HSE total energy calculations for initial image and saddle point to obtain barriers consistent with the HSE defect formation energetics.

Results & discussion

Fig. 2 compares defect concentrations as a function of temperature obtained from sequential quenching (SQ) and equilibrium (EQ) calculations assuming equilibrium with elemental Cd. Fig. 2(a) shows the calculations performed with a fixed [As] concentration of 10¹⁷ cm⁻³ – corresponding exactly to the Cd-solvent THM growth of Nagaoka and Scarpulla^4–6 – while Fig. 2(b) shows the calculations in which µ_As is determined by equilibrium with Cd₃As₂. Both SQ scenarios are calculated using a characteristic distance of 1 cm and a cooling rate of 0.05 K s⁻¹, which correspond to the growth parameters adopted from Nagaoka et al.⁶

Fig. 2 Defect concentrations as a function of temperature for CdTe equilibrium with Cd, comparing sequential quenching (SQ) for a characteristic distance of 1 µm & γ = 0.05 K s⁻¹ and equilibrium (EQ) calculations. Two doping scenarios are considered: variable µ_As determined (a) by fixing As concentration of 10¹⁷ cm⁻³, and (b) by equilibrium with the secondary phase Cd₃As₂. Solid lines and filled symbols correspond to equilibrium (EQ) calculations, while dashed lines and open symbols correspond to sequential quenching (SQ) calculations.

In equilibrium calculations (EQ), defect concentrations are allowed to adjust continuously at every temperature, leading to non-monotonic concentration trajectories during cooldown. For example, under equilibrium with Cd conditions, the As_Te-Cd_i complex remains at concentrations of ∼5 × 10¹⁵–10¹⁶ cm⁻³ between 1200 K and ∼1000 K, but increases to ∼4 × 10¹⁶ cm⁻³ at ∼700 K till cooling down. This behavior reflects the implicit assumption of equilibrium modeling that defect populations can fully reconfigure even at temperatures where long-range diffusion would be strongly kinetically limited.

In contrast, in sequential quenching (SQ) framework, defects with sufficiently high migration barriers cease to respond directly to further cooling, while more mobile species may continue to evolve. As a result, changes in defect concentrations at lower temperatures can still occur through shifts in charge balance and chemical potentials driven by mobile defects, particularly Cd_i under equilibrium with Cd conditions. This behavior is especially evident under 2nd phase-limited conditions, where variations in Cd_i concentration at low temperature can induce corresponding changes in V_Cd populations and defect complexes, even though less mobile defects remain frozen. SQ therefore constrains the active degrees of freedom during cooldown rather than enforcing complete kinetic arrest of all defect populations.

Notably, the qualitative SQ trends are consistent with prior kinetic simulation approaches developed for dopant activation and defect interactions in CdTe (e.g., kinetic models beyond purely thermodynamic treatments).^18,26 While SQ does not solve spatially resolved reaction–diffusion–Poisson equations, it reproduces the same physically expected directionality: activation decreases when mobile donor-like species remain available to compensate during cooldown, and increases when those species are eliminated or kinetically arrested at higher temperature.

The temperature ranges over which individual defects deviate from equilibrium further reflect their relative mobilities. Fast-diffusing species such as Cd_i remain mobile to lower temperatures, whereas defect complexes such as As_Te-Cd_i, having a large migration barrier, freeze in at higher temperatures. Hierarchy of freeze-in behavior is a central feature of SQ and cannot be captured within either equilibrium or full-quench limits.

These differences in defect evolution result in pronounced variations in predicted carrier concentrations. In the fixed-[As] case shown in Fig. 2(a), SQ predicts hole concentrations on the order of ∼10¹³ cm⁻³ at room temperature, whereas EQ predicts a much lower value of approximately 10⁸ cm⁻³. This discrepancy arises primarily from differences in the population of Cd_i donors, whose concentration is strongly influenced by the freeze-in of As_Te-Cd_i complexes at elevated temperatures.

Similar qualitative trends are observed in Fig. 2(b), where the µ_As is constrained by the presence of the Cd₃As₂ second phase. EQ predicts that As-related defects diminish rapidly with decreasing temperature and have negligible influence on room-temperature carrier concentrations. In contrast, SQ freezes As-related defects at higher temperatures, where AsTe-Cd_i dominate compensation. Nevertheless, the presence of a 2nd phase constraint imposes a strong thermodynamic limit on As incorporation, and SQ alone is insufficient to produce p-type behavior under these conditions.

To clarify the behavior of V_Cd in Fig. 2(b), its increase with decreasing temperature in SQ does not arise from continued defect formation at low temperature. Instead, it reflects how charge neutrality is enforced when only a subset of defects remains mobile. At each temperature, the Fermi level is determined self-consistently, and both band carriers and defect charge states contribute to charge balance. However, as cooling proceeds, many defects (including dominant negatively charged species such as As⁻_Te and As_Te-Cd_i) freeze at higher temperatures and can no longer adjust their total concentrations. In contrast, more mobile donor-like species such as Cd²⁺_i remain open to lower temperatures and continue to evolve under the imposed thermodynamic constraints. As Cd²⁺_i adjusts, the charge neutrality condition must still be satisfied. While electrons and holes also contribute to charge balance, their concentration is constrained by the position of the Fermi level and available states, and cannot always fully compensate the evolving charges. In this regime, V²⁻_Cd, which remains open over a broader temperature range, provides the necessary compensating charge. Therefore, the increase in V_cd reflects the fact that, at that stage of cooldown, it is one of the few remaining defects able to respond to charge neutrality constraints, rather than an indication of unrestricted defect formation. This behavior arises naturally from the sequential freeze-in of defects with different mobilities.

Fig. 3(a) and (b) illustrates SQ-predictions of radial defect and carrier profiles for CdTe under two thermodynamic boundary conditions; (a) CdTe in equilibrium with elemental Cd (equilibrium with Cd) and (b) CdTe in equilibrium with elemental Te (equilibrium with Te). For each boundary condition, results are shown for fast (0.05 K s⁻¹) and slow (10⁻⁴ K s⁻¹) cooling rates. The total As concentration is held fixed at 10¹⁷ cm⁻³ in all cases to examine the influence of As activation on characteristic distance.

Fig. 3 Sequential quenching (SQ) predictions of radial defect and carrier profiles in CdTe under two cooling rates and two thermodynamic boundary conditions. Panels (a) and (b) show SQ results for equilibrium with Cd and Te conditions, respectively, for cooling rates of 0.05 K s⁻¹ (fast cooling; dashed lines and open circles) and 10⁻⁴ K s⁻¹ (slow cooling; solid lines and filled circles). Faster cooling shifts defect freeze-in to smaller characteristic distances, leading to enhanced donor compensation and the onset of n-type behavior at reduced characteristic distances under equilibrium with Cd conditions. Panels (c) and (d) show the corresponding SQ results after removal of all mobile Cd_i following cooldown, representing the maximum achievable dopant activation when interstitial donors diffuse to sinks during post-growth processing. Under equilibrium with Cd conditions, substantial changes in carrier concentrations are observed due to the dominant compensating role of Cd_i, whereas under touching Te conditions, the profiles remain unchanged, reflecting negligible contribution of Cd_i in this regime. In all cases, the total As concentration is fixed at 10¹⁷ cm⁻³.

In both boundary conditions, faster cooling shifts defect freeze-in toward smaller characteristic distances, resulting in n-type behavior in equilibrium with Cd scenario beginning at lower characteristic distances (characteristic distance starting n-type behavior changes to ∼2 cm from ∼50 cm). Increasing characteristic distance further promotes early freeze-in, as longer diffusion lengths are required for defects to reach sinks during cooldown.

Under equilibrium with Cd-rich conditions, Cd_i plays the dominant role in compensation. For fast cooling and large characteristic distances, high concentrations of Cd_i freeze-in (∼10¹⁷ cm⁻³), driving strong donor compensation and n-type behavior. Slower cooling allows Cd_i to remain mobile to lower temperatures, reducing its frozen-in concentration (∼4 × 10¹⁶ cm⁻³) and thereby decreasing compensation. These results highlight the strong sensitivity of electrical behavior to both cooling rate and sample geometry under Cd-rich boundary conditions.

In contrast, under equilibrium with Te-rich conditions, compensation is governed primarily by As_Cd, which dominates at small characteristic distances. At larger characteristic distances, the concentration of As_Cd remains relatively low concentration of 10¹⁴ cm⁻³, resulting in reduced compensation. At small characteristic distances, however, the As_Cd concentration increases to 5 × 10¹⁶ cm⁻³ and becomes As_Cd = As_Te, leading to strong mutual compensation and a corresponding reduction in hole concentration.

Fig. 3(c) and (d) show the corresponding SQ results after modifying the thermodynamic environment following cooldown, allowing mobile Cd_i defects to diffuse out of the sample at low temperatures. This scenario represents an idealized upper bound on dopant activation that may be achieved if free interstitial donors escape to surfaces or grain-boundary sinks during post-growth processing, such as annealing or exposure to oxidizing environments.

Under equilibrium with Cd-rich conditions, removal of Cd_i leads to a substantial reduction in compensation, causing the hole concentration to approach the As_Te concentration and yielding activation levels approaching ∼98%. Residual As_Te-Cd_i (4–8 × 10¹⁶ cm⁻³) complexes remain frozen in the bulk and maintain a small degree of compensation. In contrast, under equilibrium with Te-rich conditions, the compensation behavior remains unchanged following Cd_i removal, as As_Cd is the dominant compensating defect. Consequently, the defect and carrier profiles before and after Cd_i removal are nearly identical in this regime.

Fig. 4 summarizes the predicted As activation in CdTe under different thermodynamic and kinetic modeling scenarios. Fig. 4(a) shows results obtained using a fixed As concentration, while Fig. 4(b) shows calculations in which the µ_As is constrained by equilibrium with a second phase (Cd₃As₂ and AsTe for the touching Cd and Te conditions, respectively).

Fig. 4 Predicted As activation in CdTe under different thermodynamic and kinetic modeling scenarios. Panels (a) and (b) show results obtained using a µ_As constrained by a fixed As concentration and by a second phase, respectively. Blue-shaded curves correspond to equilibrium with Cd boundary conditions, while red-shaded curves correspond to equilibrium with Te conditions. Equilibrium with Te conditions consistently yield higher As activation due to the absence of strong interstitial donor compensation, whereas activation under equilibrium with Cd conditions remains significantly limited, with a maximum activation of approximately 30%. For comparison, experimental activation data reported by Nagaoka et al.⁶ for bulk material are included in both panels with olive circles, demonstrating qualitative agreement with the trends predicted by sequential quenching. 2% activation reported for thin film CdTe by Mallick et al.⁶⁴ is also included in (a).

In the fixed [As] case shown in Fig. 4(a), sequential quenching (SQ) calculations under equilibrium with Cd-rich conditions yield a maximum activation of approximately 30% for an As concentration of 10¹⁶–10¹⁷ cm⁻³ at a characteristic distance of 0.5–1 cm under fast cooling. Activation varies systematically with characteristic distance, reflecting the sensitivity of defect freeze-in to diffusion length. For a given characteristic distance, fast cooling generally results in higher activation than slow cooling due to reduced donor equilibration at lower temperatures. In all equilibrium with Cd-rich cases, As activation decreases markedly with increasing dopant concentration. This trend closely resembles experimental activation data reported by Nagaoka et al.,⁶ demonstrating qualitative agreement between SQ predictions and measured behavior in CdTe crystals.

We also include one citable activation data reported by Mallick et al. for CdSeTe based thin film solar cell absorbers (shown as the orange triangle). Anecdotally, we know from wide experience of many groups involved in CdTe thin film solar research and off-the-record communications from First Solar that activation of P, As, and Sb is typically very low, however, there are few published quantifications. It is important to note that this result from Mallick et al. is for films subject to significantly more complexity than the bulk crystal growth results. These were for Se containing alloyed layers after CdCl₂ treatment and also subject to potential oxygen incorporation from the CdCl₂ environment and possibly from reactions with the front TCO layers (SIMS data we have been shown by FSLR, but is not published, has found correlations of O and As). Lastly, it is known that grain boundaries and interfaces of Cd(Se)Te usually exhibit some degree of Fermi level pinning, which means that depletion widths and built-in voltages exist at GBs and surfaces. These would provide a confining potential to charged defects, further complicating the question of whether or not defects could reach microstructural sinks/sources. Additionally, our SQ, as presented herein, only considers mass transport to/from sinks/sources as the only mechanism for eliminating excess defects. In reality (and in planned future work) complexation kinetic reactions between defects, and even precipitation of defect clusters and 2nd phases would also need to be included. In the bulk crystal growth cases, these other elements were intentionally avoided.

For these reasons, we do not believe that our SQ modelling in this paper, which included none of these possible effects, should be quantitatively compared to the Mallick et al. data point, which would be subject to all the listed complications (and possibly more). In the big picture, we focus on the fact that the SQ framework captures a physically-expected effect that to remain equilibrated, defects must either transport to sinks/sources or react to form some other species. This captures common observations like the fact that divacancies and voids may form in the center of large boules, whereas the regions closer to the surface are more perfect crystals (e.g., early experience of vacancy condensation in Ge boules grown at 5N Plus Semiconductor). No simple, non-PDE-solving defect modelling framework has provided this previously.

In contrast, equilibrium with Te-rich conditions yield substantially higher activation over a wide concentration range. Near-unity activation is obtained up to approximately 10¹⁷ cm⁻³ [As] doping at large characteristic distances, followed by a decrease in activation at higher dopant concentrations. At large characteristic distances, activation decreases to approximately 2% for slow cooling and 15% for fast cooling at a high As concentration of 10¹⁹ cm⁻³. At small characteristic distances, activation remains low, with maximum values of only a few percent even under fast cooling conditions.

The enhanced activation under equilibrium with Te-rich conditions arises from the absence of Cd_i donors and the relatively low concentration of the dominant compensating defect, As_Cd. At higher As concentrations, activation decreases due to the stabilization of compensating defects, including Te_Cd-related defect complexes, which limit the number of electrically active acceptors.

Fig. 4(b) shows activation behavior when the µ_As is constrained by equilibrium with a second phase. In this case, equilibrium calculations under equilibrium with Cd conditions exhibit trends similar to experimental data at elevated temperatures. Sequential quenching results under second-phase-limited conditions predict nearly complete activation under equilibrium with Te-rich conditions, whereas activation under equilibrium with Cd-rich conditions remains extremely low for both fast and slow cooling at a characteristic distance of 1 cm. at smaller characteristic distances, activation is strongly suppressed in both boundary conditions, indicating that second-phase constraints impose a dominant thermodynamic limit that cannot be overcome by kinetic freeze-in alone.

The SQ results highlight the central role of Cd_i in controlling compensation under Cd-rich conditions. While Cd_i is highly mobile in the CdTe lattice, its ability to diffuse out of the bulk depends not only on diffusion kinetics but also on the electrostatic landscape within the material.

In polycrystalline CdTe, built-in electric fields associated with grain boundaries and contact depletion regions may significantly influence Cd_i transport. In the dark, depletion regions near grain boundaries and interfaces can create potential barriers that inhibit long-range diffusion of positively charged Cd_i²⁺, leading to accumulation within grain interiors. Under illumination, partial flattening of band bending may temporarily increase the accessible volume for Cd_i, allowing redistribution within the CdTe lattice. Notably, the binding energy of the As_Te– Cd_i complex (0.7 eV)³⁰ exceeds typical grain-boundary electrostatic barriers (∼0.1–0.2 eV),⁶⁵ indicating that neither equilibrium band bending nor its partial flattening under illumination is sufficient to dissociate the complex. Consequently, only the population of free Cd_i, formed at high temperatures, is expected to participate in redistribution under operating conditions.

Two limiting scenarios may therefore exist for the free Cd_i: (1) Cd_i redistributes reversibly within the CdTe grains with increased volume in response to changing electrostatic conditions, or (2) Cd_i reaches grain boundaries and becomes chemically stabilized or annihilated through reactions with impurities or ambient species. The SQ framework captures the first one through characteristic distances and the second one through the removal of Cd_i. This could lead to a pathway for long-term defect evolution that may contribute to light-induced metastability and “wakeup” phenomena observed experimentally.

To our knowledge, no direct experimental evidence of the symmetry-lowering distortion of gr-V acceptors resulting in the AX configuration has been presented to date.^66,67 In both well-known papers claiming to have observed AX centers, the assignment of AX was, in hindsight, based on less-than-definitive evidence. In Ref.,⁷ the predicted transition energy for the AX center matched an observed 400 meV hole thermal emission energy not present in As-free crystals. The certainty with which this assignment was stated is, in our opinion, not warranted in crystals that were in the same paper reported to contain 2nd phase Cd₃As₂ precipitates visible with infrared microscopy. Additionally, it is now understood that the computed valence band was misaligned, thus all of the computed transition levels would need to be corrected – meaning that the 400 meV signal would correspond to a different defect. In ref. 68 the evidence for AX instability was the observation of persistent photoconductivity (PPC). Two quenching-induced doping decay modes were also detected, suggesting multiple defect processes. Photoexcitation of the AX centers to the symmetric substitutional donor would be expected to exhibit PPC, but are not the only possible origins. All that is necessary is that photoexcitation creates a long-lived metastable defect + carrier configuration, which could be a local relaxation as for DX or AX, or a dissociation and reconfiguration of defects or complexes. For example, PPC in BaTiO₃ and KTaO₃ or photodarkening of Cu-doped β-Ga₂O₃ are all posited to be caused by hydrogen reconfiguring its location.^69–71 It is important to note that forming complexes like (As_Te-Cd_i)⁺ from pre-existing As_Te⁻ and Cd_i²⁺ has no impact on net carrier density – the defect quasichemical reaction must also involve carrier capture/emission (similarly for dissociation, whether thermal or photoexcited).

Conclusions

In this work, we developed a sequential quenching (SQ) framework to model defect evolution under finite cooling rates and applied it to As-doped CdTe. By incorporating defect-specific migration barriers, cooling rate, and characteristic distance to effective sinks, SQ provides a physically grounded description of non-equilibrium defect populations that extends beyond the limiting assumptions of equilibrium and full-quench models. This approach enables the prediction of spatially varying defect and carrier concentrations under experimentally relevant processing conditions.

Application of SQ to CdTe highlights the central role of kinetic limitations in determining dopant activation. Fast-diffusing donor defects, particularly Cd_i remain mobile to low temperatures and can freeze in at large characteristic distances, leading to strong compensation and n-type behavior under equilibrium with Cd-rich conditions. In contrast, slower cooling and reduced characteristic distances limit donor retention and promote higher p-type activation. These results demonstrate that dopant activation in CdTe is governed not only by defect-formation energetics but also by the time-temperature trajectory and the length scales over which mobile defects can diffuse during cooldown.

The analysis further indicates that post-growth redistribution or removal of mobile Cd_i defines an upper bound on achievable activation, while residual defect complexes impose a fundamental compensation limit. In this context, the presence of mobile interstitial donors following growth provides a natural pathway for metastable behavior, where changes in electrostatic conditions, such as those induced by illumination, may enable redistribution of Cd_i without requiring changes in defect formation. Such processes may contribute to experimentally observed “wake-up” phenomena in CdTe, in which carrier concentrations evolve under post-processing conditions.

Taken together, these results establish sequential quenching as a practical and predictive framework for interpreting non-equilibrium defect populations and their evolution under changing thermal and electrostatic environments. More broadly, SQ provides a general methodology for connecting cooling conditions, sample geometry, and defect kinetics to dopant activation and metastability in CdTe and related semiconductor materials.

Conflicts of interest

There are no conflicts to declare.

Data availability

The defect properties, including HSE-computed formation energies and migration energies, are available in the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d6tc00828c.

The KROGER code is available at https://github.com/mikescarpulla/KROGER.

Acknowledgements

This work was based [in part] on research sponsored by the U.S. Department of Energy Office of Energy Efficiency and Renewable Energy Solar Energy Technologies Office agreement number 37989 through National Laboratory of the Rockies, operated under Contract No. DE-AC36-08GO28308. IC and AJ were supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Award Number DE-SC0025506 and made use of computational resources from the National Energy Research Scientific Computing Center (NERSC), a Department of Energy Office of Science User Facility, through the NERSC award BES-ERCAP 0034471 (m5002), and the DARWIN computing system at the University of Delaware, which is supported by the NSF Grant No.∼1919839. The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government. The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for U.S. Government purposes.

References

1. K. Harth and R. Stannarius, Front. Phys., 2020, 8, 112.
2. C. G. Van De Walle and J. Neugebauer, J. Appl. Phys, 2004, 95, 3851–3879.
3. K. A. Arnab, M. Stephens, I. Maxfield, C. Lee, E. Ertekin, Y. K. Frodason, J. B. Varley and M. A. Scarpulla, Phys. Chem. Chem. Phys., 2025, 27, 11129–11143.
4. A. Nagaoka, S. K. Swain and A. H. Munshi, IEEE J. Photovolt., 2024, 14, 397–413.
5. A. Nagaoka, K.-B. Han, S. Misra, T. Wilenski, T. D. Sparks and M. A. Scarpulla, J. Cryst. Growth, 2017, 467, 6–11.
6. A. Nagaoka, D. Kuciauskas, J. McCoy and M. A. Scarpulla, Appl. Phys. Lett., 2018, 112, 192101.
7. T. Ablekim, S. K. Swain, W.-J. Yin, K. Zaunbrecher, J. Burst, T. M. Barnes, D. Kuciauskas, S.-H. Wei and K. G. Lynn, Sci. Rep., 2017, 7, 4563.
8. M. A. Scarpulla, B. McCandless, A. B. Phillips, Y. Yan, M. J. Heben, C. Wolden, G. Xiong, W. K. Metzger, D. Mao, D. Krasikov, I. Sankin, S. Grover, A. Munshi, W. Sampath, J. R. Sites, A. Bothwell, D. Albin, M. O. Reese, A. Romeo, M. Nardone, R. Klie, J. M. Walls, T. Fiducia, A. Abbas and S. M. Hayes, Sol. Energy Mater. Sol. Cells, 2023, 255, 112289.
9. D. Krasikov and I. Sankin, J. Mater. Chem. A, 2017, 5, 3503–3513.
10. K. Ahn, Y. K. Ooi, F. Mirkhosravi, J. Gallagher, A. Lintereur, D. Feezell, E. K. Mace and M. A. Scarpulla, J. Appl. Phys, 2021, 129, 245703.
11. R. F. Davis, A. M. Roskowski, E. A. Preble, J. S. Speck, B. Heying, J. A. Freitas, E. R. Glaser and W. E. Carlos, Proc. IEEE, 2002, 90, 993–1005.
12. D. P. Bour, H. F. Chung, W. Götz, L. Romano, B. S. Krusor, D. Hofstetter, S. Rudaz, C. P. Kuo, F. A. Ponce, N. M. Johnson, M. G. Craford and R. D. Bringans, MRS Online Proc. Libr., 1996, 449, 509–518.
13. M. Miyachi, T. Tanaka, Y. Kimura and H. Ota, Appl. Phys. Lett, 1998, 72, 1101–1103.
14. P. Penning, Phillips Res. Repts., 1959, 14, 337–345.
15. N. D. Rock, H. Yang, B. Eisner, A. Levin, A. Bhattacharyya, S. Krishnamoorthy, P. Ranga, M. A. Walker, L. Wang, M. K. Cheng, W. Zhao and M. A. Scarpulla, APL Mater., 2024, 12, 081101.
16. Y. K. Frodason, P. P. Krzyzaniak, L. Vines, J. B. Varley, C. G. Van De Walle and K. M. H. Johansen, APL Mater., 2023, 11, 041121.
17. W. Patrick, E. Hearn, W. Westdorp and A. Bohg, J. Appl. Phys, 1979, 50, 7156–7164.
18. X. Xiang, A. Gehrke, Y. Tong and S. T. Dunham, ACS Appl. Energy Mater., 2025, 8, 1248–1265.
19. I. Sankin and D. Krasikov, Phys. Status Solidi A, 2019, 216, 1800887.
20. P. Calado, I. Gelmetti, B. Hilton, M. Azzouzi, J. Nelson and P. R. F. Barnes, J. Comput. Electron., 2022, 21, 960–991.
21. S. Anand, M. Y. Toriyama, C. Wolverton, S. M. Haile and G. J. Snyder, Acc. Mater. Res., 2022, 3, 685–696.
22. D. Krasikov, A. Knizhnik, B. Potapkin and T. Sommerer, Semicond. Sci. Technol., 2013, 28, 125019.
23. A. G. Squires, S. R. Kavanagh, A. Walsh and D. O. Scanlon, Nat. Rev. Mater., 2026, 1–18.
24. V. Kosyak, N. B. Mortazavi Amiri, A. V. Postnikov and M. A. Scarpulla, J. Appl. Phys, 2013, 114, 124501.
25. V. Kosyak, A. V. Postnikov, J. Scragg, M. A. Scarpulla and C. Platzer-Björkman, J. Appl. Phys, 2017, 122, 035707.
26. D. Krasikov and I. Sankin, Phys. Rev. Mater., 2018, 2, 103803.
27. L. Diaz, H. P. Hjalmarson, J. J. Lutz and P. A. Schultz, arXiv, 2025, preprint, arXiv:2511.02976 DOI:10.48550/arXiv.2511.02976.
28. L. Diaz, H. P. Hjalmarson, J. J. Lutz and P. A. Schultz, arXiv, 2025, preprint, arXiv:2504.15459 DOI:10.48550/arXiv.2504.15459.
29. I. Chatratin, I. Evangelista, B. McCandless, W. Shafarman and A. Janotti, J. Phys. Chem. C, 2025, 129, 1398–1407.
30. I. Chatratin.
31. P. C. Bowes, G. H. Ryu, J. N. Baker, E. C. Dickey and D. L. Irving, J. Am. Ceram. Soc., 2021, 104, 5859–5872.
32. S. D. Sordo, L. Abbene, E. Caroli, A. M. Mancini, A. Zappettini and P. Ubertini, Sensors, 2009, 9, 3491–3526.
33. M. A. Berding, Phys. Rev. B, 1999, 60, 8943–8950.
34. S.-H. Wei and S. B. Zhang, Phys. Rev. B, 2002, 66, 155211.
35. H. Lott, A. T. Mathew, M. R. Young, A. Goldstone, A. D. Rice, M. O. Reese and E. Colegrove, APL Mater., 2025, 13, 071131.
36. K. Biswas and M.-H. Du, Appl. Phys. Lett, 2011, 98, 181913.
37. J.-H. Yang, W.-J. Yin, J.-S. Park, J. Burst, W. K. Metzger, T. Gessert, T. Barnes and S.-H. Wei, J. Appl. Phys, 2015, 118, 025102.
38. E. Colegrove, J.-H. Yang, S. P. Harvey, M. R. Young, J. M. Burst, J. N. Duenow, D. S. Albin, S.-H. Wei and W. K. Metzger, J. Phys. Appl. Phys., 2018, 51, 075102.
39. D. Krasikov and I. Sankin, J. Mater. Chem. A, 2017, 5, 3503–3513.
40. D. Krasikov, D. Guo, S. Demtsu and I. Sankin, Sol. Energy Mater. Sol. Cells, 2021, 224, 111012.
41. B. McCandless, W. Buchanan, G. Sriramagiri, C. Thompson, J. Duenow, D. Albin, S. A. Jensen, J. Moseley, M. Al-Jassim and W. K. Metzger, IEEE J. Photovolt., 2019, 9, 912–917.
42. I. Chatratin, B. Dou, S.-H. Wei and A. Janotti, J. Phys. Chem. Lett., 2023, 14, 273–278.
43. X. Xiang, A. Gehrke, Y. Tong and S. T. Dunham, ACS Appl. Energy Mater., 2025, 8, 1248–1265.
44. S. Kavanagh.
45. GitHub – KROGER: Program for computing full and partial equilibria of point defects in solids, https://github.com/mikescarpulla/KROGER, (accessed September 15, 2024).
46. A. T. Mathew, H. Lott, E. Colegrove, M. R. Young, D. Kuciauskas, C. A. Wolden and M. O. Reese, J. Appl. Phys, 2025, 137, 115702.
47. M. A. Scarpulla, III-Mn-V ferromagnetic semiconductors synthesized by ion implantation and pulsed-laser melting, 2006.
48. H. Foll, Defects in Crystals, https://dtrinkle.matse.illinois.edu/MatSE584/, (accessed January 28, 2026).
49. J. J. McCoy, S. K. Swain, J. R. Sieber, D. R. Diercks, B. P. Gorman and K. G. Lynn, J. Appl. Phys, 2018, 123, 161579.
50. S. Qiao, Y.-N. Wu, X. Yan, B. Monserrat, S.-H. Wei and B. Huang, Phys. Rev. B, 2022, 105, 115201.
51. A. Nagaoka, K. Kimura, A. K. R. Ang, Y. Takabayashi, K. Yoshino, Q. Sun, B. Dou, S.-H. Wei, K. Hayashi and K. Nishioka, J. Am. Chem. Soc., 2023, 145, 9191–9197.
52. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864–B871.
53. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133–A1138.
54. Projector augmented-wave method | Phys. Rev. B, https://journals.aps.org/prb/abstract/10.1103/PhysRevB.50.17953, (accessed March 9, 2026).
55. G. Kresse and D. Joubert, Phys. Rev. B, 1999, 59, 1758–1775.
56. G. Kresse and J. Hafner, Phys. Rev. B, 1994, 49, 14251–14269.
57. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
58. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys, 2003, 118, 8207–8215.
59. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys, 2006, 124, 219906.
60. R. W. Birkmire and E. Eser, Annu. Rev. Mater. Res., 1997, 27, 625–653.
61. C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti and C. G. Van De Walle, Rev. Mod. Phys., 2014, 86, 253–305.
62. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys, 2000, 113, 9901–9904.
63. G. Henkelman and H. Jónsson, J. Chem. Phys, 2000, 113, 9978–9985.
64. R. Mallick, X. Li, C. Reich, X. Shan, W. Zhang, T. Nagle, L. Bok, E. Bicakci, N. Rosenblatt, D. Modi, R. Farshchi, C. Lee, J. Hack, S. Grover, N. Wolf, W. K. Metzger, D. Lu and G. Xiong, IEEE J. Photovolt., 2023, 13, 510–515.
65. C.-S. Jiang, H. R. Moutinho, R. G. Dhere and M. M. Al-Jassim, IEEE J. Photovolt., 2013, 3, 1383–1388.
66. M. Ishii, Y. Yoshino, K. Takarabe and O. Shimomura, Phys. B Condens. Matter, 1999, 273–274, 774–777.
67. J. Mäkinen, T. Laine, K. Saarinen, P. Hautojärvi, C. Corbel, V. M. Airaksinen and J. Nagle, Phys. Rev. B, 1995, 52, 4870–4883.
68. A. Nagaoka, D. Kuciauskas and M. A. Scarpulla, Appl. Phys. Lett, 2017, 111, 232103.
69. J. Jesenovec, C. Pansegrau, M. D. McCluskey, J. S. McCloy, T. D. Gustafson, L. E. Halliburton and J. B. Varley, Phys. Rev. Lett., 2022, 128, 077402.
70. M. M. Santillan, I. Chatratin, A. Janotti and M. D. McCluskey, Phys. Rev. Mater., 2024, 8, L111601.
71. C. Pansegrau and M. D. McCluskey, J. Appl. Phys, 2022, 131, 095701.

---