A stochastic model of solid state thin film deposition: Application to chalcopyrite growth

Quantitative understanding of solid state reactions involved in film deposition and growth is important for improved processing in a number of industries including microelectronics, photovoltaics, energy storage, and pharmaceuticals. While there are many useful classical modeling approaches (from simple mass action kinetics^1,2 to detailed Avrami-type modeling^3–6) the increasingly complex material systems used in modern manufacturing require more sophisticated methods. Most modern techniques, however, employ microscopic level or ab initio approaches (e.g., Refs. 7 and 8), which promote fundamental knowledge but, because of current computational hardware limitations, are incapable of providing larger-scale property or composition predictions. In this paper, we present a mesoscopic thin film growth model capable of predicting film-scale composition.

In order to examine and characterize the lateral heterogeneity that can arise during film growth, the model we present is stochastic, rather than deterministic. Our approach is related to the stochastic simulation algorithm first developed by Gillespie,⁹ which assumes a uniformly mixed system with no mass transfer limitations. In solid state systems, however, mass transfer effects are always important and often rate-limiting. Gillespie’s method has recently been extended to include diffusion, mostly used for micro-scale modeling of biomolecular systems. Erban and Chapman¹⁰ discuss two general approaches: on-lattice methods and off-lattice methods. On-lattice approaches restrict the position of molecules to discrete locations or compartments where each compartment contains multiple occupant species, while off-lattice methods allow movement on a continuous domain, usually through Brownian motion. However, these approaches are designed for fluid systems where species density may vary. In this work, we apply similar concepts and expand the capability of stochastic models for reactions in the solid state. Our approach, discussed in Section II, is similar to on-lattice methods, but with the additional restriction that lattice occupancy is always exactly one. Instead of interacting with co-occupants, species interact with adjacent lattice points. We show how this approach allows stochastic simulation of mesoscale systems of solid, crystalline species, where unit-cell level (1 Å) simulation would be impractical for complete thin film (1–10 μm) systems. Although we assume the lattice is square, with a coordination number of four, our approach is easily generalized to allow for simulation of advanced, non-isotropic materials such as graphene, carbon nanotubes, and other materials with complex microstructure. We use the chalcopyrite Cu(InGa)(SeS)₂ system as an example to demonstrate how one can simulate film growth rate, composition profile, and agglomeration using the proposed stochastic approach.

Polycrystalline chalcopyrite CuInSe₂-based materials are commonly used as the absorber layer in thin film solar cells. Devices using these materials have demonstrated efficiencies exceeding 20%.^11–13 In order to increase voltage and improve efficiency, gallium and sulfur are alloyed with CuInSe₂ to form a continuous solid solution: Cu(InGa)(SeS)₂. The most common industrial process for producing these absorbers involves two steps in which a metal precursor (Cu-In-Ga) is deposited first and then reacted with gas-phase H₂Se and/or H₂S.¹⁴ The reaction step can result in heterogeneous films with steep through-film composition gradients¹⁵ and spatially-confined agglomerations.^16,17 The spatial heterogeneity resulting from the selenium and sulfur reactions makes this process an ideal system to demonstrate our solid state reaction model.

In the remaining sections, we present a novel stochastic model for solid state reaction kinetics, with emphasis on the ability to predict the composition profile and other spatial heterogeneities. In Section II we develop the model, describe an efficient solution algorithm, and explain the relationship between the model parameters and physical properties. Then, in Section III we describe a reaction mechanism for Cu(InGa)(SeS)₂ production, show how to apply the model using this mechanism, use the model to predict composition profiles and agglomeration statistics, and compare model predictions and experimental results. Finally, we offer conclusions and suggestions for future applications.

The system in question is a thin film in which solid state reactions occur, species interdiffuse, and the film interacts with its environment by adsorption and desorption of volatile species. We represent the film with a two-dimensional square lattice, where each point contains a species or a vacancy. The model is mesoscopic; so that each lattice point does not represent an individual atom, molecule, or unit cell, otherwise the lattice would be too large to be computationally tractable. The lattice is therefore a coarse-grained approximation of the actual film; each lattice point is a finite volume element small enough such that it is accurately approximated as phase-pure.

Our approach is to recast Gillespie’s stochastic simulation algorithm⁹ for spatially heterogeneous solid state systems with approximately constant mass density and number density. In Gillespie’s method, a random number is selected at each time step to determine which reaction occurs. Here, we generalize reaction events to “lattice” events, which take place at interfaces between lattice points and are classified as reaction, diffusion, adsorption, or desorption events. The probability and the rate of occurrence of each lattice event are governed by an intrinsic parameter called the propensity constant. The propensity of a given event is the product of its propensity constant and the number of interfaces at which the event can occur.

