Modelling of nanoparticle coagulation and transport dynamics in dusty silane discharges

This paper reports on a self-consistent one-dimensional (1D) hydrodynamic model that investigates the formation, the subsequent growth and transport of nanoparticles in a parallel plate capacitively-coupled radio-frequency silane (SiH4) discharge. A fully coupled description of the first two stages of particle formation (nucleation and coagulation) is attained by the development of a model which treats the electron kinetics, the gas phase chemistry, the particle formation mechanisms, the nanoparticle charging and transport dynamics, and the coagulation phase together with a self-consistent determination of the plasma properties. In the present paper, we focus on the fast coagulation stage, incorporated by making a self-consistent coupling between the 1D fluid model and an aerosol dynamics model in which an evolution of the nanoparticle size domain is obtained by utilizing a sectional approach to solve the general dynamic equation (GDE). During each coagulation step, the effect of nanoparticle charging and transport is included and solved with the same temporal resolution. The calculated density and charge distribution profiles are presented for particles ranging in size between ∼1 and 50 nm. The concerted action of particle charging and transport is found to severely affect the location of nanoparticle growth due to coagulation. Heating of one of the electrodes immediately induces a thermophoretic force that can be considered as a useful means to control particle contamination.


Introduction
Because of the ubiquitous presence of dust particles in every plasma environment, complex plasma systems or dusty plasmas have received much attention in recent years and are an active research field in both space and under terrestrial conditions. Nano-and microscopic particles are encountered in numerous space systems, such as planetary rings, interstellar clouds, comet tails, and diffuse nebulae [1]- [3]. In terrestrial laboratory devices, dust formation has especially been an important issue in industrial plasma processing, where charged nanometre or micrometre sized dust grains appear as a result of chemical polymerization reactions in the ionized gas phase or from plasma-surface interactions [4,5].
In the microelectronics manufacturing process many production steps involve plasmas, where the formation of dust particles is generally deemed as an unwanted process. The particles can be deposited on the wafer during plasma operation and tend to cause voids, delamination and interconnect shorts that will influence device topography, performance and reliability [6]. Therefore, dust is often considered as a limiting factor, which can drastically degrade the product quality and reduce the manufacturing yield. As the typical feature size of the integrated circuitry elements becomes smaller, the presence of contaminants will become an even greater concern, since the critical killer particle size, defined as one-half of the gate length, is projected to go down below the 50 nm range in the next couple of years. 4 Notwithstanding this harmful aspect, the presence of fine dust particles can also be beneficial in certain material science applications, where the structural incorporation of nanoparticles, produced by the use of plasma technology, or trapped in the plasma after external injection, can be used for the manufacture of advanced nanomaterials and electronic devices [7,8].
Better performance of photovoltaic cells can for example be achieved by operating at plasma 3 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT conditions in which fine nanocrystalline silicon particles are formed and incorporated in the intrinsic amorphous hydrogenated silicon (a-Si:H) layer. These novel so-called polymorphous silicon (pm-Si:H) thin films in particular display superior electronic properties that can lead to a significant increase of the product quality. Therefore, these layers are considered as an interesting candidate for the production of high-efficiency solar cells that exhibit good stability against the Staebler-Wronski effect, which normally reduces the cell's efficiency under light exposure [9,10]. Recently, a new exciting research field adopts reactive silane-based plasmas for the gas phase synthesis of silicon quantum dots used for the fabrication of light-emittingdiodes (LEDs), ultrahigh-density optical memories and data storage, and chemical and biological sensing [11].
Hence, due to their various technological implications, small particles become increasingly important. To this effect, control of nanoparticle size and deposition rate is obviously crucial if one wants to optimize the above-mentioned applications. A proper selection of the reactor design and operating conditions can not only lead to an increase in growth rate and optimization of film properties, while still minimizing particle contamination in the semiconductor industry, but can also be used for the development of novel methods for the assembly of nanostructured materials.
In general three phases can be distinguished in gas phase particle formation [12,13]: 1. Nucleation, i.e. the gas species-to-particle conversion process which creates small particles (nuclei) that typically consist of hundreds of atoms and which is largely dominated by plasma chemistry. 2. Microscopic particle growth from a few nm to 50-60 nm by fast coagulation. 3. Further growth of microscopic particles by the attachment of radicals on the surface of existing particles until micrometre sized particles are obtained. Finally gravity will pull the particles out of the discharge.
Although the behaviour of particles in the micrometre regime is relatively well understood, so far a comprehensive understanding of the formation, growth, transport and charging mechanisms of the smaller nanoparticles still remains elusive, essentially due to the experimental difficulties in detecting and monitoring such small particles. Only a few experimental diagnostics are nonintrusive and especially small particles during the nucleation and early coagulation stages are difficult to detect, as the particle size is still below the detection limit of the applied diagnostic. In order to obtain more information on these initial particle generation phases, computer simulations are considered as a powerful tool that can help in clarifying the underlying reaction mechanisms and the processes occurring at these earlier stages. Previously [14,15], we developed a detailed chemical kinetics mechanism for gas phase nucleation of hydrogenated silicon particles that incorporates silicon hydrides containing up to 12 silicon atoms (Si n H m with n 12). A sensitivity analysis enabled us to determine the dominant reactions and species underlying the earliest phase of dust formation. Anion-induced polymerization reactions appeared to be the most prominent pathway to consider, as 90% of the silicon dust growth was triggered by the reactive SiH − 3 precursor in our simulations [14].
In the present study, we focus our interest on the second stage of particle formation, i.e. fast coagulation. Several experimental studies [13,16,17] have indicated that the agglomeration phase only starts when the number density of primary clusters exceeds a critical threshold, typically of the order of 10 10 -10 11 cm −3 . This condition is usually referred to as the onset of particle generation. During coagulation the particle number density quickly decreases over 4 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT several orders of magnitude, whereas the particle size strongly increases. Once the particle size exceeds 50-60 nm in diameter, further coagulation is prevented by Coulomb repulsion between the large negatively charged grains. After agglomeration the so-called α-γ transition occurs [17,18], which is characterized by a strong increase of electron temperature (from 2 to 8 eV) and a drop of the electron density due to the enhanced attachment of free electrons on the nanoparticle's surface.
In order to explain these experimental observations considerable theoretical effort has been made on the dynamics of particle coagulation in plasmas. Kim and Ikewaga [19] analysed the evolution of particle formation with a simple plasma chemical reaction scheme for a silane plasma reactor, but assumed a uniform electron number density and electron energy inside the reactor and neglected the time-varying effect of the electric field. The charging of particles was also not taken into account. Courteille et al [20] developed a simple model to describe silicon particle agglomeration by Brownian free molecular coagulation (i.e. particles collide by random Brownian motion in the regime where Kn 1 with Kn representing the Knudsen number defined as Kn = λ/r d where λ is the mean free path of gas molecules and r d the particle radius) based on the equality of the electron and ion currents and on the plasma charge neutrality. In their model, the agglomeration starts from a given radius and a given number density of particles and assumes a constant positive ion density, and electron and ion temperature. No negative ions have been taken into account, which are however important in electronegative silane plasmas. A comparison with their experimental measurements showed that at least the beginning of the particle coagulation phase can be reasonably well described with their developed model, although no charge fluctuations due to statistical variations in electron and ion fluxes to the nanoparticle surfaces have been taken into account. Schweigert and Schweigert [21] did treat the distribution function of the particle's charge during particle coagulation in a low temperature argon (Ar) plasma, but assumed a fixed electron temperature of 4 eV. From their results, it was found that most of the particles are charged negatively. Lee and Matsoukas [22] calculated the coagulation rates of charged particles analytically taking into account the statistical distribution of the particle charge and showed that these fluctuations can increase the rate of agglomeration, especially for smaller particles. Kortshagen and Bhandarkar [23] considered particle coagulation via the solution of the general dynamic equation (GDE) with the inclusion of the particle charge distribution and suggest that the positive ion density constitutes the critical density for particle agglomeration. The influence of previously neglected effects, such as UV photodetachment of electrons, electron detachment due to quenching of excited atoms at the particle surface and secondary electron emission, on the particle charge are also investigated.
However, despite the intensive investigation of earlier studies, one important aspect has practically not been covered so far, i.e. the influence of the dust particles on the discharge. Experimental observations nevertheless seem to warrant their explicit consideration as a significant change of plasma properties is observed during the coagulation stage, which also affects the rate of particle formation, charging, transport and growth. So far only the work of Kortshagen and Bhandarkar [23] incorporates the self-consistent determination of the plasma properties in their global (0D) model, but without any consideration of the particle generation routes. Their simulations start with a monodisperse distribution of primary particles of typically 2 nm in size and investigate the temporal evolution of particle coagulation at a single discharge position.
In the present study, we extend our previously developed one-dimensional (1D) fluid model that describes particle formation and growth by nucleation to include further growth 5 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT by agglomeration in a parallel-plate reactor via a coupling to an aerosol dynamics model. The agglomerate phase will be accounted for using a sectional approach, which is advantageous in that there are no a priori assumptions regarding the nature of the particle size distribution. An important extension compared to the model of Kortshagen and Bhandarkar [23] is the inclusion of the first stage of particle formation, i.e. the nucleation. By including the growth kinetics of the particles and even considering a separate electron energy distribution function (EEDF) module for the determination of the electron-neutral reaction rate and electron transport coefficients as a function of average electron energy, we are able to determine the absolute densities of primary nuclei self-consistently as a function of position in the plasma which then form the starting point of our coagulation module. Furthermore, in contrast to their global model, we are able to determine the location of particle growth by including a transport equation for every incorporated particulate. Hence, this also allows us to observe the effect of forces such as neutral drag, ion drag, thermophoresis and gravity, which besides the electrostatic force can play a dominant role in the transport of nanoparticles. In this study, we analyse the absolute density profiles of nanoparticles grown via particle agglomeration under typical plasma enhanced chemical vapour deposition (PECVD) conditions, with or without the inclusion of a thermophoretic force arising from an asymmetrical variation of the electrode temperatures. Particle charging due to the collection of plasma ions and electrons is calculated for every nanoparticle size using the Orbital Motion Limited (OML) theory. The additional complication of stochastic charge fluctuations, which results in an increase in the coagulation rate of the small nanoparticles, is, however, at this stage not taken into account.
In the next section (section 2), the different parts of the developed model are described. First the governing equations for the 1D fluid model are stated, followed by the description of the particle nucleation and coagulation module. The analytical expressions of particle charging and transport are also briefly reviewed. Section 3 discusses the numerical approach. Results are presented and analysed in section 4, whereas the conclusions are given in section 5.

