Atomistic Modelling of Si Nanoparticles Synthesis

Silicon remains the most important material for electronic technology. Presently, some efforts are focused on the use of Si nanoparticles—not only for saving material, but also for improving the efficiency of optical and electronic devices, for instance, in the case of solar cells coated with a film of Si nanoparticles. The synthesis by a bottom-up approach based on condensation from low temperature plasma is a promising technique for the massive production of such nanoparticles, but the knowledge of the basic processes occurring at the atomistic level is still very limited. In this perspective, numerical simulations can provide fundamental information of the nucleation and growth mechanisms ruling the bottom-up formation of Si nanoclusters. We propose to model the low temperature plasma by classical molecular dynamics by using the reactive force field (ReaxFF) proposed by van Duin, which can properly describe bond forming and breaking. In our approach, first-principles quantum calculations are used on a set of small Si clusters in order to collect all the necessary energetic and structural information to optimize the parameters of the reactive force-field for the present application. We describe in detail the procedure used for the determination of the force field and the following molecular dynamics simulations of model systems of Si gas at temperatures in the range 2000–3000 K. The results of the dynamics provide valuable information on nucleation rate, nanoparticle size distribution, and growth rate that are the basic quantities for developing a following mesoscale model.


Introduction
The synthesis of nanoparticles in the gas phase by a bottom-up process is emerging as a promising technology compared to traditional production techniques.This is because, even though the wet phase synthesis has a better control of the nanoparticle morphology, it requires solvents and cannot be easily used for industrial production.In contrast, the milling techniques are relatively simple and can be adapted to large scale productions, but do not provide a good control of size and morphology of the nanoparticles.Furthermore, some bottom-up processes do not use solvents and are therefore more eco-friendly.In the case of hard materials that are stable at high temperatures, the synthesis from gaseous phase in plasma reactors [1], flame reactors [2], and hot-wall reactors [3] is especially interesting.The operating conditions-in terms of density and high temperature (2000-3000 K) of the gas-in these types of reactors make the environment in which the synthesis occurs extremely complex, difficult to monitor and investigate with experimental techniques.This is due to the non-equilibrium thermodynamic processes, the kinetics that may develop in the range of ns or even of ps, and the chains of chemical reactions between the precursors and any contaminant.
For such complex systems, computational modelling might in principle provide valuable insights into the basic mechanisms that affect the synthesis process.These are, for example, nucleation, growth of the embryo-particles, and their surface chemical activity.Of course, numerical modelling has its objective difficulties too, due to the size of the models and the time scale of the simulations required to obtain results comparable to experiments.The basic approach in this case is to define-when possible-different modeling regions of spatial and temporal dimensions that are contiguous and partially superimposed.In the specific case of the plasma synthesis, the modelling range extends from the atomic size of the primary species up to the macroscopic continuous reactor and from the femtoseconds of molecular dynamics up to the milliseconds of diffusion and turbulence.Chemical engineers have developed and already employ numerical models on the macroscopic scale of the reactors, based on empirical parameters which, however, are not easily measurable (as in the case of the reactors here considered).This has led to the development of meso-scale modelling for a better understanding of the systems at a microscopic level.However, such models also depend on a set of basic parameters derived from thermodynamical data obtained on bulk material that can generally have properties very different from that of a nanoparticle.Atomistic modelling can provide insight about the basic processes that occur on the scale of nanometers and picoseconds, and then drive the subsequent models on larger scales to a better description of the system.
Atomistic modelling is the subject of the present report concerning the bottom-up synthesis of nanoparticles from Si low temperature plasma.The main goal is to estimate-by means of first-principle and molecular mechanics calculations with force fields (FF)-the morphology of embryo-particles (small size clusters) that form through homogeneous nucleation of the monoatomic gas and their growth during the first steps of the nucleation process.The numerical procedures described here can be easily extended to other hard materials.They are based on density functional theory (DFT), basin hopping [4] (DFT-BH) calculations, and molecular dynamics simulations based on reactive force fields.These force fields can describe bond breaking and formation during the evolution of the systems, which is fundamental in the cases examined in this study.
After a brief overview of the features of the reactive force field [5] and of the procedure to optimize its parameters for the case of Si nanoparticles, the results of the dynamics in terms of inception time of the nanoparticles as a function of their size and of the environmental conditions (temperature and density of the primary gas) will be presented.These results will be critically discussed and compared to the predictions provided by the Classical Nucleation Theory (CNT) [6], which, in different forms, is largely used to describe homogeneous nucleation.

