Effective chemical potential for non-equilibrium systems and its application to molecular beam epitaxy of Bi2Se3.

First-principles studies often rely on the assumption of equilibrium, which can be a poor approximation, e.g., for growth. Here, an effective chemical potential ([italic small mu, Greek, macron]) method for non-equilibrium systems is developed. A salient feature of the theory is that it maintains the equilibrium limits as the correct limit. In application to molecular beam epitaxy, rate equations are solved for the concentrations of small clusters, which serve as feedstock for growth. We find that [italic small mu, Greek, macron] is determined by the most probable, rather than by the lowest-energy, cluster. In the case of Bi2Se3, [italic small mu, Greek, macron] is found to be highly supersaturated, leading to a high nucleus concentration in agreement with experiment. This journal is © The Royal Society of Chemistry.

The concept of chemical potential is one of the most fundamental quantities associated with thermodynamic equilibrium. It generalizes the effect of the environment on the energetics of a single type of particle. Through the use of atomic chemical potentials, we can compare the energies of systems consisting of different numbers of atoms. In defect physics, this plays a central role in the calculation of the formation energies of defects and allows for a statistical determination of defect concentrations which determine the electronic properties of semiconductors. Further, surface reconstruction and even the shape of nanocrystals are intimately tied with the availability of different chemical species, with the surface energy generally depending on the atomic chemical potentials of the various species. While kinetic Monte Carlo (kMC) or molecular dynamics (MD) can be used to investigate system evolution, when comparing the energies of systems containing different numbers of particles, calculation within the grand canonical ensemble requires knowledge of the chemical potentials of the atomic species. The most ubiquitous approximation used to determine the range of chemical potentials is that of so-called near equilibrium growth (NEG),[1,3] in which the stoichiometric sum of the atomic chemical potentials is assumed to be equal to that of the bulk crystal. This assumption, however, is an outstanding issue in the theoretical community as growth itself is a highly non-equilibrium process.

The chemical potential, μ, is a thermodynamic quantity, defined by , where G is the Gibbs free energy and n is the number of the nth species.[4] For equilibrium or near-equilibrium growth (NEG), there is a well-defined relationship between the μ′s of the species in the system.[4] Consider, for example, a AB binary, we have

As a result of eqn (1), only one μ is an independent variable, i.e., one can use μB for this purpose and μA will be given by μA = (μAmBn − nμB)/m. Should μ for a species i exceeds its upper bound μ, defined as the chemical potential of the corresponding elemental solid or gas, this species will precipitate,[1,5] thereby preventing any further increase of μ. In other words,

Using the μ′s, the state of the system, whether it is a defect, a surface, or a crystal structure (denoted here as X) is determined by its formation energy with respect to bulk,where Etot is the total energy.[1]

This NEG theory has been applied to numerous materials behaviors with ample successes, ranging from surface reconstruction[5,6] and defect physics,[1,7] to nanostructures.[8,9] For systems in which a subsystem of concern can be isolated from the rest due, for example, to high reaction barriers, one can also use the NEG theory, provided that one can establish an approximate equilibrium within the subsystem.[10-12] However, when such a division is not obvious and the system is far away from equilibrium, the NEG theory can run into problems. For instance, recent interest in the physics of topological insulators (TIs) has spurred considerable activity to grow high-quality TIs such as Bi2Se3 with controlled properties.[13-18] Molecular beam epitaxy (MBE) is widely accepted as an effective method to grow high-quality films.[15-18] However, the defect density is still far from being satisfactory, due in part to the small size of the islands. Strictly speaking, MBE is not an equilibrium process;[19] its outcome may rely heavily on the dynamics of the source atoms/molecules. On the other hand, MBE is relatively simple for its slow growth rate and the simplicity of the sources. This raises the question: can one incorporate kinetics into the chemical potential-based approach for non-equilibrium process?

Note that considerable theoretical efforts have been made in the past to study growth. Analytical rate equation approach, based on empirical kinetic parameters, has been employed to study island formation and growth.[20-26] Statistical approaches such as kinetic Monte-Carlo simulation and molecular dynamics have been used to study film growth.[27-32] Molecular dynamics based on first-principles calculations could provide a vast amount of information, however, it is generally limited to smaller systems and shorter times due to the high computational cost. kMC offers a great advantage for efficiency and capability by simulating the time evolution of the system via statistically evaluating the probability of selected events. kMC relies on knowing the rate and mechanism of all the relevant transitions from a given initial state, with loss of reality. Moreover, it appears that these non-equilibrium approaches do not use the concept of chemical potential, which had served as a basis for NEG, except perhaps for the thermodynamic theory of nucleation.[2] Recently, a renewed interest in the chemical potential has emerged, as it can help choosing lowest-energy paths for growth,[33] without knowing relevant transitions in advance during a significant time scale.