Description of the model
A radio-frequency capacitively coupled plasma in a parallel-plate discharge geometry is described. The plasma chemistry and nucleation are incorporated by a detailed chemical kinetics scheme in the 1D hydrodynamic model that self-consistently determines all plasma properties. A sectional model simulates further nanoparticle growth by coagulation and is iteratively coupled to the fluid model. Transport and charging are considered for each agglomerate size distribution in the nanoparticle domain. An overview of all different submodules is given below, whereas the numerical scheme, describing the coupling between the submodules, is discussed in section 3.

Fluid equations
The fluid model describes the discharge by a combination of continuity and momentum equations for the electrons, positive and negative ions, radicals, molecules and nanoparticles.
For each species, a particle balance equation is constructed 6 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT in which n j represents the particle's density and S j describes the different source and sink terms of species j, including volume reactions, gas inlet and pumping (where the latter is incorporated by an average residence time). The flux terms j of small species (molecules, ions and electrons) are estimated by the drift-diffusion approximation where µ j is the mobility and D j the diffusion coefficient. Ions are too slow to follow the instantaneous electric field E due to their much lower momentum transfer frequency. Therefore, an effective electric field E eff,j is adopted in the drift-diffusion approximation for the ions that accounts for their inertia effects [24]. A separate flux equation is needed to describe the transport of the nanoparticles, as will be explained in subsection 2.4. The electric field E and the potential V are calculated from Poisson's equation which accounts for the number of elementary charges residing on the nanoparticles (Z d ) as well as the number of positive and negative ions and free electrons. n + , n − , n e , n d represent the positive ion, negative ion, electron and nanoparticle densities, respectively. The electron energy balance is also solved in order to obtain the average electron energy density w e = n e dn e dt with the average electron energy. The first term on the right-hand side (RHS) represents the Ohmic heating of electrons, whereas the second term describes the loss of energy due to electron impact collisions, including ionization, dissociation, excitation and electron attachment, as well as the recombination of ions with electrons on the nanoparticle's surface. A separate EEDF module is incorporated to compute the EEDF from the Boltzmann equation in the two term approximation, which is needed to obtain the reaction rate coefficients of every electron-neutral collision, as well as the electron mobility and diffusion coefficients, as a function of the average electron energy, . The electron energy density flux w is given by Other plasma species (ions, neutrals and nanoparticles) are assumed to have the local gas temperature and, hence, no energy balance equation has to be considered.