Model and Computational Details
The atomistic model adopted to simulate the early stages of the gas phase condensation of Si nanoparticles consists of classical molecular dynamics simulations based on a reactive force field (ReaxFF) originally developed by van Duin [5], where the parameters were recently optimized [7] for Si-containing systems.ReaxFF is flexible, computationally efficient, and has been successfully applied to a large variety of material modelling problems [8][9][10][11][12][13][14][15].For these reasons, it has been selected to model the early stages of the plasma-assisted growth of Si nanostructures.
For the simulation of this kind of process in Si-containing materials, two formally simpler potentials-namely, Tersoff and Stillinger-Weber potentials-have been frequently employed.Among the most recent applications there are the computational studies for the simulation of the condensation, growth, and crystallization of Si nanoclusters obtained by magnetron sputtering [16][17][18].In this technique, the fragments of Si (mainly embryo clusters but also isolated atoms) produced by sputtering at very high temperatures (even 6000 K) pass through a condensation region occupied by an inert gas (Ar) at a much lower temperature.The sample is obviously very far from equilibrium, and has been simulated by molecular dynamics for an ensemble with constant number of atoms, volume and total energy (NVE) with a time gradient of the temperature that is taken as large as 30 K per ns.In our case we intend to refer to environmental conditions much closer to equilibrium and then better represented by simulations for an ensemble with constant number of atoms, volume and average temperature (NVT) with large fluctuations induced by the Berendsen's thermostat adopted.These conditions correspond locally to those that are present in the thermal plasma reactors of industrial type employed for large-scale production of nanoparticles.
ReaxFF is expressed by an extended set of analytical functions that reflect the complex picture of the atomic interactions and depend on a number of different parameters.These parameters have been derived, in the present case, from quantum chemistry calculations on model systems representative of the species present in the Si plasma.The complex analytic form of the force field is, anyway, much less computationally expensive than the Hamiltonian formulation of the quantum molecular dynamics (QMD) methods, and as a consequence, classical simulations are the only practical choice for describing large aggregates containing up to one million atoms.In comparison to most of the force fields available in the literature, which essentially provide a physical picture of the system, ReaxFF allows the simulation of bond breaking and formation.This aspect is very relevant in order to properly simulate not only interactions, but also reactions, and then dynamically characterize the various processes that determine and regulate the early formation of nanomaterial aggregates from the plasma.The potential energy in ReaxFF is calculated as the sum of a number of energy terms depending on the bond distance/bond order and bond order/bond energy relationships.The algorithm appropriately describes under-and over-coordination and dismisses the energy terms depending on bond distance, bond, and dihedral angles upon bond breaking.Non-bonded interactions (Coulomb and van der Waals forces) are calculated between atom pairs with a shielding factor [19] to dampen close-range interactions; atomic charges are calculated using the electronegativity equalization method (EEM) [20] and long-range Coulomb interactions are calculated by means of a taper function with an outer cut-off radius.

