Atomistic mechanisms of orientation and temperature dependence in gold-catalyzed silicon growth

Semiconductor nanowires (NWs) are promising candidates for building innovative nano-devices,¹ due to their unique electrical and optical properties. For example, silicon (Si) NWs have a diameter-dependent band gap that is tunable by impurity doping,² and low-dimensional Si nano-structures have a higher optical absorption than bulk Si.³ The vapor-liquid-solid (VLS) growth is a commonly used technique for semiconductor NW synthesis.^4,5 The necessity to grow long, reproducible and defect-free NWs has motivated extensive research on the VLS NW growth mechanism, both theoretically^6–10 and experimentally.^11–14 However, dynamic observations of VLS nucleation and growth at the atomistic level have not been fully achieved, due to the limitations of current experimental characterization techniques. Thus, atomistic simulations, such as molecular dynamics (MD), can serve as a useful tool to integrate the findings from experiments, and to construct bridges to continuum models.^15–19 In the literature, the capability of MD simulations have been demonstrated for studying the gold (Au)-catalyzed Si growth via the VLS mechanism.^8,9 But limited by the computing power of the day, the simulation time in these studies was relatively short (several ns).^8,9 Thereby, the simulation conditions were usually in a deeply saturated and high temperature regime.⁸ 

To systematically study the VLS growth of the Si-Au system, we perform MD simulations under various orientation, Si fraction, and temperature conditions, with much longer simulation time (100 ns) than that reported in the earlier studies. From these simulations, the Si growth velocities are obtained, based on which the relationship among the growth velocity, temperature, and Si supersaturation is discussed for both $〈110〉$ and $〈111〉$ orientations. We explain the orientation dependence of the Si crystal growth by distinguishing different growth mechanisms between the two substrate orientations. Particularly, we identify that the growth on the {111} substrate is nucleation controlled, but nucleation is faster than island expansion. This is likely caused by the Au atoms adhering to the liquid-solid interface when a new island is nucleated. For the temperature dependence of the growth velocity, we find that it is more controlled by the interface atom activity rather than diffusion in the liquid.

The rest of the paper is organized as follows. A detailed simulation setup will be provided in Sec. II. Key simulation results and discussions, regarding the growth, the nucleation process, and their orientation and temperature dependences, will be presented in Sec. III. A conclusion will be given in Sec. IV.

An MEAM potential for the Si-Au system²⁰ is adopted for this study. Figure 1(a) shows the Si-Au binary phase diagram predicted by the potential together with the experimental phase diagram.²¹ On this diagram, the temperature and Si fraction conditions that have been studied previously are marked with “+”⁸ and “*.”⁹ The MD simulation conditions in this work, which are relatively closer to the eutectic point, are marked as “▲,” “▪,” and “▼.”

The “▲” markers represent the simulation conditions chosen for investigating the Si crystal growth speed, specifically, at temperatures of 900, 950, and 1000 K, and Si fractions of 40%, 45%, and 50%.²³ For each point on this phase diagram, six MD simulations are performed, three for Si growing on the {111} substrate and the other three for the growth on the {110} substrate. For observing the nucleation process and identifying the growth patterns, additional simulations are performed at 850 K with the initial Si fraction in liquid ranging from 27% to 29% (the liquidus composition is 26.5% at 850 K), as labeled by the “▪” symbol and also shown in the inset. To reveal the atoms' behavior at solid-liquid equilibrium, further simulations are performed at conditions along the liquidus line, marked by the “▼” symbol in Fig. 1(a). Table I summarizes all the above simulation conditions and the sections in which they are discussed.

For simulations on the {111} substrate, as an example, Fig. 1(b) provides a 3D view of the initial configuration at 900 K and 40% Si fraction. The random mixture of yellow colored Au atoms and grey colored Si atoms forms the liquid alloy on top of the crystalline Si substrate. The simulation cell consists of 26 144 atoms in total, with the orientations of the box axes chosen to be $[1 2¯ 1]$ in x, [1 1 1] in y, and $[1¯ 0 1]$ in z, respectively. At 900 K, the dimensions of the simulation box are 85.47 Å × 110.66 Å × 81.86 Å. As the periodic boundary conditions are applied to x and z, thermal expansion of the crystalline Si substrate is considered, such that the simulation box size at a given temperature is adjusted based on the calculated expansion coefficient of Si. It should be noted that the total number of atoms is kept the same for all Si compositions. Therefore, the ratio of the Au atoms to the Si atoms is changed in the liquid to reach the specified Si composition.