Particle nucleation
The nucleation module deals with the growth of small gas species to larger molecules that typically consist of hundreds of atoms. The generation of dust in chemically active silane plasmas is an in situ phenomenon that appears to take place by a series of chemical reactions in the ionized gas phase, known as gas phase polymerization [25,26]. Negative ions are generally believed to play a crucial role in the initial stages of particle formation [27]- [29], as they remain 7 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT electrostatically trapped in the plasma due to the action of the sheath potentials, and hence have a longer residence time in the plasma which favours their further growth. Mass spectrometry measurements also validate the role of anions in particle generation, as higher mass anions are observed in the time-resolved mass spectra, whereas neutrals and cations are only detected at lower mass ranges [30]. In our model two different anionic pathways comprising successive reactions between anions and silane molecules start primarily from SiH − 3 and SiH − 2 [14,15], which polymerize into larger silyl (Si n H − 2n+1 ) and silylene (Si n H − 2n ) anions, respectively ) are also incorporated, since these molecules carry sufficient internal energy to overcome the barriers in endothermal reactions [17]. A total of 69 different species (electrons, ions, neutrals, and radicals) and 180 reactions, comprising electron impact reactions with silane and hydrogen molecules, neutralneutral reactions (hydrogen abstraction, SiH 2 insertion, de-excitation of vibrationally excited species), and ion-neutral reactions (anion-silane reactions, and mutual neutralization) is taken into account. More details on the complete chemical kinetic scheme can be found in [14]. Note that besides a detailed description of the gas phase chemistry, the plasma-wall interaction is taken into account by the introduction of a sticking model. In the sticking model, the surface reaction probability of each neutral species is stated and it thus ensures that the loss of species due to deposition is described.
Polymerization reactions in the plasma chemistry module are eventually stopped at anions containing 12 silicon atoms (Si 12 H − x ), since it is not possible to describe the plasma chemistry for unlimited number of plasma species, and these anions then form the starting point of the second stage of particle formation, i.e. coagulation. A direct coupling between the nucleation and coagulation module is made by taking the production rates of all species with more than 12 silicon atoms in the anionic pathways as the source term for the smallest volume section in the coagulation module (see below). As soon as the agglomeration starts, the particle behaviour changes from ionic (determined by gas phase chemistry) to a floating probe behaviour. Therefore, the particles are no longer considered as 'species' but as 'nanoparticles' that can not only grow further by coagulation, but also acquire a charge which depends on the particle size and they are found at a certain region of the discharge where the various forces on the particles are balanced. Before discussing the incorporation of fast coagulation, the electrostatic charging and transport will be briefly reviewed.

