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

Developing high fidelity quantitative models of solid state reaction systems can be challenging, especially in deposition systems where, in addition to the multiple competing processes occurring simultaneously, the solid interacts with its atmosphere. In this work, we develop a model for the growth of a thin solid film where species from the atmosphere adsorb, diffuse, and react with the film. The model is mesoscale and describes an entire film with thickness on the order of microns. Because it is stochastic, the model allows us to examine inhomogeneities and agglomerations that would be impossible to characterize with deterministic methods. We demonstrate the modeling approach with the example of chalcopyrite Cu(InGa)(SeS)$_2$ thin film growth via precursor reaction, which is a common industrial method for fabricating thin film photovoltaic modules. The model is used to understand how and why through-film variation in the composition of Cu(InGa)(SeS)$_2$ thin films arises and persists. We believe that the model will be valuable as an effective quantitative description of many other materials systems used in semiconductors, energy storage, and other fast-growing industries.


I. INTRODUCTION
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][4][5][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 largerscale 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 9 , 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 10 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 materi-als with complex microstructure. We use the chalcopyrite Cu(InGa)(SeS) 2 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 2 -based materials are commonly used as the absorber layer in thin film solar cells. Devices using these materials have demonstrated efficiencies exceeding 20% [11][12][13] .In order to increase voltage and improve efficiency, gallium and sulfur are alloyed with CuInSe 2 to form a continuous solid solution: Cu(InGa)(SeS) 2 . 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 2 Se and/or H 2 S 14 .
The reaction step can result in heterogeneous films with steep through-film composition gradients 15 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) 2 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.

II. MODEL DEVELOPMENT AND THEORY
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 9 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 i N 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.

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.

A. Solution Algorithm
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 ∈ N (N ×M ) 0 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) ∈ N 0 , 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 (l 1 adj , l 2 adj ), and (2) define a mapping from indices in L to index in X, or k = g((i 1 , j 1 ), (i 2 , j 2 )), where ((i 1 , j 1 ), (i 2 , j 2 )) are the indices of (l 1 adj , l 2 adj ). Each adjacent pair of lattice elements constitutes an interface, which can be represented by a single value {x|x ∈ N 0 , x < S 2 }, where there are S 2 "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 1 , j 1 ), (i 2 , j 2 )) 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: 4. Define an initial condition L 0 .
5. Use Equations 2 and 3 to map L 0 → 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": 7. Calculate the total propensity array, A ∈ R ((rows(D)+rows(N))×1) from each event in D and N, which has the elements: 8. Determine the time step, τ , drawn from the probability mass function: 9. Determine which lattice event occurs using probabilities: 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 18 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 7 time steps) run in several minutes on a typical personal computer.

B. Relation to Physical Properties
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. We now show how to derive physical property values from the value of the propensity constants.

Reaction Propensity
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−1 c 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 ℓ 2 , 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.

Adsorption Propensity
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

Desorption Propensity
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, 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 ℓ 2 , 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:

Diffusion Propensity
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 s 1 ,s 2 is the diffusivity and R(t ′ ) is the velocity autocorrelation function of species s 1 surrounded by species , denoting the inner product of v(t) and v(t + t ′ ). In our system, consider a single lattice point with value s 1 surrounded by an infinite lattice species s 2 . In this case, because there are four s 1 , s 2 interfaces, the propensity for diffusion is a s 1 ,s 2 = 4p s 1 ,s 2 , which means that the average time until the occurrence of a diffusion event is: Therefore, the velocity magnitude of species i is 4p s 1 ,s 2 ℓ. 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 1 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 nonisotropic effects or explicitly consider the presence of discrete crystalline grains. Therefore, our approach will strictly only be valid if (1)  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 2 Se 3 ), and (2) chalcopyrite formation (CuInSe 2 , CuInS 2 , 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 16 . This asymmetry in reaction rates leads to mostly CuInSe 2 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) 2 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 21 . 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.

B. Parameter Reduction
The reaction mechanism presented in Table I involves 12 reaction propensities, 4 adsorption/desorption propensities, and, in principle, 14 2 = 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) 2 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: 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 −14 m 2 s −1 . 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.

Slow Chalcopyrite Interdiffusion:
One feature commonly observed in Cu(InGa)(SeS) 2 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) 2 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. 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.

C. Composition Profile Prediction
The simulation algorithm from Section II was applied to the chalcopyrite growth model in Table I For a first case, we considered a process with only H 2 Se, for which the propensity for adsorption of H 2 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][16][17] .
For a second case, the precursors are reacted simultaneously in equal concentrations of 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 µ = pr 1−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, [26][27][28] . Now, if the formation of an agglomeration can be (f) MLE estimates (points) and fitted response surface for negative binomial random variable r.
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 Specifically, from the partial derivative of the response surface ( ∂r ∂x 1 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.

E. Model Predictions vs. Experimental Results
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 2 Se concentration of 1%, and varying H 2 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).
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 2 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 2 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 2 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 25 reported more substantial gallium homogenization than we observed, which may be explained by the lower H 2 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 2 Se before reaction is complete (the authors of 25 use a flow reactor, where H 2 Se is continuously replenished).

IV. SUMMARY AND CONCLUSIONS
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  Table I 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) 2 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 2 , for Cu 2 ZnSnS 4 , 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. 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 χ 2 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 χ 2 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 2 should approximate a χ 2 (ν) random variable with ν = m − 3 degrees of freedom. Thus, we calculate the probability that C 2 > χ 2 (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 χ 2 tests are shown in Table II and there is no evidence to reject the null hypothesis in any of the cases.

Appendix B: Response Surface Methods for Agglomeration Distribution
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 1 is film thickness and x 2 is Ga/(In+Ga). We employ a face-centered cubic response surface design and the results from each run are shown in Table II

Appendix C: Experimental Methods
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.8 Ga 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 −3 Pa to remove gas impurities, and then charged with H 2 Se, H 2 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.