The modeling approach is as follows:

1. A square lattice is initialized with the starting species. If adsorption/desorption events are included, the lattice should contain vacancy points above the species. If the lattice is represented as an N × M array, row 0 and row N are considered boundaries with no interactions above row 0 points or below row N points. Column 0 is considered adjacent with column M (cf., periodic boundary conditions in a boundary value problem involving a partial differential equation).

2. Propensity of each lattice event is calculated as: a_i = p_iN_i, where p_i is the propensity constant and N_i the number of interfaces associated with lattice event i.

3. Probability of each lattice event is proportional to the propensity of an event; a random number is generated to determine which lattice event occurs.

4. The time, τ, until the next time step is selected from an exponentially distributed random variable.

5. The reaction chosen in Step 3 occurs at one possible interface. For example, if a reaction takes place between species A and B, then one of the A—B interfaces is selected at random and updated accordingly. For reaction events, the final orientation (that is, the relative position of the product species) is random; it is fixed for adsorption, and diffusion events.

6. The lattice is updated and steps 2–5 are repeated until an exit condition is met.

The conceptually simplest algorithm for implementing this modeling approach is to store the lattice in a 2D array that is updated at each time step. Although straightforward, this method is inefficient and will be computationally tractable only for relatively small lattice sizes. There are two possible bottlenecks that can lead to very slow solutions: (1) counting the number of interfaces of each kind, and (2) choosing one of these interfaces at random for a reaction event. To address these issues, we do not store and update the lattice itself at each time step; instead we track each interface and its position in an array. First, we define the simple 2D representation of the lattice; then we demonstrate how to convert this to the 1D representation that is used in the algorithm.

Consider a film discretized to a lattice and represented by the 2D array: $L∈N0(N×M)$ with elements l_(i,j). The indices of each element in L represent the position of each lattice point (i.e., volume element) in physical space, and the value of that element represents its occupant species. If there are S unique species, including vacant elements, and the value of each array element corresponds to its occupant species, then the domain of l_(i,j) is given as: {l_(i,j)|l_(i,j) ∈ ℕ₀, l_(i,j) < S}. Now, consider the set of adjacent points in L. For this model, we define the set of adjacent elements to be:

where the first two cases are trivially adjacent points and the third is analogous to applying periodic boundary conditions to a PDE to reduce edge effects. The periodic boundary condition is applied only to the columns of L, as the rows represent the full thickness of the thin film. Next, we convert the 2D representation, L, to a 1D representation X, with elements x_k by: (1) defining a mapping from species pairs to interface kinds, or, $x=f(ladj1,ladj2)$⁠, and (2) define a mapping from indices in L to index in X, or k = g((i₁, j₁), (i₂, j₂)), where ((i₁, j₁), (i₂, j₂)) are the indices of $(ladj1,ladj2)$⁠.

Each adjacent pair of lattice elements constitutes an interface, which can be represented by a single value {x|x ∈ ℕ₀, x < S²}, where there are S² “kinds” of interfaces (observe that interfaces of different orientation are considered distinct). We can then map pairs of species to interface kinds:

where x is the interface kind.

Next, we define a mapping for the indices for interface location in L to the index for interface location in X. Here, we assume that the number of columns, M, is an even number (a similar procedure applies if M were odd-valued). With ((i₁, j₁), (i₂, j₂)) as indices of adjacent elements in L and k as the index of X:

Using the mappings defined in Equations (2) and (3) to translate L → X, the simulation algorithm may now be written as follows:

1. Define a set of allowable species, {0, 1, 2, …(S − 1)}.

2. Define a set of directional lattice events, D. Directional events will preserve the orientation of the interface and are used to represent diffusion, adsorption, and desorption events. D is a (N_D × 3) matrix, where N_D is the number of directional lattice events. The columns of D correspond to {p, x₀, x_f} for the propensity constant, initial interface kind, and final interface kind, respectively (interface kind is the unique identifier for a pair of species defined in Equation (2)).

3. Define a set of non-directional lattice events, N. Non-directional events will not preserve the orientation of the interface and are used to represent reaction events. N is a (3 × N_N) matrix, where N_N is the number of non-directional lattice events. The columns of N correspond to {p, x₀, x_f} for the propensity constant, initial interface kind, and final interface kind, respectively. Directionality refers to the orientation of the products, not the reactants. For example, if there are species A, B, C, and D, then the reaction A + B ⟶ C + D is distinct from B + A ⟶ C + D.

4. Define an initial condition L₀.

5. Use Equations (2) and (3) to map L₀ → X

6. Count the number of interfaces of each kind in array X, saving the results in array Y, referred to as the “interface count array”:

   $ym=∑m∈X1,form∈N0,m<S2$

   (4)