Particle charging
Similar to floating Langmuir probes, dust grains suspended in a plasma usually tend to become negatively charged, mainly because the initial electron mobility is substantially larger than the ion mobility. Therefore the electrons collide much more frequently on the nanoparticle's surface, resulting in a negative equilibrium charge. In contrast to ions, the charge on a particle is closely related to its size and greatly depends on the local plasma conditions. Every particle attains a certain floating potential that ensures equality of the electron and positive ion currents 8 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT towards the dust particle surface. For a spherical nanoparticle with radius r d , these currents can be calculated by means of the OML theory when the condition r d λ L applies, where For the collection of Maxwellian electrons and ions, these orbit-limited currents are given by [31] with n i , n e the ion and electron densities, T i , T e the ion and electron temperatures and m i , m e their masses, respectively. k B is the Boltzmann constant and V fl is the floating potential of the particle relative to the average potential of the background gas.
In order to account for the ions drift velocity v i in the plasma sheaths, the mean energy E i is used to replace k B T i in equation (7) and is calculated by [32] with T gas representing the gas temperature. The above equation is obtained by calculating 1 2 m i v 2 s , where v s is the mean ion speed expression of Barnes et al [33], that includes an average thermal component, v th,i = 8k B T gas /πm i , and a drift component, v i , respectively.
Equalizing I e = I i generates the particle's floating potential. The nanoparticle charge Q d is related to the floating potential by with Z d the number of elementary charges and C the capacitance of the particle in the plasma. Note that also the recombination of ions and electrons on the surface of the dust particle is taken into account by means of a recombination rate [32], representing an additional important loss mechanism for plasma ions and electrons.

Particle transport equation
Besides the electrostatic force given by nanoparticle transport is dominated by other important forces, including neutral drag, ion drag and thermophoresis. Gravitational forces can be neglected, since we are dealing here with 9 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT submicrometre particles. Below, the analytical expressions of the forces important to nanoparticle transport will be briefly reviewed. More information can be found in [32,34]. In the present model neutral drag, resulting from collisions with neutral gas molecules, only acts as a damping force, since no advection of the background gas is taken into account, resulting in [35] where m n and n n are the background gas mass and number density, respectively, v th the average thermal velocity of the background gas, and v d the drift velocity of the nanoparticle. Momentum transferred from positive ions accelerated by the electric field causes the creation of an ion drag force that consists of two components, i.e. a collection (F c i ) and orbit force (F o i ), which are given by [33] The collection force accounts for the momentum transferred from collected ions with the mean speed of v s and the collection parameter b c whereas in the orbit force deflected ions transfer their momentum via Coulomb interactions with b π/2 the impact parameter for 90 • deflection. Recent improved theories on the collection force [36] suggest that collisions of ions beyond the Debye sphere need to be taken into account, i.e. the so-called large angle scattering, that results in a revised form of the Coulomb logarithm adopted in equation (15) compared to the original approach of Barnes et al [33] which considered a maximum radius or cutoff radius equal to the linearized Debye length. Usage of the linearized Debye length seems to result in an underestimation of the observed ion drag force by an order of magnitude [37]. In our calculations, the electron Debye length has therefore been used instead of the linearized Debye length for the screening in the Barnes approach since recent simulations have indeed shown that this simple approximation provides for the same effect as the improved theories. The error caused by this simplified approach is at most 25% in a limited part of the discharge [38]. Furthermore, calculations have shown that the Barnes approach with the electron Debye length and the solution of the Khrapak theory give similar results for small particle sizes [39], suggesting that the Khrapak approach should mainly be adopted for the description of the collection force of highly charged dust grains in the micrometre regime.
The additional inclusion of a thermophoretic force arises when a thermal gradient is induced by cooling or heating of the electrodes, that pushes the nanoparticles towards colder regions of the gas discharge. For a gas temperature gradient ∇T gas , the thermophoretic force is found to Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT be [40] F th = − 32 15 where κ T is the translation part of the thermal conductivity and α is the thermal accommodation coefficient, taken equal to 1.
Assuming that the damping neutral drag balances the sum of all other forces, a drift diffusion equation for the nanoparticle's flux can then be derived with µ d = Q d /m d ν md and D d = µ d k B T gas /Q d the nanoparticle's mobility and diffusion coefficient, respectively [34,41]. By adopting this technique the nanoparticles can be treated with the same numerical procedures as the other charged species. Particle transport is included for each section considered in the nanoparticle domain (see below).