In this paper, we derive a generic effective chemical potential () theory, which is suited for non-equilibrium physics. Taking the MBE growth of Bi2Se3 as an example, the use of allows us to simplify the overall complex process. We divide the non-equilibrium growth (nonEG) into three stages: pre-nucleation, nucleation, and island growth. For pre-nucleation, we determine by explicitly solving rate equations. First-principles calculation shows that is solely determined by the most probable, rather than by the lowest-energy, cluster(s) on the surface. For nucleation, it is found that the nucleation barrier vanishes when the critical size of the clusters is only a few molecule large due to the highly supersaturated nature of . This results in a high concentration of nuclei – a conclusion that contradicts the NEG model, but is in qualitative agreement with experiment.

One can define chemical potentials for an individual cluster α, irrespective of equilibrium,[4,34]where c is the areal concentration, G is the Gibbs free energy, E0 is the energy of the cluster on the surface relative to the isolated constituent atoms, NS is the number of sites per area that can hold the cluster (which, for simplicity, has been assumed to be much larger than the number of clusters), and k and T have the standard definitions. Here, considering a binary AB system (α = AB) as in NEG, we postulate that that enters eqn (3) is given bywhere the sums over s and t run over all possible AB clusters. As the concentrations of the clusters depend on the kinetics of the system, this weighted average incorporates both the energetics and kinetics in a simple way. When the system contains only a single element, eqn (5) is reduced towhere the subscript 0 denotes single-element quantity, as in eqn (2). For a binary system, as long as the formation of AB from A and B clusters is exothermic (i.e., with ΔE < 0), eqn (6) cannot hold except when species A (or B) starts to precipitate as pure clusters. Hence, it sets an upper bound for A. The same is true for species B. Therefore,

It is necessary to point out that eqn (5) and (7) are analogous to eqn (1) and (2), respectively.

To calculate , we solve the following rate equations,[2]where F is the deposition rate, kdes is the desorption rate, k− and k+ are the dissociation and association rates by emitting and absorbing a cluster i, respectively, D is the diffusion coefficient of cluster α and i[35] with a being the in-plane lattice constant, and ν ∼ 1013 s−1 is the vibrational frequency. The use of eqn (11) above assumes that the process is diffusion limited [see discussion in ESI and Fig. S1† for the general case]. For simplicity, in the following we consider only a flat surface, so the effect of the Ehrlich–Schwoebel barrier[2] can be ignored. We also assume a steady-state, i.e., letting ċ0 in eqn (8).

Note that in the limit when all the barriers in eqn (8)–(12) can be ignored, the system is dominated by its lowest-energy form, namely, the concentrations for all clusters, as well as for all finite-sized islands diminish, leaving only the AB bulk. In this case, eqn (5) and (7) become eqn (1) and (2) respectively, i.e., the system approaches its NEG limit. Also note that chemical potential is not an easily measurable quantity. People may, instead, use the off-stoichiometry of the system, which can be directly measured, as an indirect indicator of μ. This is true not only for the NEG model[36] but also for the nonEG model here.

First-principles calculations were performed using the VASP code[37] to determine the relevant cluster size, structure, and energy barrier. The projector augmented wave (PAW) method[38] and the local density approximation (LDA) to the exchange–correlation functional[39] were used. Plane waves with a kinetic energy cutoff of 210 eV were used as the basis set. Integration of the Brillouin zone was done with sufficient k-point sampling[40] such that the numerical error is less than 0.01 eV. All atoms were fully relaxed until forces on the atoms were less than 0.01 eV Å−1. The improved tangent finding method,[41] within the framework of the nudged elastic band (NEB) method, was used to determine the energy barriers. Experimental growth temperature of 500 K (ref. 15–17) was assumed in the rate calculations.