7. Calculate the total propensity array, A ∈ ℝ^{((rows(D)+rows(N))×1)} from each event in D and N, which has the elements:

   $ai={p∈Di}yx0∈Diandai+rows(D)={p∈Ni}yx0∈Ni$

   (5)
8. Determine the time step, τ, drawn from the probability density function:

   $f(τ|A)=∑aiexp(−τ∑ai)$

   (6)
9. Determine which lattice event occurs using probabilities:

   $Pi=ai∑aj$

   (7)

10. From a uniformly random distribution, choose the reactive interface, i.e., the allowable interface at which the reaction takes place.

11. Update X and Y. The reactive interface (an element in X) and the interfaces adjacent to the reactive interface will be updated. The interface count array Y is updated by subtracting the previous the values of X and adding the new values of X to their corresponding elements in Y. (Refer to the code at the URL below for a complete definition of adjacent elements in X).

12. Check if exit condition is reached. If yes: continue; else: return to Step 7. The exit condition should be determined on a case-by-case basis. In this work, a specified simulation time is used and alternatives include a specified number of time steps or composition (though specifying a composition as an exit condition will not ensure that the composition will ever be reached).

13. Using the inverses of Equations (2) and (3), calculate L. End.

It should now be clear that the solution of the model is equivalent to sampling from a Markov chain where the elements of X define the system state; further, each step of the Markov chain will be assigned a step time in continuous time by assuming that each step is a Poisson process with intensity equal to the sum of all event propensities. The algorithm has been written in Python using Numpy¹⁸ and the code is available at www.bitbucket.org/rlovelett/stochastic_solid_state. Moderately sized simulations (about 1000 elements in L, 100 elements in D and N, and 10⁷ time steps) run in several minutes on a typical personal computer.

Determining the appropriate lattice size (i.e., the dimensions of L) requires a trade-off between the computational cost (reduced by using small lattices) and the variance of the model predictions (reduced by using large lattices). With the exception of edge effects that are present when simulating very small lattices, the dimensions do not affect the average model predictions (i.e., the model “accuracy”). However, changing the dimensions of L might inadvertently change the physical system that the model describes. For example, adding rows to L, without changing the other model parameters, would model a film with a greater thickness. Therefore, in the subsequent section, we show how physical properties are functions of both the propensity constants and the size of an individual lattice site. To ensure that model accuracy is not sacrificed when the number of rows in L is changed, the propensity constants must be rescaled using the functions that we derive in the next section (however, if the number of columns is changed, no other changes are necessary).

One of the advantages of a spatially distributed model is its ability to decouple the rates of different processes—reactions, adsorption, desorption, and diffusion. The physical properties that govern these rates can be related to the propensity constants used to advance the model and to the dimensions of the lattice. We now show how to derive physical property values from the value of the propensity constants and the size of lattice sites.

In its original implementation, the Gillespie algorithm’s propensity constant, c_i, is related to the macroscopic reaction rate constant, k_i according to k_i = V^b−1c_i, where V is the reactor volume and b is the reaction order. Our system is analogous, except that our reactions occur at surfaces. We assume that every formula unit on a lattice site is available for reaction at the surface, i.e., that diffusion within a lattice site is rapid. The characteristic surface area for reaction is ℓ², where ℓ is the length of a lattice site. Therefore, we can write the rate constant as:

where ρ_n,s is the area specific number density (units of inverse area) of formula units. With enough data at varying temperature, activation energies of the reactions can be estimated using the Arrhenius equation.

The rate of adsorption can be described by the product of the rate of collisions between a surface and a gas, and a sticking coefficient, S_i, or the probability that the species will stay on the surface. In some cases, like the chalcopyrite growth system we examine in Section III, the probability of dissociation of a gas phase species should also be included. For simplicity, however, probability of hydride gas dissociation will be lumped with the sticking coefficient. First, from kinetic theory, we can determine the rate of collisions with the surface per area, $Fi,a=Pi/2πMikbT$⁠, where F_i,a is the adsorptive flux of species i (m⁻²s⁻¹); P_i, the partial pressure; M_i, the molecular weight; k_b, Boltzmann’s constant; and T, the absolute temperature. Therefore, we can determine the rate of adsorption on a single lattice site using the area of that lattice site: r_ads = F_i,aSℓ². However, this rate is in molecules per time, not lattice sites per time; thus, it should be scaled by the number of molecules in a lattice site, or ρ_nℓ³, where ρ_n is the number density of the species. Since 1/p_i is the average time until an adsorption event occurs on a single site, p_i is the stochastic equivalent of the rate of adsorption on a single lattice site. If the propensity constant is known, then we can derive an expression for S_i as:

The adsorption propensity is an important adjustable parameter that can be used to model several common situations. First, adsorption may be “selective”, that is, an adsorbent is more likely to adsorb on one particular substrate than on others. Selective adsorption can be modeled by adjusting the adsorption propensity for each pair of possible adsorbent and substrate species. Second, the concentration and/or flux of adsorbing species may be time-dependent for some manufacturing processes. The adsorption propensity can be made a function of time in order to model these dynamic deposition processes.

Desorption is physically equivalent to evaporation. Unlike other lattice events, however, the rate of evaporation is not a function of thermodynamic properties only; it depends also on system-specific parameters such as gas flow rate and reactor geometry. Therefore, the rate of evaporation is best captured by a mass transfer coefficient, k_m,i = F_i,e/Δf_i, where F_i,e is the evaporative flux, and Δf_i is the difference in the species i fugacity between the adsorbed phase and gas phase. Assuming that the gas phase is ideal, the fugacity difference reduces to (P_vap,i − P_i). Similar to the adsorption case, we can determine the rate of evaporation from a single lattice site, r_evap = F_i,eℓ², scale it by the number of molecules in a lattice site, and set the scaled evaporation rate equal to the propensity constant. Solving for k_m,i yields:

Multicomponent systems with significant diffusion limitations can be challenging to model appropriately. However, because our modeling approach involves binary interactions between species, we can relate the propensity constants for diffusion events to binary diffusivities. We will use a well-known result from statistical physics to obtain the diffusivities, where diffusivity can be written as:

Here D_s1,s2 is the diffusivity and R(t′) is the velocity autocorrelation function of species s₁ surrounded by species s₂, defined as R(t′) = 〈v(t) ⋅ v(t + t′)〉, denoting the inner product of v(t) and v(t + t′). In our system, consider a single lattice point with value s₁ surrounded by an infinite lattice species s₂. In this case, because there are four s₁, s₂ interfaces, the propensity for diffusion is a_s1,s2 = 4p_s1,s2, which means that the average time until the occurrence of a diffusion event is:

Therefore, the velocity magnitude of species i is 4p_s1,s2ℓ. Choosing the current time to be immediately before a diffusion event, and recognizing that our model consists of discrete time steps, we can rewrite Equation (11) as:

where v_k is the velocity of species s₁ during time step k, Δt_k is the duration of time step k, and initial element of the series (k = 0) is moved outside the summation. Considering that the direction of the velocity vector v_k will be chosen randomly from two orthogonal unit vectors and their inverses at each time step, we conclude that the sum will converge to zero, and that only the average value of the first term should be retained. Therefore, by using the average velocity magnitude and average time until the occurrence of an event given in Equation (12), the binary diffusivity is obtained as:

Finally, we should note that our approach is unconventional for describing diffusion in solids. Typically, the crystallinity of the material should be taken into account explicitly and non-isotropic effects should be considered. Our analysis, however, does not include non-isotropic effects or explicitly consider the presence of discrete crystalline grains. Therefore, our approach will strictly only be valid if (1) the crystal does not show preferential orientation and (2) grain interior diffusion (such as interstitial or vacancy-mediated diffusion) occurs at a rate similar to or greater than grain boundary diffusion. The second condition could be very restrictive, as grain boundary diffusion usually dominates in polycrystalline films. However, if this assumption were to be violated, the result would be that the diffusion propensities in the model would correspond to effective diffusion coefficients (via Equation (14)). In this case, the diffusion coefficient would not correspond to the energy of a fundamental reaction step or obey the Arrhenius relation that is typical of solid state diffusivities. However, in cases where defect-mediated diffusion is the dominant mechanism of mass transfer, information about the defect concentration can be incorporated into the model for more accurate estimation of the diffusion propensity.

To apply the model to the reaction of Cu-In-Ga precursors with H₂Se and H₂S, we require a reaction mechanism. Because the model is driven by interactions between pairs of species, we restrict the mechanism to only 2-species reactions. Note that for most solid state reaction systems, this constraint is not very limiting. Several groups^6,15,19,20 have suggested plausible reaction pathways. The specific phases or species involved vary somewhat among groups, suggesting that the reaction pathway may be process dependent. However, some elements are common to each mechanism. First, there are at least two stages to the reaction: (1) metal chalcogenide formation (e.g., InSe, InS, Ga₂Se₃), and (2) chalcopyrite formation (CuInSe₂, CuInS₂, etc.). Second, although reaction rates are often not determined quantitatively, it is suggested that the reaction of Se and In is faster than the reaction of Se and Ga.¹⁶ This asymmetry in reaction rates leads to mostly CuInSe₂ near the front of the film and Ga-containing species accumulated at the back. Next, we note that the reaction mechanism we propose implies that the Cu(InGa)(SeS)₂ films are stoichiometric and that the discrete nature of the model enforces the stoichiometry. In practice, however, the precursors (and consequently the final films) are deposited as copper-deficient with the ratio Cu/(In+Ga) ≈ 0.9.^14,21 We do not treat copper deficiency explicitly, though the diffusion coefficients (and, in this model, diffusion propensities) would change if stoichiometric or copper rich-film films were modeled because indium and gallium may diffuse through copper vacancies.^22,23 Finally, for the initial condition, we assume that a mixture of CuIn and CuGa binary species is an adequate simplified representation of an actual metal precursor, which typically contains more complex Cu-Ga-In alloys and elemental In.²¹ Consequently, we propose that the system can modeled using the mechanism in Table I and the values given for propensity constants are discussed in the following section.