For {110} growth simulations, as illustrated in Fig. 1(c), the initial configuration setup is similar to that for {111} growth simulations, except that in this case, the orientations of the box axes are [0 1 0] in x, [1 0 1] in y, and $[1¯ 0 1]$ in z, containing 25 184 atoms in total.

The simulations are performed under the canonical (NVT) ensemble, using the Nosé-Hoover thermostat to control the temperature. At each given temperature and initial liquid Si composition, MD simulations are run for 100 × 10⁶ steps using a time step to 1 fs with LAMMPS.²⁴ The bottom two layers of the substrate are fixed to prevent the center of mass of the simulation cell from drifting. The wall-clock time for each MD simulation is around 10 days using 48 CPUs. The simulation results presented in this study are obtained from 70 independent MD runs with a total of ∼800 000 CPU hours.

In order to distinguish Si atoms in the crystal from those in the liquid, we use an order parameter q₃.²⁵ Atoms with q₃ > –0.5 are identified as in the crystalline state. (The selection of the threshold q₃ value is discussed in Sec. S1 in the supplementary material.) The crystal growth speed can thus be quantified by plotting the number of crystalline Si atoms N_Si as a function of time t. For instance, Figs. 2(a) and 2(b) present the N_Si -t relationships for Si growth along $〈110〉$ and $〈111〉$ orientations, at 950 K with initial Si fractions of 40%, 45%, and 50%, respectively. The curves with the same color correspond to the three simulations at the same temperature and Si fraction. In Fig. 2, all the curves initially rise up quickly, and then enter a quasi-steady-state growth regime at around 30 ns. The ramping-up of the curve at the beginning stage may be an artifact of the initial condition, in which the substrate is perfectly flat as shown in Fig. 1, so that some time is required for the liquid-solid interface to reach an equilibrium morphology. Therefore, only the data after 30 ns are used for the growth velocity analysis.

The atom attaching rate $dNSidt$ can be obtained by extracting the slope of the curves via a linear fitting, which is performed at 60, 80, and 100 ns, using the data over a time span of around 10 ns, marked as three rectangular regions in Fig. 2. During the simulation, as the crystal grows, the liquid gradually becomes depleted of Si. Therefore, the actual Si fractions in liquid at 60, 80, and 100 ns are separately measured from the MD configurations. (Detailed discussion on how to reliably extract Si fraction in the liquid is provided in Sec. S2 in the supplementary material.)

The crystal growth velocity can be estimated by the following equation:

where a₀ is the lattice constant of Si, n is the number of atoms per unit cell (n = 8 for the diamond-cubic crystal structure), and A is the cross section area of the solid-liquid interface.

Following the above procedure, the growth velocities v(Δx_Si, T) are evaluated for simulations at 900, 950, and 1000 K. The velocities and their corresponding temperature and Si supersaturation are plotted in Fig. 3 for both $〈110〉$ and $〈111〉$ growths, to examine the relationship among the crystal growth velocity, temperature, Si supersaturation, and orientation.

For $〈110〉$ growth, the data points in Fig. 3(a) show that at a given temperature, the growth velocity increases almost linearly with Si supersaturation Δx_Si. When keeping Δx_Si a constant, the growth velocity has an Arrhenius-like dependence on temperature. (See Fig. S3 in the supplementary material.) Therefore, a function form described by Eq. (2) is proposed to fit the simulation data

Setting the values of the two parameters to a = 34.13 m/s and b = 0.83 eV, the fitting curves match well with the simulation data in Fig. 3(a).

For $〈111〉$ growth, observing that the crystal growth velocity has a clearly non-linear relationship with Si supersaturation, we adopt a function form proposed by Brice,²⁶ as expressed in Eq. (3). (The validation of this choice will be given in Sec. III B.)