Description of the particle coagulation model
In the agglomeration phase, the particles from the nucleation module quickly grow in time from a few nanometres to several tens of nanometres, due to the collision of two smaller nanoparticles resulting in the production of one larger one. Coagulation can generally be treated by using techniques from aerosol physics, where the phenomenon has been studied in fairly great detail. In our code, this second stage of particle generation has been incorporated by making a selfconsistent coupling between our 1D fluid model and an aerosol dynamics model [23,42]. In this model, the temporal evolution of the particle number density in a volume range (υ, υ + dυ), n(υ), can be described by the general dynamic equation (GDE) for an aerosol [42,43] ∂n(υ) ∂t where the first term on the RHS describes the formation of particles in the (υ, υ + dυ) volume range from coagulation of two smaller particles. The factor 1/2 must be introduced as collisions are otherwise counted twice in the integral. The second term accounts for the loss of particles from this volume range due to coagulation with particles of any volume. The formation rate of new particles with volume v 0 by nucleation is described by J 0 (v) in the third term, since δ(v − v 0 ) equals unity when v = v 0 and is zero otherwise. β(u, υ) represents the coagulation frequency between the two interacting particles with a volume u and υ in the free molecular regime (i.e. in the regime where the Knudsen number Kn 1, which is applicable to particles as large as 1 µm in size for the low-pressure plasma range currently under investigation, with Kn being the ratio of the mean free path of the gas molecules to the particle radius) [42] β(u, υ) = 3 4π with ρ d the mass density of the silicon particles, i.e. 2.3 × 10 3 kg m −3 . Note that in the formulation of β(u, υ) resides the implicit assumption that particle charging does not affect the coagulation rates and that nanoparticles do not undergo any charge fluctuations. However, due to the discrete nature of collection of electrons and ions, strong charge fluctuations can especially occur for small nanoparticles, which obey Q d / Q d = 0.5 Q d /e −1/2 [44], and can lead to oppositely charged grains that enhance particle coagulation, as shown in [45,46]. Electrostatic repulsion between large negatively charged grains, on the other hand, has the effect of inhibiting particle coagulation. Methods for correcting the collision frequency function β(u, υ) for particle charging and charge fluctuations will be included in future study. Equation (21) also assumes that the clusters keep a spherical form with solid density in contrast to fractal models that adopt an appropriate fractal dimension corresponding to the degree of irregularity of the particles [43].
Since the GDE (equation (20)) is a nonlinear partial integro-differential equation a sectional approach, schematically illustrated in figure 1, has been adopted to predict the evolution of the particle size distribution within reasonable computation times. In the sectional model, the nanoparticle domain is divided on a logarithmic scale in several volume sections (υ, υ + dυ) each with its own specific average mass, radius and charge. For every section, a general aerosol property q(υ) = aυ b n(υ) is introduced that is assumed to be constant in each section. Depending on the selection of the parameters a and b, it can be used to describe the particle number density n(υ) (a = 1 and b = 0), the volume υ (a = 1 and b = 1) or the surface area distribution (a = π 1/3 6 2/3 and b = 2/3). Here, we have applied a volume-based sectional model, since it has been experimentally shown that the overall particle volume does not change very drastically during particle coagulation [8] and computations have furthermore shown that this model is very good in predicting the particle size distribution for systems in which particles grow by coagulation [47]. By integrating q(υ) over each section the total volume Q k becomes the constant quantity for any section k with V L the lower volume limit and V U the upper volume end of section k. By dividing the entire particle size domain into sections and only dealing with the total volume in each section for which the upper volume limit is a constant factor of the lower volume range, the number of conservation equations required is only equal to the number of sections.
In this computation, we have divided the volume domain logarithmically in 78 sections, comprising particles between 0.80 and 50 nm in diameter. This combination facilitates the solution for particles covering a volume range of five orders of magnitude, in which the initial sections are much smaller than the latter ones. A coupling with the plasma nucleation module is made by taking the production rate of the largest nucleation species, i.e. Si 12 H − x or the so-called stable 'nuclei', as the source term for the smallest volume section. Hence, in order to directly link the coagulation module with the nucleation module, the lower limit of the initial, smallest volume section is taken equal to 0.24 nm 3 , which corresponds to the particle diameter of the stable nuclei Si 12 H − x , i.e. 0.77 nm, when the bulk silicon density of 2.3 × 10 3 kg m −3 is applied. The average volume of the largest volume section is 6.42 × 10 4 nm 3 , corresponding to an average diameter of ∼50 nm. In the code, the coagulation frequencies β(u, υ) need to be determined only once for all possible coagulation combinations. During each coagulation step, the charge and transport of each volume section is determined simultaneously. Hence, nanoparticle charging, transport and coagulation are solved with the same temporal resolution, as will be explained below. Figure 2 shows a schematic representation of the adopted numerical approach and the selfconsistent coupling between the different sub-modules used to describe the dusty silane discharge. All governing fluid, charging, transport and sectional equations are spatially resolved on a uniform mesh comprising 128 grid points in the axial direction, i.e. the direction normal to the parallel plates. A time-stepping procedure has been implemented by adopting two separate calculation cycles, each with its own specific time-step. In the first calculation cycle, i.e. the plasma module, the coupled set of particle balance, momentum and electron energy equations are solved together with the Poisson equation for the self-consistent determination of the plasma properties. The rates for electron collisions k e and the electron transport coefficients µ e and D e are calculated as function of average electron energy in the separate EEDF module by solving the EEDF using the Boltzmann equation in the two-term approximation, and they serve as input values for the plasma module. The first stage of particle generation, i.e. the nucleation, is implemented in this plasma module by a detailed chemical kinetics scheme that starting from SiH 4 leads to the growth of species that contain a maximum of 12 silicon atoms. The time-step within the RF cycle (50 MHz) is set to 2.5 × 10 −10 s (i.e. the RF cycle is divided into 80 time-steps). In order to speed up the calculation, a longer time-step of 10 −5 s is adopted for the neutral-neutral chemistry (not shown). Due to the large difference in timescale between the nanoparticle motion and the RF period, a second calculation cycle with a larger time-step is introduced in the nanoparticle module, which includes the computation of the nanoparticle growth in the coagulation module and the nanoparticle transport and charging module. Hence, for every step in the coagulation module, nanoparticle transport and charging is solved with the same temporal resolution in each of the 78 considered sections. Since nanoparticle motion influences the location of the particle growth, the newly obtained dust density profiles are introduced in the next coagulation cycle.