The reaction mechanism presented in Table I involves 12 reaction propensities, 4 adsorption/desorption propensities, and, in principle, $142=91$ unique diffusion propensities (though many are neglected). Furthermore, the lattice size may affect the model results (see Section II B) and will greatly affect the computation time. With such a large parameter set and a computationally intensive model, conventional parameter fitting is impractical. We present three simplifying assumptions and heuristics and show how we can use them to guide us in determining physically meaningful estimates for the model parameters.

1. Parabolic Film Growth: Results from the literature^5,6 suggest that Cu(InGa)(SeS)₂ films produced via reaction of metal precursors follow a parabolic growth mechanism, referring to a solid state reaction process where there is a planar, advancing reaction front, rather than a nucleation and growth mechanism. Invoking this mechanism suggests that the gas phase reactants, Se and S, can diffuse through reacted species, but not through the original CuIn and CuGa species. Therefore, diffusion propensities are set such that no species can diffuse with precursors—except for precursors themselves (CuIn and CuGa), which may diffuse with each other.

   Recognizing that the diffusion of Se and S are rate limiting in the parabolic mechanism, we can estimate the magnitude of the diffusion coefficients using the characteristic diffusion time:

   $τD=L2Ds1,s2$

   (15)

   Since the reaction takes place on the order of minutes and the film thickness is approximately 2 μm, a reasonable estimate for the diffusivity of Se is 6.7 × 10⁻¹⁴m²s⁻¹. After selecting a lattice size, the diffusion coefficients can be used to estimate the propensity constant for diffusion from Equation (14). In the simulations presented below, the lattice size is 100 nm, suggesting 1.7 (which we truncate to 1.0) is a reasonable estimate of selenium diffusion propensity.

2. Slow Chalcopyrite Interdiffusion: One feature commonly observed in Cu(InGa)(SeS)₂ films produced by reaction of metal precursors is a persistent gradient in gallium. That is, a gallium gradient forms and does not quickly anneal out to yield a uniform film.^{15,17,19,23,24} However, because Cu(InGa)(SeS)₂ is a continuous solid solution, the gradient is not a thermodynamic limitation, but must be limited by mass transfer. Therefore, we set the diffusion propensities of fully reacted species with each other to be zero, as their diffusion time is longer than the time scale of a typical reaction.

3. Fast Precursor Interdiffusion: In contrast to fully reacted chalcopyrite species, the unreacted species must interdiffuse relatively quickly. The fast interdiffusion of CuIn and CuGa (at least faster than the time scale of the diffusion and reaction of Se) is required for the reaction asymmetry to cause composition gradients.

The three heuristics presented above were used to guide parameter selection, especially the diffusion propensities for all sulfur-free lattice events (sulfur-containing lattice events are discussed in the next section when sulfur-containing models are presented). From the first heuristic, to ensure a reasonable time scale, the diffusion propensities of selenium with chalcogenides and chalcopyrites were set to 1.0. Based on the second heuristic, most of the remaining diffusion propensities were set to zero, except for interdiffusion of precursor species (CuIn and CuGa). The third heuristic compels us to select a diffusion propensity for CuIn and CuGa that is larger than that of selenium and reacted phases; therefore we selected 20.0 for the baseline. As discussed earlier, the rate of reaction of Se and In species is much faster than Se and Ga, therefore a baseline estimate for the propensity constant is 50.0 for InSe formation and 1.0 for GaSe formation. The remaining reaction propensities (chalcopyrite formation reactions) did not have a substantial effect on the composition profile and were set to 0.1. One possible approach for estimating adsorption and desorption propensities would be to use first principles and the equations derived in Section II B. However, unless reasonable estimates for sticking coefficient are available (in many cases, its order of magnitude is much less than 1), we recommend choosing a value of similar order to the propensity constants of the other processes to reduce computational burden. For our case, because we do not expect significant accumulation of selenium in its elemental form (due to its high vapor pressure and the stability of H₂Se molecules), the propensity of adsorption should be less than propensity for desorption. Therefore, for selenium, we set adsorption propensity to 0.1 and desorption propensity to 5.0. These propensities were used for the baseline simulations presented in the next section and are shown in Table I.