Optimization of the ReaxFF Parameters
Starting from a parametrization of the FF for Si based materials (VD) available in the literature [7] and designed to describe the bulk structure of SiC and SiO 2 , we re-optimized the FF for the specific systems under study in order to consider surface effects of small Si clusters.These effects can be fundamental in the nucleation and growth processes, leading to nanoparticles with maximum diameter of about 1 nm.The FF parameters have been optimized in order to reproduce a training-set (t-set) of properties (i.e., energies and geometries) derived from quantum chemical (QC) calculations on selected model systems.Few experimental measurements of dissociation/ionization energies and mobilities of the ionized structures are available in the literature [21][22][23][24], whereas geometries and electronic structures of small silicon clusters have been studied extensively by QC approaches over the last few years [25][26][27][28][29][30][31][32][33][34][35][36][37][38].
The first-principle quantum method we have adopted is based on DFT [39], which is widely used for studying medium-large molecular systems.Calculations have been performed by the ESPRESSO package [40] using ultrasoft pseudo potentials [41] for an approximate description of the core electrons that do not directly participate in the formation of molecular bonds, whereas the electronic structure of the valence electrons (4 for Si) is described by the GGA-PBE [42] functional.Energy spin-restricted optimizations were performed in a cubic box with side between 10 an 20 Å (depending on the size of the cluster) and periodic boundary conditions in order to effectively project the wave function on a set of plane waves with cut-offs on the wave-function/electronic density of about 140/1100 eV, respectively.Other details of the DFT calculations are: gaussian smearing for the occupation of the Kohn-Sham single-particle states of 0.03 eV; Brillouin zone explored at the Gamma point only; single local minimization with a threshold on energy/forces of 10 −3 /10 −4 au.
A small Si cluster may adopt various conformations corresponding to slightly different energy values; at the plasma temperatures (>2000 K), a great number of such conformations is accessible.For the optimization of the force field, we selected the most populated minimum energy configurations of the clusters that have been identified by a specific algorithm based on geometry optimization and Basin-Hopping [4].This is one of the several global optimization methods that are used to find the minimum of a multi-variable Potential Energy Surface (PES).As any other similar method, this approach does not guarantee a successful identification of the absolute energy minimum in a finite time, but provides an efficient algorithm for exploring the PES and accumulating a voluminous structural database.In our approach, for each cluster size, we performed a search for the minimum by running several BH calculations (from three to five), composed of 500 steps each.At each step, the cluster geometry is slightly modified randomly and a DFT geometry optimization is performed; the step is accepted or refused according to the Monte Carlo/Metropolis criterion [43] with a kT value in the range 0.5-1.0eV [44].The lowest-energy structures resulting from this search are summarized in Figure 1 for Si clusters composed of an even number of atoms with a size in the range 6 to 36.These structures were considered as model systems for computing the benchmark data forming the training set.To refine the VD FF, we adopted the ReaxFF standalone code (version 2.0).The t-set adopted was given by geometries and energies of a number of model systems selected among the structures optimized by the DFT-BH searches, namely: Si 2 , Si 3 , Si 4 , Si 5 , Si 6 , Si 8 , Si 10 , Si 12 (two structures), see Figure 2, to which the basic unit of silicon diamond cubic crystal was added.Starting from the original parametrization (VD) [7], we optimized only the relevant parameters responsible for the flexibility of the Si clusters at high temperature against our t-set: (i) Si-Si bond (14 parameters); (ii) Si-Si-Si valence angle (5 parameters); (iii) Si-Si-Si-Si torsion angle (5 parameters).These latest parameters were not present in the original parametrization, but were added in this new force field to get a better agreement with the t-set data.Starting values of the parameters were taken from carbon in VD.
The ReaxFF code adopts a sequential optimization procedure of the parameters by the minimization of an error function (EF) which is calculated in the following way: where f FF,i and f QC,i are, respectively, the running FF and the benchmark QC values of the i-th entry of the t-set and w(i) is its weight.The sum is extended to all the geometrical/energetic entries of the t-set.
The ReaxFF optimization procedure works in the following way: at each step, a single parameter is considered and its starting value is varied by adding and subtracting a given quantity; for these three values of the considered parameter, the total EF is estimated and a parabolic extrapolation is achieved by finding (in correspondence of the minimum of the parabolic expression) the optimal value of the parameter within the chosen window; if this new predicted value corresponds to a lower value of the EF, the change is accepted and the algorithm moves to the successive parameter.As the parameters are strongly entangled, several sequential optimization cycles have to be repeated in order to achieve a real reduction of the EF against the t-set.It should be added that at each step the QC geometries optimized by the DFT-BH procedure are compared to those optimized by using the running FF; this is done in view of using the FF for an effective and realistic search of the different stable conformations of the clusters.From the point of view of this performance it is convenient to analyze the predicted lowest-energy structures by using the two synthetic structural descriptors: R and θ.The first one refers to an average value of the interatomic distances: where N is the number of atoms in the cluster, each i−atom has n i first-neighbours, and r i,j is the distance between atoms i and j; the second one refers to an average value of the molecular angles:

Reactive Molecular Dynamics
Classical Molecular Dynamics simulations of the Si + Ar plasma have been carried out by using the LAMMPS code [45] with the presently optimized reactive force field (in the following named ND) in the NVT ensemble with the Berendsen's thermostat in order to have temperature fluctuations around a constant average temperature.In our approach, such simulations-performed at different compositions and temperatures-intend to model the environment present in thermal plasma reactors employed for large-scale production of Si nanoparticles.The results are employed to investigate nucleation processes and estimate specific parameters and mathematical expressions that can be used to feed mesoscale models for the simulation of much larger systems.
On the basis of the good performance of the ND FF on the exploration of the PES at 0 K, extensive MD simulations have been carried out by using this FF; nevertheless, a comparison with VD on a specific set of thermodynamic conditions will be discussed in the next section as well, in order to show how the two different parameterizations influence the nucleation and growth processes.
The computational details of the nucleation simulations are the following: • Simulation box: an appropriate value of the length of the cubic simulation box has been chosen in order to reduce the effect of the periodic boundary conditions.The number of events that occur on the border becomes less and less relevant than the number of total events when the size of the box is increased.In the limit of a box with an infinite size, the biasing constraint of the translational periodicity disappears.According to a series of simulations in which the volume of a cubic unit cell has been progressively increased (from a minimum of 64 ×10 3 nm 3 to a maximum of 13824 × 10 3 nm 3 ), the volume chosen was 4096 ×10 3 nm 3 .• Temperature: each simulation has been carried out at fixed average temperature.The temperature range between 1800 and 2400 K has been explored by increasing the temperature by steps of 100 K in each simulation.• Atom density: the number of atoms within the simulation box has been chosen on the basis of particular experimental conditions: total gas pressure within the plasma reactor around 1 atmosphere (standard ambient pressure).On the basis of this latter value, the density of the simulation gas (which is composed of Si and Ar atoms) is about 3.6 × 10 24 atoms/m 3 , which corresponds-in the chosen cubic box-to an ensemble of about 12,000 atoms.• Composition of the mixture: the content of the (Si, Ar) mixture has been varied according to the following percentages: (i) 25% Si and 75% Ar; (ii) 50% Si and 50% Ar; and (iii) 75% Si and 25% Ar. • Time step: The time step of the MD simulations has been set to 2 fs, having verified that a smaller time step of 0.5 fs reproduces the same results at both the qualitative and quantitative levels.
We underline that the adopted models refer to environmental conditions close to equilibrium corresponding, locally, to those that are present in the thermal plasma reactors of industrial type for large-scale production of nanoparticles.