Setting a = 2.30 × 10⁴ m/s, b = 1.26 eV, and c = 16.81, this equation generates the fitting curves in Fig. 3(b), which match the v-T-Δx_Si behavior predicted by the MD simulations reasonably well.

Based on the classical theory of crystal growth,^26–28 our velocity data for $〈110〉$ growth reflect a continuous growth mechanism on a rough surface, while the velocity data from the {111} simulations should correspond to an island growth mechanism on a smooth surface.

This substrate orientation dependence of Si growth can be verified by identifying the difference in the growth interface morphology between the {110} and {111} substrates. Figure 4 compares the final configurations for simulations performed at 1000 K and initial Si fraction of 50%. In Fig. 4, the grown atomic layers are distinguished from the initial substrate region by a white dashed line. As the color of the atom is assigned based on its height, the atoms which have the same color are within the same plane. Figure 4 also provides the horizontal cross sectional views, revealing the atomic arrangements at different substrate heights. For both {110} and {111} simulations, although some defects start to occur at high temperature and large Δx_Si, the grown layers inside the solid are single crystalline and relatively completed. At the growth interface, it is noteworthy that the atoms on the {110} substrate arrange themselves randomly, leading to a rough surface, while the atoms on the {111} substrate tend to cluster into flat islands. These identified features of the interface morphologies for {110} and {111} substrates are consistent with our knowledge from the classical crystal growth theories.^26–28 

If the island growth mechanism applies to {111} growth, the nucleation process is expected to be a critical step for initializing the growth, especially when the growth driving force is small. In order to capture the nucleation process and further confirm the crystal growth mechanism suggested above, a series of MD simulations are performed at 850 K, with small initial Si supersaturation Δx_Si of 0.5%, 1.0%, 1.5%, 2.0%, and 2.5%, respectively, for both {110} and {111} substrates. Figure 5 presents the final interface morphologies of the (101) and (111) substrates after 100 ns for each Δx_Si. Removing all the atoms in liquid and coloring the substrate atoms grey, the Si atoms deposited during the simulation are highlighted with orange. Our simulations show that with a very small Δx_Si of 0.5%, a reasonable amount of atoms have already been able to grow on the {110} substrate, in contrast to the {111} substrate where no obvious growth is observed until Δx_Si reaches 2.0%. This clearly suggests that the $〈111〉$ growth follows the island growth mechanism, where the nucleation rate is low at small supersaturation. Additionally, compared with Fig. 4, Fig. 5 provides a more clear view of the growth mechanism: the atoms grown on the {110} substrate are rather scattered, roughening the growth interface; on the other hand, the adatoms show a preference to form flat islands on the smooth {111} surface.

Figure 6 shows a series of simulation snapshots at different times for growth on the {111} substrate at 850 K with initial Δx_Si of 2%, in order to clarify the dynamic process of the growth. Here, only the crystalline atoms at the interface are displayed, and the atoms within the same cluster are assigned the same color. From the fact that new nuclei appear on every simulation snapshot of Fig. 6, we deduce that the nucleation event occurs rather frequently in this simulation. It is interesting to note that the large island observed at the end of the simulation is actually formed by the merging of several small islands, instead of expanding from a single nucleus. This indicates that compared with the nucleation rate, the island growth speed is lower. That is to say, the $〈111〉$ growth in our MD simulations is expected to be in the regime where small islands form rapidly but large islands expand slowly. This growth pattern matches the description of the polynuclear growth model, which has been discussed in the literature.^29–32 Based on this model, Eq. (3) was derived²⁶ and used in this study, to fit v-T-Δx_Si predicted by the {111} simulation data.