The simulation algorithm from Section II was applied to the chalcopyrite growth model in Table I. The initial lattice, L₀, is 23 rows × 100 columns, with rows 0 to 11 specified as vacancy elements, and rows 12 to 22 specified as a 0.25:0.75 mixture of CuGa and CuIn elements (this ratio was chosen for industrial relevance). Because Cu(InGa)(SeS)₂ films for photovoltaic cells are typically about 2 μm thick, the 23 rows result in lattice length ℓ ≈ 100 nm. The algorithm was applied to advance the lattice until 40 minutes simulation time elapsed. A useful way to visualize the results from the simulations is to plot the Ga/(In+Ga) and S/(S+Se) ratios as functions of depth.

For a first case, we considered a process with only H₂Se, for which the propensity for adsorption of H₂S was therefore set to zero. The results are shown in Fig. 1, where a large amount of gallium accumulates near the back of the film, which is a feature that has been well-documented experimentally.^15–17 

For a second case, the precursors are reacted simultaneously in equal concentrations of H₂Se and H₂S. The simultaneous reaction process is modeled by setting adsorption of sulfur lower than adsorption of selenium, but diffusion of sulfur faster than diffusion of selenium (the effect of sulfur adsorption is discussed in more detail in Section III E). Through-film profiles of gallium and sulfur are shown in Fig. 2. Consistent with experimental results from Ref. 25, gallium is distributed more homogeneously than in a process with H₂Se only.

The gallium and sulfur profiles in Figs. 1 and 2 result from complex interactions of the different propensity constants. One way to understand how these profiles arise is to examine the dynamics of the process. A simplified representation of the dynamics is shown in Fig. 3, where composition is spatially averaged through the entire film and shown as a function of time. In the first case, because the propensity for reaction with CuIn is much faster than with CuGa, most CuGa remains in its precursor form until nearly all of the CuIn is converted to chalcopyrite. During this time, CuGa diffuses toward the back contact leading to the profile in Fig. 1. The observation that Ga-containing species diffuse toward the back contact, while In-containing species are incorporated in the chalcopyrite phase (CuInSe₂), is supported by experiments from the literature. For example, in Ref. 16, x-ray diffraction patterns indicate that Cu₉(In_.2Ga_.8)₄ is present on the back contact after a 10 minute reaction in H₂Se, but only Cu₉Ga₄ is present after 30 and 90 minute reactions. For the second case, with H₂Se and H₂S, there is first a rapid increase in CuInSe₂ because of the relatively faster adsorption of Se compared to S; however, the faster diffusion of sulfur, combined with its preference to react with CuGa, leads to more incorporation of CuGa earlier in the process than if there were H₂Se alone.

One of the advantages of applying stochastic simulation at the mesoscopic scale is that it allows one to study fluctuations in local composition explicitly. The fluctuations result in some spatially confined features that can be described statistically in the example system, Cu(InGa)(SeS)₂. One of the important features observed experimentally is agglomerations of Ga-rich species at the back contact.¹⁶ In principle, there are two causes of these agglomerations: thermodynamic phase separation, and random fluctuations in species positions. The modeling approach cannot account for thermodynamic phase separations, but we can use it to examine the effect of random fluctuations. Although thermodynamic driving forces are likely important in Cu(InGa)(SeS)₂ films (evidenced by formation of Cu₉Ga₄ phase near the back contact^16,17), we still present the model results to demonstrate our approach for characterizing agglomeration size distribution.

For this purpose, simulations with large lattice sizes (2500 columns) were run with no adsorption or reaction events allowed and with precursors (i.e., L₀ matrices) of varying composition and thickness. We define any isolated collection of CuGa elements to be an agglomeration, regardless of size. Thus, a single, isolated CuGa element is considered a size-1 agglomeration. Figs. 4(a)–4(d) show histograms of the size distribution of agglomerations for different film thicknesses and compositions (which models later stages of the reaction, where the precursor film gets thinner as the front of the film is converted to chalcopyrite).

We consider the negative binomial probability distribution function as an appropriate model for quantifying the effects of composition fluctuations on agglomeration formation for the following reasons: The negative binomial random variable is, by definition, the number of Bernoulli trials with probability p until r “successes” are reached. Further, r can be thought of as a “waiting time” parameter and need not be integer-valued. More generally, the negative binomial random variable is an overdispersed form of a Poisson random variable with mean $μ=pr1−p+1$⁠, but with gamma-distributed intensity. In this form, the r parameter is often referred to as an “aggregation” or inverse dispersion parameter, where as r → 0, σ² → ∞, where σ² is variance.^26–28 Now, if the formation of an agglomeration can be considered as a series of trials, where p is the probability of adding another species to the agglomeration, then the agglomeration size should follow a negative binomial distribution with r = 1 (i.e., a geometric distribution). However, because of its improved handling of highly dispersed data, we use the more general negative binomial form allowing both p and r to be estimated. Intuitively, an increasing p or increasing r, will move the distribution rightward, toward a tendency for larger agglomerations.