Results and Discussion
From the point of view of assessing the performance of the reactive FF in exploring the PES of the small Si clusters, it is interesting to compare the structures predicted by the VD and the presently refined parametrization ND of the reactive FF against DFT results.In Figure 3, the lowest-energy structures predicted at sizes 12, 20, 30, and 36 are reported, together with a structural and energetic analysis of the conformations obtained.It should be noted that the lowest-energy structures shown in Figure 3 have been obtained by adopting the same BH search protocol described above, with energies and forces calculated at the FF level, but the energy of such structures has been then computed at the DFT level for a correct comparison with the DFT-BH results.
From the structural point of view, we can observe two main features that characterize the geometries of small Si clusters: (i) a compact arrangement of atoms with coordination numbers and coordination angles much different from bulk silicon, where the coordination number is strictly four and the angles are perfectly tetrahedral; in the DFT-BH optimized structures, instead, the number of nearest neighbours can be as large as 6 and even 7, and the bond angles are remarkably different from tetrahedral; (ii) cluster surfaces are not smooth, but sharp with the presence of convex re-entrances, leaving the surface atoms with a very low coordination of one, or maximum two, neighbours.2) and ( 3), respectively, while σ refers to their root-mean-square dispersion.DFT: density functional theory; ND: present optimization of ReaxFF; VD: parametrization for Si-based materials available in the literature [7].
By inspecting the results of the two FFs, we can see that the presently optimized parametrization ND using the t-set shown in Figure 1 better reproduces both the structural and energetic DFT results.Indeed, the two structural features discussed above can be found in the geometries predicted by the ND FF, and the relative FF energies are closer-compared to the VD values-to the benchmark DFT ones.The limitations exhibited in the present application by the VD FF could be due to its apparent tendency to favor Si atoms with "bulk-character".As a consequence, in the structures predicted for the small clusters, the average coordination number is about four, as evidenced by the high number of red atoms in Figure 3, which are tetrahedrally shaped, and by very smooth surfaces; moreover, the DFT energy of the VD FF optimized structures is higher than that predicted by the ND FF [44].
A direct comparison of our present MD results with those available in the literature [16][17][18]-using the Tersoff or Stillinger-Weber potentials in order to evaluate the performance of different force fields-would be biased by the large differences in the model systems adopted and in the characteristics of the MDs employed.An accurate comparison would require calculations in which the features of the simulations are more homogeneous, which is outside the scope of the present work.
Two typical results of the nucleation MD simulations at different temperatures are shown in Figure 4. On a short time scale, the dynamics is dominated by the formation of smaller aggregates (mainly Si dimers, which are quickly formed); the concentration of the smaller clusters reaches a plateau and then starts decreasing due to the formation of the larger aggregates.As can be noted in Figure 4a, a further peak is found at size 8, which corresponds to a remarkably stable structure of the Si system.As the simulation proceeds, larger aggregates start to form.At lower temperature (Figure 4a), several Si 30 clusters can be found at the end of the run (about 300 ns), whereas at higher temperature (Figure 4b), the particles growth is slowed down due to the larger thermal agitation and only few aggregates of size around 30 can be obtained at the end of the simulation.This picture is different from that provided by the classical nucleation theory, in which, out of nothing, a nanoparticle of critical mass is produced without a previous history.In our opinion, this phenomenological model does not apply to the cases where-as for Si-the interaction between atomic components leads to the formation of covalent bonds with bond energy much higher than the available average thermal energy.In Table 1, the characteristic inception times for particles of size 20, 30, and 40 are reported as a function of the temperature for the three mixture compositions specified above.The inception time of a particle composed of N atoms is estimated as the time when the first particle of that size appears in the simulation box.For a specific nuclearity, the inception time increases when increasing the temperature due to the higher thermal agitation that slows down the aggregation process.As could be expected, the inception time decreases at increasing Si content when comparing the same particle size and temperature.The Ar gas present in industrial thermal plasma reactors provides a more effective coupling with the electromagnetic field applied for the formation and maintenance of the plasma; the case of a pure Si gas is practically not realized.Our simulations (not included in this study) of such a system show that the processes of condensation and growth of Si nanoparticles remain substantially unchanged; the presence of Ar is not necessary for their implementation.When using the force field VD, the results of the nucleation and growth processes are completely different: at T = 2000 K with a plasma made of 50% Si and 50% Ar (within a time scale of around 160 ns), only a very low percentage of atoms form aggregates, as it can be evinced by inspecting Figure 5a; moreover, only dimers are formed (about 50 during the simulation, see Figure 5b) with the sporadic appearance of trimers, which soon dissociates again in dimers and single atoms (see Figure 5c).This behaviour could be due to an instability of the trimer that originates a bottleneck in the growth process [44], as this nuclearity is a mandatory step in the growth of larger aggregates, due to the fact that the formation of a tetramer due to the aggreagtion of two dimers is a very unlikely event in the present conditions.To conclude, we compare the results of the present MD simulations with the ND FF and the only-to our knowledge-work in the literature reporting the values of partial pressures of small Si N clusters (N = 1-6) in the vapour phase in equilibrium with the Si liquid phase calculated on the basis of classical thermodynamics [46].In order to compare with the reported data, in Figure 6 we have plotted the ratio between the numbers of dimers and single atoms vs. the ratio between the number of trimers and dimers as a function of the simulation time.In this case, an NVT dynamics calculation has been performed at a total pressure of 1.22 × 10 −2 MPa (1 atm = 0.1 MPa) and a temperature of T = 3000 K, matching one of the conditions considered in the experiment.The choice to operate at a higher temperature (with respect to the previous simulations) has been done to work with a higher value of Si pressure, and hence collect more statistical information about the investigated system.The simulation box contained 2860 atoms, and time-step has been set to 2.0 fs; the Berendsen's thermostat has been employed, in agreement with the methodological procedure applied for the estimation of the nucleation rates.According to the results reported in Figure 6, it can be seen that the dimers/atoms ratio is quite stable and reaches a plateau around the value 0.05; on the other side, the trimers/dimers ratio is less regular and is characterized by remarkable oscillations, mainly due to statistical reasons (lower absolute number of trimers in the simulation box during the dynamics); nevertheless, we can estimate that the latter ratio oscillates around the value 0.10 during the simulation.The obtained values of the ratios are in qualitative agreement with the data of Reference [46], which reports a dimers/atoms ratio of about 0.11 and a trimers/dimers ratio of about 0.25 at the considered values of temperature and pressure.