Numerical approach
An iterative procedure couples both calculation cycles. First the computation of the discharge dynamics is carried out over a number of RF cycles, during which the nanoparticles are assumed to be immobile. Then nanoparticle coagulation, transport and charging are started during the second calculation cycle, using the time-averaged electric field, and electron and positive ion fluxes. Note that the time-averaged negative ion fluxes do not have to be considered, since the negative ions do not contribute to the charging process of the nanoparticles, as the negative ions do not have enough kinetic energy to overcome the negative floating potential of the nanoparticle. Coagulation, transport and charging are run until the relative change in the updated positive ion density distribution is more than 0.01%. Indeed, in order to avoid instabilities in the solution of the Poisson equation and the electron transport, due to the artificially space charges created from adopting time-averaged immobile electron and positive ion densities during the second calculation cycle, the positive ion density distributions are first corrected to the new nanoparticle profiles by moving the total positive ion density with the same amount as the negatively charged nanoparticles before switching back to the next RF cycles. This technique reduces the computational effort and makes sure that interplay of the particle and discharge dynamics is taken into account in a self-consistent way.
The numerical method used to treat the system of nonlinear coupled differential equations is based on an implicit finite-difference technique using the Scharfetter-Gummel exponential scheme for the spatial discretization of the balance equations [48]. Convergence of the fluid model is finally reached when the relative change of the discharge parameters at the beginning of two succeeding RF cycles is less than 10 −6 . More details concerning the applied algorithms and the discretization schemes can be found in [48].