Understandably, the Si growth requires the separation of Si and Au atoms in the liquid. This could lead to a temporary increase of the local concentration of Au atoms near the {111} growth interface. Thus, the Si atoms in the liquid may be blocked from attaching to the newly formed island, resulting in a low island expansion rate. Figure 7(a) supports this hypothesis by showing the atomic arrangement in liquid close to the {111} surface. Interestingly, in Fig. 7(a), a mono-layer of Au atoms, represented by opaque yellow spheres, is formed on the crystalline Si cluster (The Si atoms in the substrate and the island on top are colored grey and orange, respectively). These Au atoms are periodically arranged, in contrast to the randomly distributed atoms (set to semi-transparent) away from the island. In our simulation, sometimes the gold mono-layer could survive for a relatively long time (several ns). This meta-stability is consistent with the fact that the Au thin film has good wettability on the Si substrate,³³ and the Si surface is an effective gettering site for Au atoms.^34,35 However, the formation of this gold mono-layer is not found in any of our {110} growth simulations. As shown in Fig. 7(b), the Si and Au atoms in liquid are well mixed upon the {110} substrate, without a significant increase of the local Au concentration. This suggests that these interface atoms are more mobile near the {110} interface. The different behavior of the atoms on {111} and {110} substrates may be due to their orientation dependent surface atom activity, which will be discussed in detail in Sec. III C.

In Eqs. (2) and (3), the fitting parameter b determines the temperature dependence of the growth, which can be interpreted as the activation energy Q_g for Si crystal growth. Apparently, the Q_g for $〈111〉$ (1.26 eV) is higher than Q_g for $〈110〉$ growth (0.83 eV). Also, under the same temperature and Si fraction condition, the $〈111〉$ growth is slower than the growth along $〈110〉$ orientation. We believe that the origin of these observations could be found by quantifying the interface atom activity.

As discussed in Sec. III A, the change of the order parameter q₃ for individual atoms can be used to characterize their activities during the simulation. The fluctuation of q₃ reflects the vibration of the atom in real-space. Note that q₃ = –0.5 is the criterion to distinguish solid from liquid, a jump of q₃ from −1 to 0 indicates an atomic transition event from the solid region to liquid region. Thus, the activities of the surface atoms can be quantified by monitoring their q₃ value as a function of time.

To eliminate the possible effect from the Si supersaturation on the interface atom activity, a set of MD simulations are performed at the liquidus composition, where the net crystal growth is zero. As a validation, Fig. S4 (see supplementary material) plots the number of crystalline Si atoms N_Si as a function of time for simulations at 1000 K. It shows that the N_Si fluctuates around its initial value for both {110} and {111} simulations, confirming that the system is in an equilibrium state.

As measures of the atom activities at 900, 950, and 1000 K, Fig. 8(a) provides the variations of q₃ with time for 50 atoms randomly picked from 397 Si atoms on the {110} interface that are initially in the crystalline state. Similarly, Fig. 8(b) plots the q₃ fluctuations of all the 1013 atoms on the {111} interface. In Fig. 8, at the same temperature, many more jumps of q₃ take place on the {110} surface than on the {111} surface, which clearly suggests that the atoms on the {110} surface are more active. Besides, for both $〈110〉$ and $〈111〉$ orientations, the jumps of q₃ become more frequent at higher temperature. This strongly indicates the positive correlation between the atom activity and the temperature. It should also be mentioned that here we only monitor the surface atoms initially recognized in the crystalline state. During the simulation, it is highly possible that when one surface atom leaves its lattice position, another atom in liquid may occupy the vacancy and become crystalline, so that the original atom is not able to go back to the substrate. This may explain the lack of q₃ jumping from 0 to −1 in Fig. 8.

To further quantify the surface atom activity, we define a transition rate ν, based on our observations from Figs. 8(a) and 8(b). Since every atom on the plots is crystalline at time 0, if an atom is found in the liquid region at time t, at least one solid-to-liquid transition should happen to this atom during this period. To determine the state of an atom at t, the time average $〈q3〉$ for a time span of 5 ns is calculated. If $|〈q3〉|<0.1$⁠, the atom is identified as in the liquid phase. Thus, the transition rate ν of a simulation is defined as the ratio of the total number of the transitions (the number of liquid atoms at t) to the simulation time t per unit area. For $〈111〉$ simulations, from Fig. 8(b), we can observe that one atom has at most one transition. Therefore, the calculation is performed at t = 100 ns to include as many transition events as possible. For $〈110〉$ simulations, since the back-and-forth jumps between 0 and 1 occur on the q₃ fluctuation curves in Fig. 8(b), the $〈110〉$ transition rate is obtained by taking the average of the calculations at t = 30, 40, and 50 ns, respectively, in order to improve the accuracy.