Conclusions
We have proposed an theoretical approach to the difficult problem of modelling the nucleation and growth of embryo-particles in a low temperature plasma, with application to the Si gas.The approach is based on reactive molecular dynamics with a FF specifically optimized using DFT calculations as benchmark for a set of small Si clusters.The presently proposed (ND) FF-optimized starting from the one (VD) available in the literature-shows appreciable improvements in the description of Si nanoparticles.The results obtained by exploring the PES at 0 K evidenced the tendency of VD FF to over-estimate the bulk behavior in Si clusters.On the contrary, ND FF was capable of more confidently reproducing small clusters and favoring the formation of compact aggregates.The results of the molecular dynamics simulations performed with the ND FF are in reasonable agreement with the very few experimental data available for this kind of system, although the experimental environment corresponding to Si vapor in equilibrium with the liquid phase is only roughly represented by our simulation model.Unfortunately, an extension of the model including a representation of the liquid phase would lead to too-CPU-demanding calculations.Our results can be viewed as describing a state of quasi-equilibrium only in the sense that the obtained ratios are quite stable in a wide time window; consequently, the agreement of theory and experiment can be expected only at a qualitative level.The results obtained show that the proposed theoretical model can be satisfactorily applied to systems where the experimental studies are very difficult to perform and an accurate description of the basic processes at the atomistic level is presently missing.

Figure 1 .
Figure 1.Lowest-energy structures of Si clusters with an even number of atoms and size in the range 6-34.

Figure 2 .
Figure 2. Si clusters included in the training set used in the FF optimization.

Figure 3 .
Figure 3. Structural and energetic results obtained by using different parametrizations of ReaxFF for Si clusters of size 12, 20, 30, and 36.The structural data R and θ are defined in Equations (2) and (3), respectively, while σ refers to their root-mean-square dispersion.DFT: density functional theory; ND: present optimization of ReaxFF; VD: parametrization for Si-based materials available in the literature[7].

Figure 4 .
Figure 4. Nucleation and growth picture at (a) 1800 K and (b) 2400 K.

Figure 5 .
Figure 5. Number of (a) unreacted Si atoms; (b) dimers; and (c) trimers as a function of time at T = 2000 K.The results from the ND and VD FFs are shown in red and green, respectively.

Figure 6 .
Figure 6.Ratio between the number of dimers and unreacted atoms (in red) vs. ratio between trimers and dimers (in green) for a set of thermodynamic variables specified in the main text.