Denoting t0 as the time when nucleation takes place, the aforementioned 3 growth stages can be restated as follows: stage 1: t ≤ t0 (pre-nucleation), stage 2: t ∼ t0 (nucleation), stage 3: t ≫ t0 (island growth). Fig. 1 illustrates schematically the evolution of the size of the clusters, and those of the nuclei and islands, and the corresponding μ′s.

Fig. 1

A schematic illustration of cluster/island distribution and the corresponding chemical potentials in three different stages of Bi2Se3 growth: (a) pre-nucleation, (b) nucleation, and (c) island growth. t0 is the characteristic nucleation starting time.

In stage 1, there are only clusters, no nuclei, on the surface of Bi2Se3, whose thickness in the calculation is taken as five atomic layers, or called one quintuple layer (QL). Because the binding between QLs is van der Waals (vdW),[13-15] adding QLs to the thickness of the substrate is unlikely to affect our results. An in-plane supercell of 6 × 6 is used to conduct the cluster-adsorption calculation to avoid interactions among periodic images. Fig. 2 depicts the clusters we considered in our ab initio calculations. They are organized according to the chart in Fig. 2(a), where dashed line denotes the cutoff size of the clusters considered in this study (see details in Fig. S2, ESI†). Our calculations show that stable clusters on the surface are the same as those in vacuum. This can be expected as Bi tends to form 3 bonds, whereas Se tends to form 2 bonds, in line with the octet rule.[42]Fig. 3 shows the calculated diffusion barrier (φdiff), desorption barrier (φdes), and dissociation barrier (φ−) for the clusters: φdiff is generally small, often a couple tenth of an eV, except for atomic Se and BiSe2. φdes is higher, especially for smaller clusters. φ− is also generally higher and depends sensitively on the reaction pathway. Both φdes and φ− exhibit a large variation over a couple eV.

Fig. 2

(a) Small clusters for Bi2Se3 growth. They are organized in terms of size: horizontal = increasing number of Se atoms; vertical = increasing number of Bi atoms. Dashed zigzag line is discussed in the main text. (b)–(d) Calculated most-stable cluster structures on a Bi2Se3 (0001) substrate: (b) pure Se, (c) pure Bi, and (d) mixed BiSe clusters. Red atoms are Se whereas green atoms are Bi.

Fig. 3

Calculated barriers for (a) diffusion (φdiff), (b) desorption (φdes), and (c) dissociation (φ−) of the clusters in Fig. 2. Cluster names for both (a) and (b) are given in (b). For cluster dissociation in (c), a cluster labeled in the horizontal axis can have multi-reaction pathways, which are labeled inside or by the vertical bars.

Fig. 4(a and b) show c and the corresponding μ, which reveal that the highest concentrations of clusters are that of atomic Se and BiSe2, respectively. At first glance, it may appear counterintuitive that some clusters with lower μ, e.g., Bi2Se2 and Bi2Se3, also have low c. This is because of the low-desorption barrier of Bi2Se3 in Fig. 3, which depletes not only Bi2Se3 but also Bi2Se2 (via an association with atomic Se, see Fig. S1†). In MBE, the Se partial pressure is much higher than that of Bi. As a result, the population of Se is determined almost entirely by pure Se clusters viaeqn (7),

Fig. 4

(a) Logarithmic concentration of small clusters in stage 1. (b) The corresponding chemical potentials by eqn (4) and 's by eqn (13) and (14), respectively. (c) Logarithmic concentration of small clusters in stage 3. Solid squares are based on the nonEG model with cisland = 109 cm−2, whereas open squares are the prediction by the NEG model.

And for Bi, one has

These values, Bi = −4.44 eV and Se = −3.82 eV, are given in Fig. 4(b) as horizontal dashed lines. Our study reveals that it is the most-probable, not the lowest-energy, clusters that determine . This conclusion is in line with the basic principles of statistical physics, irrespective of the computational details, whereby it reinforces the notion that using c as the weighting factor for in eqn (5) is physically correct. It has been experimentally realized that the observed clusters were not the lowest energy ones in silicon cluster anions under non-equilibrium conditions.

Note that there is a large (3.08 eV) difference between the current model (2Bi + 3Se = −20.34 eV) and the NEG model (2μBi + 3μSe = −23.42 eV). The fact that the former is significantly higher than the latter gives rise to supersaturation during the growth. In other words, eqn (7) sets a new set of boundaries for chemical potentials as the prime reason that account for the differences between nonEG and NEG. Importantly, this finding allows one to reexamine defects in Bi2Se3, as it suggests that native defects may be more easily formed than predicted by NEG theory,[44] highlighting experimental challenges in growing high-quality epitaxial films.[45]