Figure 9 plots the logarithm of the normalized transition rate as a function of the inverse of temperature, where ν₀ is the transition rate at 900 K. For {111}, $ν0111=0.001(nm2·ns)−1$⁠, and for {110}, $ν0110=0.012 (nm2·ns)−1$⁠, roughly 12 times larger than $v0111$⁠. This is consistent with our observations in Fig. 3 that given the same growth condition, the {110} growth is faster. In Fig. 9, the relationship between $log (ν/ν0)$ and 1/T can be well described by a linear fitting curve for both {110} and {111} simulations. This suggests that the transition rate's temperature dependence follows the Arrhenius relationship, $v∝ exp (−QakBT)$⁠. The activation energy Q_a can be estimated from the slope of the fitting curve. For {111}, $Qa111=1.28 eV$⁠, and for {110}, $Qa110=0.74 eV$⁠.

It is of interest to compare Q_a (activation energy of surface atom activity) and Q_g (activation energy of growth) with the activity energy Q_d for diffusion in liquid. We note that diffusion, which transports the Si atoms to the interface, is a necessary step in VLS growth. To investigate the temperature dependence of Si diffusivity, separate simulations of the Si-Au binary alloy are performed, and the value of Q_d is obtained in Table II (Simulation details are given in Sec. S5 in the supplementary material).

Table II shows that the crystal growth velocity's dependence on temperature is positively correlated with both the surface atom activity and the diffusivity in liquid, i.e., all rates become larger at higher temperature. However, since the liquid diffusivity of Si does not have an orientation dependence, the diffusivity alone cannot explain the difference of Q_g between {110} and {111} surfaces. In addition, Q_d is significantly lower than Q_g for {111} and {110} substrates, which implies that the strong temperature dependence of the Si crystal growth must come from somewhere else. The diffusivity of Au is also estimated (Sec. S5 in the supplementary material), whose value is close to that of Si in their liquid alloy for all temperature ranges considered here. This means that for the Si-Au system, when far away from the interface, the Si diffusion in liquid should not be obstructed by the Au atoms. All these findings suggest that the diffusion in liquid is not the rate limiting step in VLS growth.³⁶ At the same time, Q_a for the {111} surface is larger than that for the {110} surface, and is largely consistent with Q_g. This indicates a strong correlation between the growth velocity and the interface atom activity for the Si-Au system.

To systematically study the Au-catalyzed Si crystal growth via the VLS mechanism, we perform extensive long MD simulations under conditions of zero, low, and high supersaturation in the liquid. The velocities for growth on {110} and {111} substrates are extracted from the simulations at various temperatures and Si fractions, based on which the fitting functions are successfully proposed to describe the relationships among the growth velocity, temperature, and Si supersaturation. Different patterns for growth on different oriented substrates are recognized: continuous growth on the rough {110} substrate, and island nucleation on the smooth {111} substrate with the island expansion rate slower than the nucleation rate. Additionally, the interface atom activity is quantified, and found to be related to the surface morphology. This helps explain both the temperature and orientation dependence of the Si VLS growth. The findings from these simulations improve our understanding on the fundamental mechanisms of the VLS process for the Si-Au system, and should inspire atomistic studies on other VLS systems, e.g., the Au-catalyzed Ge growth.³⁷ 

See the supplementary material for the diamond cubic structure identification, calculation of the Si fraction in liquid, raw simulation data showing the temperature dependence of Si growth velocity, MD simulation data at liquidus composition, and calculation of Si and Au diffusivities in their liquid alloy.

Adriano Santana would like to acknowledge the computational support from the Beijing Computational Science Research Center (CSRC). The computer simulations are performed on Sherlock cluster operated by the Stanford Research Computing Center, and Tianhe2JK cluster operated by CSRC. Yanming Wang and Wei Cai acknowledge support from the National Science Foundation Division of Materials Research Program No. DMR-1206511.