The expression for the negative binomial probability distribution,

is used to characterize agglomeration size distribution by estimating the two parameters, p and r, using Maximum Likelihood Estimation (MLE) (see  Appendix A) from the simulation data.

In order to uncover the effects of geometry (i.e., film thickness) and composition (Ga/(In+Ga)) on the parameters, (p, r), we constructed a face-centered cubic response surface design, and used it to develop an empirical model. In the response surface model, geometry and composition are considered factors that affect the response variables, p and r. Nine simulations were run with varying lattice thickness and Ga/(In+Ga) ratios; p and r were estimated with MLE; and the estimates were used to fit 2^nd order polynomial models of p and r as a function of lattice thickness and composition (see  Appendix B). The fitted response surfaces, i.e., the polynomial models, are shown in Figs. 4(e)–4(f), from which we conclude that the probability of success parameter, p, is a function of only composition. Thus, all else being equal, when the fraction of CuGa sites increases, the agglomerations of CuGa tend to grow larger. The waiting time/aggregation parameter, r, has a more complicated response surface that is strongly affected by the interaction of geometry and composition. Specifically, from the partial derivative of the response surface (⁠$∂r∂x1$ of Equation (B2)), films with a gallium fraction less than 0.33 form larger agglomerations (due to larger r) when thickness is greater; the inverse is true for films with a gallium fraction greater than 0.33.

In order to compare our model predictions to experimental data, we produced a series of samples in the laboratory with the following process conditions: reaction temperature of 550 °C, H₂Se concentration of 1%, and varying H₂S concentration in Ar gas. (See  Appendix C for a summary of experimental methods). Sample compositions were measured using Energy Dispersive X-ray (EDX) spectroscopy with results for Ga/(In+Ga) shown in Fig. 5. While sulfur was detected in each sample, the composition was too low (less than 1 atom percent) for accurate quantitative measurement. The EDX spectroscopy measurements have a depth sensitivity of approximately 0.7 to 0.9 μm (or about one half of total film thickness). Therefore, EDX measurements can be treated as spatial averages of composition in the front half of the film. Each precursor sample has an average Ga ratio of 0.25; therefore, if the measured Ga ratio is less than 0.25, then the Ga is segregated toward the back of the film.

Without time- and depth-resolved data, conventional parameter fitting is impossible. Further, even if the data were available, parameter fitting would be very computationally intensive because it requires numerous simulation runs. Instead, we compare the effect of varying H₂S concentration in the gas phase to the equivalent change in our model, that is, changing the propensity for adsorption of sulfur. We expect that varying H₂S concentration will affect the through-film gallium profile, and that this effect will be captured in experimental results and by our model. However, rather than absolute convergence (which would require more precise parameter estimates), similar trends should be observed in experiments and simulations.

Fig. 5 shows the through-film gallium profiles from the model using the baseline simulation parameters, except for adsorption of sulfur, which varies between simulations. Also displayed in the figure is the measured gallium ratio of the films produced with varying H₂S concentration in the gas phase. In both simulation and experiment, when the sulfur (propensity in simulations or concentration in experiments) is low, gallium fraction increases near the front surface with increasing sulfur, but the effect is diminished as sulfur increases further. However, our simulation overestimates the effect of sulfur on gallium fraction. Several mechanisms may explain the discrepancy; for example, the reaction between selenium and indium or the adsorption rate of selenium may be faster than our estimate, which limits the tendency of sulfur to increase gallium homogenization. Other authors²⁵ reported more substantial gallium homogenization than we observed, which may be explained by the lower H₂Se concentration (0.35%) in their experiments. However, we cannot use such low concentrations in our reactor, which operates in batch mode, because the gas phase would be depleted of H₂Se before reaction is complete (the authors of Ref. 25 use a flow reactor, where H₂Se is continuously replenished).

We have presented a novel method for modeling thin film growth using a stochastic simulation method. We demonstrated an algorithm that makes the model computationally tractable for a typical desktop PC. In particular, we showed how the model applies to the industrially relevant case of thin film Cu(InGa)(SeS)₂ growth using a precursor reaction process. The model explains how the complicated, experimentally observed through-film profiles in Ga and S arise from complex interactions of reaction rates, adsorption rates, and diffusion limitations. We show that the stochastic nature of the model allows it to be used to understand lateral inhomogeneities, such as agglomerations, that would otherwise be ignored or averaged out in continuum approaches.