Results and discussion
In this section, the calculated results of the fully coupled 1D model are presented for typical operating conditions in a PECVD reactor that lead to powder formation in silane-based plasmas [24]. The simulated reactor is operated at a pressure of 40 Pa, a radio-frequency of 50 MHz, and a power of 5 W. The inter-electrode gap is set at 3 cm and 20 sccm of pure SiH 4 is fed into the discharge.
In the first subsection (subsection 4.1), the absolute density profiles of the nanoparticles grown via particle coagulation are analysed at a uniform gas temperature of 400 K. The average nanoparticle charge distribution due to the collection of plasma ions and electrons is discussed for several of the 78 considered coagulation sections (i.e. nanoparticle sizes) in subsection 4.2. Finally, the effect of thermophoresis on the location of nanoparticle growth is presented in subsection 4.3, where a thermophoretic force is induced due to the heating of the grounded electrode to a temperature of 500 K, whereas the powered electrode remains at 300 K. Figure 3 shows the spatial profiles of the nanoparticle number densities grown via particle coagulation at the discharge conditions mentioned above. The results are depicted for a quasisteady-state particle size distribution, where the growth of a particular particle size is compensated by their losses as a result of coagulation and transport. This quasi-steady state is also established in other aerosol studies. The size indicated on each plot corresponds to the average particle diameter of the respective coagulation section considered in the particle coagulation module. Several specific nanoparticle profiles are depicted in figures 3(a)-(h).

Density profiles of the nanoparticle agglomerates
The figure clearly illustrates that the density profiles of the charged nanoparticles are influenced by the balance of several forces acting on the particles. In these computations no thermophoretic force is applied, since we assume a uniform gas temperature of 400 K over the entire reactor. Therefore, the nanoparticle transport will be dominated by the concerted action of the electrostatic force and the ion drag force as we assumed that the neutral drag force only acts as a damping force on the nanoparticles. Indeed, charged particles are affected by the electric field in the plasma and in the sheath regions in two separate ways. The sheath fields first repel the negatively charged particles from the plasma sheaths and effectively trap the nanoparticles in the bulk of the discharge. On the other hand, the electric field accelerates positively charged particles out of the discharge leading to a momentum transfer from positive ions in the vicinity of the dust grains that creates an ion drag force which will drive nanoparticles in the direction of the net positive ion flux, i.e. towards the plasma boundaries. The competition between these two different forces determines the nanoparticle location in the discharge reactor and thus also the position of nanoparticle growth due to coagulation. Relatively high concentrations of small nanoparticles ( 1 nm diameter) are observed in the bulk of the discharge (figure 3). When going to larger nanoparticle sizes (figures 3(b) and (c)), their specific spatial distribution becomes narrower, as the electric field forces higher negatively charged particles more to centre of the plasma. Particles with diameters larger than 27 nm (figures 3(f)-(h)), however, begin to experience the influence of the ion drag force that pushes the particles towards the wall surfaces at both sides of the reactor vessel and decreases the particle density at the plasma centreline. Eventually for particles exceeding the 30 nm range two peaks of equal magnitude on both sides of the sheath boundaries are formed, where the electrostatic force and the ion drag force are balanced (cf [34]).
Coagulation also causes a decrease in concentration as particles grow to larger sizes. Similar to experimental observations the calculated density of the clusters seems to dramatically drop by several orders of magnitude. This is clearly shown in figure 4, where the volume-integrated densities of the various volume sections are depicted.

The particle charge distribution
Particle charging due to the collection of plasma ions and electrons is calculated for every volume section considered in the coagulation module using the OML probe theory. Due to the high mobility of electrons compared to ions, the nanoparticles generally become negatively charged. The computed charge for several particle sizes is shown in figure 5. It is found that the number of elementary charges collected on the nanoparticle's surface is not only dependent on the particle size but is also very sensitive to the location in the discharge, as the charge of a specific nanoparticle varies significantly through the discharge according to the local plasma conditions, i.e. the ion and electron densities and the electron temperature. In the centre of the discharge, the weak electric field impedes electrons to overcome the repulsive Coulomb force arising from the dust particle's similar polarity. This results in a minimum in the grain charge at the plasma bulk of the discharge. In the sheath regions, where large electric fields are present, electrons gain energy to overcome the repulsive Coulomb barrier, causing an increased electron attachment and thus a maximum in the nanoparticle's charge. However, close to the electrodes practically no electrons are present, resulting in a sharp decrease of the particle charge. Furthermore, the average electron charge on a particle is proportional to the particle diameter as large nanoparticles are charged more negatively than small nanoparticles. This is also shown in figure 6, where the number of elementary charges is depicted as a function of nanoparticle size at two different plasma locations, i.e. in the plasma centre, where a minimum charge is attained, and at the maximum charge near the sheath edges. Hence, due to the larger negative charge near the sheath edges, an electrostatic barrier is formed that effectively traps negatively charged particles in the bulk of the discharge, where subsequent growth by particle coagulation can occur.
From figure 5, we can also deduce that small nanoparticles have a greater probability of becoming neutral or even positively charged, as the charge on particles below 5 nm is less than three elementary charges (see figure 5(b)). Therefore, these particles can undergo strong charge fluctuations due to the statistical variations in electron and ion fluxes towards the particle's surface. These charge fluctuations can not only enhance particle coagulation between oppositely charged grains, but also provides a means for nanoparticles to escape from the discharge when they become neutral or positively charged, by diffusion (for neutrals) or drift (positively charged particles). The larger negatively charged nanoparticles, on the other hand, in figure 5(a), become trapped in the discharge by the sheath electric fields and grow to larger sizes. The stochastic fluctuations on the average nanoparticle charge will be added to the model in the near future and implemented in a more sophisticated coagulation scheme, where particles of opposite charges will then be able to coagulate very fast, whereas particles of the same charge will be repelled and do not collide.