Before finishing up the discussion of stage 1, we would like to point out that our consideration of the availability of clusters is important, but only a first, step in improving upon NEG theory. Within the framework of the current development, higher order effects (such as the kinetics associated with the incorporation of available species into the growth front) can also be incorporated, which however will be deferred to future work. For stages 2 and 3, i.e., nucleation and island growth, due to continued deposition, we do not expect the concentrations of the small clusters to change considerably (as will be shown later).

In stage 2, owing to the formation of nuclei, bifurcation of μ will take place – one for the nuclei (μnuc) and one for the small clusters (), as depicted in Fig. 1(b). Bifurcation reflects the fact that once nuclei are formed, it is energetically favorable for them to grow further by consuming available smaller clusters on the surface, rather than consuming themselves. This should be contrasted to NEG where bifurcation of μ is strictly forbidden. Once for clusters are determined, one may approach the nucleation problem in different ways. For example, to study the temperature T dependence of nucleus density, one could apply eqn (4)–(12) to potential nuclei [which should be among the larger clusters in Fig. 2 and beyond (=embryo nuclei)] to determine μnuc and c, respectively. The non-negligible association barriers [see, e.g., Fig. S1(a)† for clusters] suggest that increasing T could help increasing the nucleus size and hence decreasing its density. However, the window for the increase can be limited, as although dissociation barriers are usually higher than association barriers [e.g., by comparing Fig. S1(a)† with Fig. 3(c)], the difference may not be that dramatic. As a result, further increasing T could also lead to the dissociation of already-formed nuclei, thereby increasing the nucleus density. Although rather qualitative, the conclusion here is already in line with experiments.[16,43]

If, however, our interest is only on the critical size of the nuclei, not on its T dependence, explicitly invoking kinetic theory here can be an overkill. Instead, we only need calculate the nucleation barrier, which is defined as the maximum of the Gibbs free-energy,[2]

After reaching the critical size, adding additional atoms/clusters to the nuclei is energetically favored. Using in stage 1, we obtain Fig. 5, showing that ΔG is remarkably small and even negative for (Bi2Se3) clusters with n > 3. Usually, ΔG for cluster of these sizes is significantly positive due to its high surface energy. Due to high and the relatively low energy of the vdW surface, however, this is no longer the case here. Fig. S3† shows that except for n = 1 and 3, the desorption barriers are reasonably large so the nuclei in Fig. 5 do grow into islands. Consequently, the island concentration is expected to be high. In other words, the small nuclei size is consistent with measured high island concentration, cisl ∼ 109 to 1010 cm−2.[16,17]

Fig. 5

Free-energy barrier ΔG for nucleation on various clusters. Shaded area is where only stoichiometric molecular clusters (Bi2Se3) are considered in the calculation.

In stage 3, the basic physics should be the same as in stage 2. Due to the presence of the islands, however, one should recalculate c and by adding a term −k+islcislc to eqn (8), where cisl is the concentration of the islands and k+isl = 2π(D). Fig. 4(c) shows the results for a typical cisl = 109 cm−2. More results on the cisl-dependence can be found in Fig. S4, ESI.† As it turns out, neither c nor is changed significantly. As a comparison, Fig. 4(c) also shows c under the NEG condition. They are many orders of magnitude smaller.

In summary, an effective chemical potential approach for non-equilibrium growth is developed. We find that for an atomic species is determined primarily by the most probable cluster in which the atom resides during growth. Because “most probable” is a balance between energetics and kinetics, our findings thus set a new criterion for most relevant events during growth. Application to MBE growth of Bi2Se3 suggests that is highly supersaturated, resulting in an exceedingly small critical size of nuclei. While a high concentration of islands is in agreement with experiment, our results also reveal that to grow better-quality films requires finding ways to stabilize the most probable clusters, thereby lowering . Our formulation is general. It may be used to study structure, surface morphology, and shape of a nanoparticle, as well as defect and impurity, in non-equilibrium-grown solids. While kinetic theory has been around for a long time, most first-principles studies today are still based on the NEG model. Our chemical potential-based development here offers a natural vehicle to transform such bulk studies into the more experimentally relevant non-equilibrium regime.

Conflicts of interest

There are no conflicts to declare.