We believe that this modeling approach can find wide application in a number of thin film or other solid state material systems. In particular, we suggest that this approach will be especially useful for vacuum-deposited Cu(InGa)Se₂, for Cu₂ZnSnS₄, and for silicon systems with impurities because of similarities to our model system, including the potential for composition profiles and lateral heterogeneity. Furthermore, because of our method’s emphasis on material adjacency, it can be applied to systems with complex geometry, such as graphene or carbon nanotubes. More fundamentally, however, the stochastic simulation that we developed is a new approach that could be generalized for any system where system evolution is governed by network connectivity (in our case, material adjacency) instead of bulk composition, and may be applied in entirely unrelated fields.

This material is based upon work primarily supported by the Quantum Energy and Sustainable Solar Technology Engineering Research Center (QESST ERC), the National Science Foundation (NSF) and the Department of Energy (DOE) under NSF CA No. EEC-1041895. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect those of NSF or DOE.

In this appendix, we use a slightly altered version of the negative binomial probability model (first shown in Equation (16)), given by

to characterize the agglomeration size distribution. In this form, the domain of x has been shifted leftward by 1 unit; thus the data should be shifted down by 1 unit (i.e., size-1 agglomerations will be considered size-0 agglomerations), but the parameter values remain unchanged. The two parameters, r and p were estimated using Maximum Likelihood Estimation (MLE). The likelihood function for a sample of agglomeration data, X is then obtained as:

and for a total of N agglomeration samples, the log-likelihood function is:

To maximize the log-likelihood function (equivalent to maximizing the likelihood function) the partial derivatives of the log-likelihood function with respect to the paramters are set to zero:

where Ψ(x) is the digamma function, or Γ′(x)/Γ(x). While the second of these equations can be solved explicitly for p, in general, there is no closed form solution to Equations (A4) for r and p. The system of equations is solved numerically to obtain the parameter estimates

In order to validate this approach, simulations were run with varying film thicknesses and gallium fractions, p and r were estimated, and a χ² test was applied to compare model predictions to data. The response surface method (see  Appendix B) was used to select the specific values of film thickness and gallium fraction. The χ² test statistic is:

f_i is the observed count of agglomerations of size i, ϕ_i is the predicted count using the negative binomial distribution with MLE parameter estimates, and m is the largest agglomeration size with at least 5 instances. If the model is appropriate, C² should approximate a χ²(ν) random variable with ν = m − 3 degrees of freedom. Thus, we calculate the probability that C² > χ²(m − 3) and if this value is less than 0.05 (a commonly-used significance level), then we reject the null hypothesis that that model is appropriate. The results of the χ² tests are shown in Table II and there is no evidence to reject the null hypothesis in any of the cases.

Response surface methodology is an experimental design approach usually used for optimizing a process with a response variable that is approximated as a 2^nd order function of several input (or factor) variables. In this work, we use the methodology to understand empirically the effects of film thickness and composition on the negative binomial random variable parameters, (p, r). We postulate that the following model is appropriate:

where y is the response variable (p or r), x₁ is film thickness and x₂ is Ga/(In+Ga). We employ a face-centered cubic response surface design and the results from each run are shown in Table II (traditionally, two or more replicates of the center point are included—we include only one replicate because of the low variance resulting from simulated, rather than experimental, data). The parameters (β_i) were estimated using ordinary least squares and insignificant effects (p > 0.05, assuming normally distributed error) were removed from the models. The resulting response surfaces (plotted in Figs. 4(e)-4(f)) are:

The films analyzed in Section III E were produced as follows. Soda lime glass substrates (2.5 cm × 2.5 cm) were coated with 700 nm of Mo followed by 700 nm of a Cu-In-Ga mixture using DC magnetron sputtering. The Cu-In-Ga layers were rotationally sputtered from Cu_0.8Ga_0.2 and In sputter targets to yield final Ga/(In+Ga) ratio of 0.25. The samples were then placed on a graphite sample holder and inserted in a 5 cm diameter quartz reactor tube. The reactor was evacuated to at least 6 × 10⁻³ Pa to remove gas impurities, and then charged with H₂Se, H₂S, and Ar at the desired concentrations. Samples were heated to 550 °C using a 1000 W quartz-halogen lamp and maintained at that temperature for 10 minutes. Each sample composition was measured using an Oxford Instruments PentaFET 6900 EDX detector in an Amray 1810 scanning electron microscope with the samples in plan view. The microscope acceleration potential was set to 20 kV which corresponds to an EDX depth sensitivity of approximately 0.8 μm, or one half of the film thickness after reaction.