Influence of a thermal gas temperature gradient
The force balance on coagulating particles can seriously be affected by the additional incorporation of a thermophoretic force arising from a thermal gradient in the gas temperature [41,49]. Hence, we also performed calculations where an additional thermophoretic force is induced by heating of the grounded electrode (at 0 cm) to a temperature of 500 K, whereas the powered electrode (at 3 cm) remains at 300 K. The thermophoretic force is added in the transport equation used to calculate the motion of nanoparticles grown via particle coagulation. The effect of thermophoresis on several of the calculated particle concentration profiles is shown in figure 7 for similar operating conditions (pressure, power, RF). At an electrode temperature difference of 200 K it is found that the particle concentrations displace towards cooler regions in the direction of the powered electrode, and are thus influenced by the thermophoretic force.
For small particles (∼2 nm size) a minimal effect on the particle density profile is observed (see figure 7(a)) as the nanoparticle growth still mainly occurs in the centre of the discharge. Their transport is thus still largely dominated by the electric field that exerts an electrostatic force in the direction opposite to the field, and thus away from the plasma sheaths. The growth towards larger particles, however, gradually starts to shift in the direction of the colder, powered electrode due to the action of the thermophoretic force. Since the thermophoretic force scales more strongly with the radius of the nanoparticle than the electrostatic force, the densities of particles larger than 3 nm in diameter start to noticeably diminish in the centre of the discharge as the electrostatic and the thermophoretic force balance further towards the sheath regions of the powered electrode. As the particle size keeps increasing the thermophoretic force is eventually able to drive all particles away from the middle of the plasma and a narrow peak in front of the powered electrode is formed for 30 nm particles. At this location, the thermophoretic and ion drag force are balanced by the electrostatic force.
Hence, we can conclude that thermophoresis is found to be a useful means for the control of particle contamination exceeding the 2 nm size. High electrode temperatures can indeed sweep particles away from sensitive areas of the discharge before they reach the processing substrate.

Conclusions
Formation and subsequent growth of nanoparticles in a low-pressure capacitively coupled silane discharge has been investigated with a self-consistent determination of the plasma properties by combining an aerosol dynamics model with a self-consistent 1D fluid model, including nanoparticle nucleation, transport and charging. This model is used to predict the number density, the charge and the transport of nanoparticle agglomerates undergoing particle coagulation. In the coagulation module, a sectional approach is adopted to solve the GDE for an aerosol by dividing the nanoparticle domain in 78 volume sections that covers nanoparticle growth over two orders of magnitude, i.e. from ∼ 0.8 to 50 nm in diameter. For every step in the coagulation module, Calculated spatial distribution of several nanoparticle densities for an asymmetric variation of the electrode temperatures. The grounded electrode (at 0 cm) is heated to 500 K, whereas the powered electrode (at 3 cm) remains at 300 K. The number density profiles of 15 and 20.5 nm particles need to be multiplied by a factor 2 in order to obtain the values of their absolute densities.
the effect of particle charging as well as external forces such as the ion drag force, electrostatic force and thermophoresis is included and solved with the same temporal resolution. The absolute densities of nanoparticles grown via particle coagulation are presented for typical conditions of PECVD process variables. As particles grow to larger sizes, coagulation will cause a decrease in the particle concentration. It is found that the location of particle growth depends on the nanoparticle size. Small negatively charged nanoparticles (below 27 nm) are excluded from the plasma sheaths by the repulsive electrostatic force and show peaks at the centre of the discharge, whereas most larger nanoparticles are located around the sheath boundaries owing to the balance of the electrostatic force from the sheath region and the ion drag force arising from the acceleration of positive ions towards the electrodes.
The particle charge, calculated via the OML theory, is very sensitive to changes in particle diameter and local plasma conditions. A minimum in the charge is attained in the bulk plasma region. An increased electron attachment is observed in the sheath regions, resulting in a maximum in the nanoparticle's charge that forms an electrostatic barrier which effectively traps the negatively charged particles.
Finally the effect of thermophoresis arising from a thermal gradient in gas temperature, induced by a 200 K temperature difference between the two electrodes, on the spatial distribution of coagulating nanoparticles has been evaluated. It is shown that thermophoresis can have a strong effect on the location of particle growth and can be used as a physical method to push particles away from the processing area.