Relating the rate of growth of metal nanoparticles to cluster size distribution in electroless deposition

Abstract

Electroless deposition on patterned silicon substrates enables the formation of metal nanomaterials with tight control over their size and shape. In the technique, metal ions are transported by diffusion from a solution to the active sites of an autocatalytic substrate where they are reduced as metals upon contact. Here, using diffusion limited aggregation models and numerical simulations, we derived relationships that correlate the cluster size distribution to the total mass of deposited particles. We found that the ratio ξ between the rates of growth of two different metals depends on the ratio γ between the rates of growth of clusters formed by those metals through the linearity law ξ = 14(γ − 1). We then validated the model using experiments. Different from other methods, the model derives k using as input the geometry of metal nanoparticle clusters, decoded by SEM or AFM images of samples, and a known reference.

---

1 Introduction

Electroless deposition is a chemical method of plating in which a metal is deposited on a surface without the support of external fields or driving voltages.^1–5 In the technique, absorption of metal ions at the surface is activated by catalytic and reducing agents added in the solution with the ions. If the substrate is made up of silicon, the substrate itself performs simultaneously as a reducing and catalytic agent and the process of deposition is greatly facilitated. Since silicon serves as a basis of the majority of micro- and nano-fabrication processes and techniques and can be easily integrated with other materials,^6–9 electroless deposition on silicon displays potential to be used for the synthesis of metal nanoparticles, metal nanomaterials or other nanostructures.

The method enables the fabrication of nanostructures by incrementally depositing smaller building blocks on a flat substrate. Pre-patterning of the substrate by optical or electron beam lithography allows site selective formation of metal nanoparticles with tight control over the final aspects of the nanoparticles.⁵ On the active sites of the silicon surface exposed to growth, metal ions aggregate into clusters where the geometrical characteristics of the aggregate will depend on a fine-tuning of size and shape of preexisting patterns on the substrate, and the parameters of electroless growth including temperature, pH, and concentration of metal ions in solution.⁵

In the absence of stirring or other convective flows, ions are transported towards the active sites of the substrate by pure diffusion. Since the mechanisms of diffusion are very well understood,^10,11 electroless deposition can be simulated using discrete models of transport in cellular spaces. In the models, cells of a grid can take definite 0 or 1 states – 1 indicates the presence of a particle (ion). A particle will move with time in a continuous time probabilistic Brownian motion (Fig. 1a), which is discretized in the space as a random walk (Fig. 1b). Particles aggregate upon contact. Displacement of a great many particles in a lattice (Fig. 1c) enables the formation of numerical aggregates with a dendritic appearance that is typical of fractals (Fig. 1d). Simulation of particle growth using iterative arrays is called Diffusion Limited Aggregation (DLA).^12–17 In previously reported analyses, we have used DLA models to examine the effects of pattern size^4,5 and pattern distance¹⁸ on the characteristics of metal nanoparticles in electroless growth. Here, we revise DLA models of particle growth to consider the effects of the sticking probability of an ion to the aggregate – p.

Fig. 1 Diffusion of metal ions in a solution is described by a probabilistic continuous time Brownian motion, in which ions change arbitrarily the direction of velocity upon collision with atoms in the solution (a). The process can be represented in a discrete grid (b). Multiple ions can be placed in a domain to reproduce electroless deposition (c). Each ion obeys the laws of diffusion following a random walk – upon contact with a line of nucleation sites, ions are incorporated into an aggregate of particles growing with time (d). The structure of occupied sites of the aggregate has a dendritic appearance typical of fractals. The density of aggregates can be modulated by changing the probability of adhesion of a particle to the aggregate p (e and f).

On changing p between 0 and 1, one can modulate the packing fraction and superficial aspect of the aggregates. Low values of p generate aggregates with high atomic packing densities (Fig. 1e), in contrast to low-density aggregates obtained from high values of p (Fig. 1f). p regulates the adhesion of ions to the substrate. If p is low (p → 0) the probability of adhesion of a particle to the aggregate is low, meaning that a particle may not immediately stick to the aggregate upon contact. Under these conditions, a particle may penetrate more deeply in the structure without being captured by the external shells of the aggregate. Higher penetration depths imply, in turn, higher final densities of the numerical clusters of particles. Notice though that, even if for small p the surface density of the aggregates is high, the growth of the aggregate is slow compared to systems with larger sticking probabilities, as demonstrated in the following sections of the paper. If p = 0, no adhesion can occur on the surface.

Using data from the simulations, we developed models that correlate the cluster size distributions and the mean cluster size with p. Since p is related, in turn, to the kinetics of metal deposition, these models enable extraction of the rate of growth of nanoparticles from topographical maps of samples, obtained by SEM, AFM or other similar techniques of imaging. We tested the performance of the model analyzing the growth of gold and silver aggregates in patterned silicon substrates. The model predicted a difference of growth rate between the two metals of a factor of 10. The growth kinetics of metal nanoparticles measured by UV spectrophotometric techniques matched the predictions of the models with a good level of accuracy.

Coupled to experimental SEM or AFM data, this scheme can be used to derive the growth rate kinetics of metals plated on a substrate by electroless methods.

2 Results

2.1 Generating numerical DLA aggregates

Using the numerical scheme described in the methods and in Fig. 1, we generated numerical aggregates of metal particles on patterned silicon substrates. Particles in the aggregate form non-continuous dendritic structures with details over multiple scales that are typical of fractals. For each imposed value of sticking probability p, we extracted the local density at different distances from the substrate r. Information on the density of the aggregate at a local scale was used to derive the density–density correlation function c(r).

The density–density correlation function c(r), also called pair correlation function, is a distribution routinely used in statistical mechanics to describe how density varies as a function of distance from a reference point. In the present configuration, we used as a reference point the lines of nucleation sites of the aggregate, i.e. the base of the aggregate. c(r) is a measure of the probability of finding a particle at a distance r from the reference. It describes the internal structure of an aggregate as a function of continuous, smoothly varying spatial coordinates – for this it can be used to link the microscopic to the macroscopic characteristics of a system: c(r) delivers the information content of an aggregate as a function of the scale of the aggregate.

In a log–log plot, c(r) is a line with a slope α (Fig. 2). Since the fractal dimension D_f of an aggregate is an index that quantifies the change in detail to a change of scale, D_f may be derived from α as D_f = 2 − α (Methods). The fractal dimension varies with the sticking probability and rapidly undergoes transition from an initial value D_f ∼ 1.9 for p = 0 to D_f ∼ 1.625 for p = 1. Notice that for sufficiently high values of p, p > 0.2, the best fit of D_f attains the steady state value D_f ∼ 1.625, similar to the theoretical limit 5/3 ∼ 1.667.¹⁴ The purpose of calculating the fractal dimension is deriving the cluster size distribution of an aggregate. Clusters are defined as naturally separated trees – a collection of particles connected to the same nucleation site through nearest neighbours forming the deposit^16,17 (Fig. 2d). From ref. 16, the average cluster size 〈S〉 of an aggregate with fractal dimension D_f and mass N is

〈S〉 ∼ N^{D_f/(D_f−1)} (1)

the total number of clusters or branches Q in an aggregate is therefore

Q = N/〈S〉 (2)

Fig. 2 Aggregates can be analyzed to derive the local density ρ(r) of the aggregate at different lengths r from the line of nucleation sites (a). The local density of the aggregate is used in turn to derive a density–density correlation function C(r) following a power law dependency. In a log–log plot, C(r) has the appearance of a line (b). The slope of C(r) is used to derive the fractal dimension D_f of the aggregates as a function of the probability of adhesion p, D_f attains a steady state value D_f ∼ 1.625 for p > 0.2 (c). The fractal dimension is correlated with the number of clusters in a numerical aggregate (d).

The cluster-size distribution I(S,N), i.e. the number of clusters in an N-aggregate as a function of cluster size S, is:¹⁹

The probability of finding an S-site tree on a nucleation site is limited by N.

2.2 Structure of aggregates and cluster size distribution

The structure of occupied sites was characterized applying relationships (1) to (3). The mean cluster size was estimated by arithmetic average after calculating the largest branch and fixing the smallest cluster containing 10 particles. Fig. 3a reports the mass N of the aggregates as a function of the time of deposition (i.e. the number of iterations) for different values of the sticking probability p. For a fixed p, the size of the aggregate increases with a quadratic function of time, which is easily explained considering that the larger the number of pixels in an aggregate, the larger the number of docking sites to which diffusing ions in solution can possibly bind. For a fixed time, the size of the aggregate increases linearly with the sticking probability p. Thus, approximately i = 6 M iterations from the beginning of the simulations, the number of particles in the aggregate is N ∼ 500 for p = 0.1, rises to N ∼ 650 for p = 0.6 and reaches N ∼ 800 for p = 1 (Fig. 3b). Remarkably, the linear form of the N(p) characteristic indicates that the kinetics of metal deposition k is proportional to the sticking probability pwhich in turn implies that p can be used to compare metal growth rates across different simulations. Fig. 3c describes the cluster size distribution of the aggregate for different values of the sticking probability p, derived from the numerical analysis of the aggregate and eqn (3). For a fixed sticking probability, the number of clusters I(S) with a certain size S decreases hyperbolically with S. Thus the probability of finding clusters with a number of pixels that is lower or equal to 10 is above the mean of the distribution and decreases below the mean of the distribution for any S > 10. The number of clusters with a size greater than 80 pixels is nearly 0 (p = 1, i = 20 M, N = 5000). This trend is maintained for values of p lower than 1, for which the overall number of clusters and the average cluster size are shifted to the left of the diagram. From the distribution of clusters in the aggregate and eqn (1) to (3), we derived the number Q of clusters as a function of time for different values of p (Fig. 3d). We observe that Q varies linearly with time and that the slope of Q can be modulated by changing p – Q is indicative of the number of separate clusters, i.e. macrograins or isolated particles that are deposited on a substrate during an electroless process. Fig. 3d suggests that the number of isolated particles resulting from an electroless deposition, easily measured using SEM or AFM techniques, can be correlated with p that in turn is associated with the rate of metal growth k through eqn (4).

Fig. 3 Numerical aggregates were examined to extract the total number N of particles in an aggregate as a function of time for different values of the probability of adhesion p (a). For a fixed time i = 0.5 M iterations, N varies linearly with p (b). Numerical aggregates are formed by clusters with different sizes S and numbers I, where I varies hyperbolically with S – cluster size distribution I(S) is determined as a function of the probability of adhesion p (c). The total number of clusters Q in an aggregate is a linear function of time – iteration number – i and probability p (d).

2.3 Connecting the rate of metal growth to the kinetics of cluster deposition

Fig. 3d, reproduced in a separate Fig. 4a for the sake of clarity, reunites in one unique diagram (i) the time evolution of the number of clusters in an aggregate Q̇ and (ii) the sticking probability p. The sticking probability is proportional to the kinetics of metal deposition p → k through relationship (4).

Fig. 4 The cluster size distribution of numerical aggregates was used to determine the ratio between the rates of growth of two materials ξ as a function of the ratio between the rates of growth of the clusters γ formed by those materials over time, the form of the ξ(γ) relationship is linear (a). The model enables estimation of the rate of metal growth k from SEM or AFM images of a sample and a known reference (b).

Q̇ is a measure of the number of separate particles that are deposited on a substrate over time. This assertion is substantiated by the following argument. A cluster is a structure, in the aggregate, with some internal correlation. An operational definition of a cluster is that it is a structure clearly distinguishable as a subsystem – a self-contained system within the larger system. For this, it is the analogue of an isolated particle in a real process of chemical deposition. Therefore, can be considered a good estimate of the number of separate particles that are deposited on a substrate in the unit time.

k is a key parameter in materials science and nanotechnology and is relevant for the rational design of processes or structures that imply the deposition, growth and self-assembly of metal nano-materials. On the other hand, Q̇ is experimentally observable and can be readily determined through conventional imaging techniques, including SEM and AFM. It may therefore be convenient to express p as a function of Q̇.

Since numerical simulations and Fig. 4a reflect metal deposition less quantitatively than qualitatively – calibration standards may be required for matching the predictions of the model with experimental data (see comments in the Discussion) – in the following we will use ratios between variables rather than the absolute values of variables to describe particle growth.

We then define the non-dimensional parameters ξ = p₂/p₁, γ = Q̇(p₂)/Q̇(p₁), the subscript i = 1, 2 indicates two different states or materials with p₂ > p₁. Then, we derive from the diagram in Fig. 4aξ as a function of γ: ξ(γ). The best fit of data and graphical representation of ξ(γ) (Fig. 4b) indicate that ξ varies linearly with γ as

ξ = 14.34γ − 13.7 (5)

we calculate a r-squared statistic r² to test whether the data in several different bands are consistent with the matching template. Values of r² near unity and of estimated variance σ near zero indicate that the signal is consistent with the model (r² = 0.99, σ = 0.092).

In separate ESI 1,† we provide a full statistical report of the linear regression, including values for the parameter table, r², adjusted r², estimated variance and ANOVA table. We used at least 25 simulations for each p to fit data and evaluate the model.

Eqn (5) puts in relation the growth kinetics of different materials expressed in different forms. ξ is the ratio between the kinetics of growth of two materials, γ is the ratio between the rates of growth of clusters emerging from the deposited volume of those materials. The proportion between ξ and γ is ∼14. Say that we have two different metals A and B, the kinetics of deposition of A is known, being k_A. Then k_B can be readily derived using an approximate form of eqn (5):

where the velocity of cluster/particle deposition Q̇ may be determined through direct optical inspection of samples.

2.4 Predicting the relative rates of growth of silver and gold

We used experimental data to benchmark the model. Using electroless deposition techniques described in the methods, we deposited clusters of (i) silver and (ii) gold nanoparticles on patterned silicon substrates. Patterns on the silicon surface are hexagonal arrays of disks penetrating into the positive resist S1813 with a size d, spacing δ and thickness h: d = 10 μm, δ = 30 μm and h ∼ 1 μm.

On comparing the experiments with the predictions of the model, we propose that the aspect ratio of the patterns exposed to growth is the same in the experiments and in the numerical DLA scheme. In this configuration, the patterns are short systems, with d > h and aspect ratio lower than one . Since the model is a comparative method of analysis, it does not necessitate a one-to-one correspondence between the geometry of the real physical prototype and the model. More sophisticated evolutions of the model that will be developed over time will enable direct simulation of real systems and the determination of absolute values of growth. As regarding spacing between patterns, δ, we have used in the experiments spacing between patterns three times larger than the pattern size, which guarantees non-interference between patterns¹⁸ and justifies the use of a single-well numerical scheme and an isolated system to simulate particle formation.

Then, we examined the structure of metal clusters at different times t from immersion in the electroless solution, t = 5, 20, 60, 120 s, using scanning electron microscopy (SEM) and atomic force microscopy (AFM) imaging techniques. We fabricated at least 5 different samples per time of deposition, and acquired more than 25 SEM images and 2 AFM images for the samples. We used SEM images for determining the rate of growth of the nanoparticles, and AFM imaging to extract the topographical details of the clusters.

For fixed values of time, we observe that the number and density of isolated particles on the silicon surface are greater for silver compared to gold deposition (Fig. 5a). This trend is maintained for all considered time frames of growth (Fig. 5b). Using standard image analysis algorithms, we extracted the number Q of isolated particles formed during silver and gold deposition over time (Fig. 5c). Q varies linearly with time (r² = 0.989 for gold, r² = 0.984 for silver; the statistical significance of the linear regression is reported in separate ESI 2†), in agreement with the theoretical model and the simulations. We elaborated information contained in Fig. 5c to derive how rapidly the number of clusters in an image changes with time: at any time deposition of silver is more rapid than that observed for gold. Numerical analysis of data yields the values for silver and for gold, and a ratio between the two γ^Ag/Au ∼ 3.45. This value of γ will be used in the model to determine k.

Fig. 5 We benchmarked the model using experiments. Using electroless deposition on patterned silicon substrates, we obtained silver and gold nanoparticle clusters, with different topographies, cluster sizes and cluster size distribution between metals for a fixed time of analysis – the bar length in the image is 100 nm (a). We derived the time evolution of the number of clusters Q in a pattern over time (d) for either silver (b) and gold (c) from SEM images of samples taken at different time intervals t = 5, 20, 60, 120 s – the bar length in the SEM images is 5 microns.

From AFM imaging (Fig. 6a) one can observe that the number of particles deposited during electroless growth of silver is larger than the number of gold nanoparticles deposited in the same time (60 s), in accordance with the results of the model and SEM inspection of samples. Few clusters of gold nanoparticles are particularly large, this is deceptive and may suggest that the cluster growth rate of gold is larger than silver – that it is not as proved by image analysis algorithms applied to AFM images that yield an estimate of Q_Ag ∼ 106 clusters in a square pattern of 2 μm for silver, and Q_Au ∼ 29 clusters for gold, with a ratio . While AFM imaging may be as accurate as SEM in deriving the topographical characteristics of a sample surface, we use here SEM results to benchmark the model because SEM imaging of samples is faster compared to AFM. SEM images are more in number, have larger formats, include more particles, are more informative and statistically significant than AFM images. While AFM can achieve ultra-high resolution and can examine sample topography at the sub nanometer level, it is not high-throughput.

Fig. 6 AFM imaging was used to investigate the topography of the silver (a) and gold (b) clusters of nanoparticles. Fourier analysis of AFM data enabled derivation of the power spectrum density function and the fractal dimension of the samples as D_f ∼ 2.65 (c and d).

Post-processing of AFM topographic data enabled us to derive the power spectrum density function for both silver and gold (Fig. 6b). From this, we derived the values of fractal dimensions for the clusters of Ag (D_f = 2.7) and Au (D_f = 2.65). Notice that, remarkably, the experimental values of fractal dimension derived in the space are about one dimension higher than the corresponding numerical values derived in the plane, D^3d_f ∼ D^2d_f + 1.

Using the ratio between cluster growth dynamics (γ) determined from experiments in the model of eqn (6), we estimate that the ratio between the rates of growth of silver (k_Ag) and gold (k_Au) is about ξ(γ) = ξ^model ∼ 34.

2.5 Model validation

Using data from electroless deposited silver and gold nanoparticle clusters, we found that the ratio between the kinetics of growth of these two metals is ξ ∼ 34. We then verified the output of the model using independent UV spectrophotometric measurements (Methods). We produced nanoporous silicon particles (NP-Si-NPs) using electrochemical porosification of bulk silicon^20,21 followed by ultrasonic disintegration of silicon into nanoparticles.²² We then suspended NP-Si-NPs into two separate aqueous solutions with (i) hydrofluoric acid (HF) and silver nitrate (AgNO₃) and (ii) HF and gold chloride (AuCl₃) (Methods). Immersion of NP-Si-NPs in solution triggered immediate reaction of deposition of gold and silver within the porous silicon matrix of the particles. We used UV/Vis spectrophotometry to analyze the resulting electroless reaction of deposition. We determined the dynamics of Au and Ag nanoparticle formation as the onset of the variation of absorbance of the solutions measured at λ₁ = 560 nm for gold (Fig. 7a) and λ₂ = 460 nm for silver (Fig. 7b). Rate constants of growth were determined as the reciprocal of the time constant τ resulting from the best fit of data with a first-order system model (Fig. 7c and d) – the time constant is defined as the time necessary for the system to reach 90% of its final value. From the fit, we determined k_Ag ∼ 0.02 ± 0.0034 (s⁻¹) (τ ∼ 3.4 ± 0.57 s) and k_Au ∼ 0.00054 ± 0.00007 (s⁻¹) (τ ∼ 120 ± 14 s), with an experimentally derived value of . Data are reported as mean ± standard deviation. Analysis was carried out over a set of 10 measurements, another example of which is reported in separate ESI 3.† The measured ratio of growth kinetics ξ^exp matches the predictions of the model ξ^model with a good level of significance and ∼93% accuracy.

Fig. 7 Independent techniques of analysis were used to derive the relative values of growth rates of gold and silver nanoparticles within porous matrices in an electroless solution. UV/Vis spectra of samples were analyzed. Growth dynamics of Au (a) and Ag (b) were determined as the variation of absorbance at wavelengths λ₁ = 560 nm for gold and λ₂ = 460 nm for silver, over time. Time evolution of absorbance values, fitted using the output of a first order dynamic system (c and d), enabled us to derive the relative ratio of metal growth rates as .

3 Discussion

The rate of growth k in electroless deposition is an index that quantifies how rapidly a metal is deposited on a surface. If ν is a certain measure of the bulk metal being deposited, then k is the change of ν to a change of time, k = dν/dt. ν can be the mass, volume, or thickness (extensive properties of a system), or concentration or density (intensive properties of a system), thus k is usually given in terms of the magnitude of ν per s. In all cases, measuring k implies estimating the magnitude or intensity of some physical quantity as a function of time. Consolidated methods for the determination of k are cyclic (CV) or linear sweep (LSV) voltammetry, or methods concerning the use of a quartz crystal microbalance (QCM). The first CV or LSV methods measure the shift of peak potentials ΔE_p that develop under conditions where an external voltage of increasing intensity is cyclically (CV) or linearly (LSV) applied to the system.^23–25 Then, ΔE_p is correlated with k. The second QCM method measures a mass variation in the unit time by measuring the change in frequency of a quartz crystal resonator.^26–28 While they are accurate in determining the rate constants of reactions involving reduction of metals, these methods (a) require the acquisition, installation and set-up of electrochemical workstations and electrodes in sometimes complicated configurations (CV and LSV methods) – or of QCM systems – that are not always available in nanotechnology or materials science laboratories; (b) may be conducted under conditions that deviate to a significant extent from the real conditions of metal deposition; (c) raise the problem of the integration of different components or parts: using hydrofluoric acid in a QCM or connecting a patterned silicon substrate in a three-electrode configuration in CV is not obvious; (d) necessitate repeated measures in wide time or frequency intervals for data fitting and the identification of the parameters (i.e. k) from those fits (definition of a calibration standard); (e) realize blind estimates of a physical quantity (i.e. potential, current, or mass) of a system that do not necessarily account for the geometrical variables of that system, i.e. size, shape, distribution, topographical or topological properties of metal nanoparticle clusters that are being formed over time.

The last point implies that, while in nanotechnology the geometrical characteristics of metal nanomaterials are relevant for their functions,^29–33 traditional methods for the determination of k ignore the geometrical form of the material assembling into nanostructures.

Here, we used mathematical modelling and numerical simulations to correlate the topological properties (Q̇) of a system grown by electroless deposition with the variation of its volume (Ṅ → p) to determine k. Different from other existing methods, our approach centers on geometry. The geometry of target structures, in turn, can be easily reconstructed using SEM, AFM, or other techniques of imaging that can be easily found in a nanotechnology laboratory and that are, on the other hand, routinely used to examine the aspects of the structures and the efficiency of a method of fabrication. Thus, incidentally, much of the SEM, AFM, TEM material produced during ordinary inspection of samples and reports can be redirected as the input of the model – generating high volumes of data (big data) that can contribute to reducing noise, reducing uncertainty, and increasing the accuracy of the estimate of k. Moreover

(i) Our model uses as input the ratio between cluster multiplication velocity of two species, . In the more general case, Q̇_A and Q̇_B should be determined as the variation of Q over time, which in turn implies determining the number of clusters Q at different times through multiple measurements. Nevertheless, assuming linearity, γ can be more simply determined as the ratio between one of the possible couples Q_A(t_a) and Q_B(t_a), , measured at any time t_a comprised between zero and the final time of growth, 0 < t_a < t_fin. This reduces the search for k to the estimate of the number of clusters (i.e. isolated metal nanoparticles) in two different SEM or AFM images, and the division between the two – that dramatically reduces the number of measurements and the time necessary to find k compared to other experimental techniques of analysis. In determining the number of clusters Q from SEM or AFM imaging of samples, it is preferable that the nanoparticles do not overlap (or overlap partially) on the substrate to assure accuracy. Images may be analyzed using algorithms developed over time,^5,34,35 which deconvolute information using maximum likelihood methods and Fourier transform decomposition/reconstruction of images. Image analysis is accurate under sparse and sub-confluent conditions. When particles on the substrate overgrow, non-dense assumption breaks down and the number of clusters/particles in the field of view can be under-estimated. To avoid miscalculations, one should limit the analysis to the early time of particle deposition, or in any case to the linear regime of particle growth, which is the case of Fig. 5d and results presented in this work. When Q(t) deviates from linearity, results of the analysis can be inaccurate.

(ii) Our model derives k using direct optical, SEM or AFM inspection of samples. The method does not require sample treatment, preparation or modification, it is not destructive, does not hamper and has no adverse loading effects in the process of growth. The measure of k is carried on the real physical prototype and not on a simplified version of it. The size and shape of the patterns on the silicon exposed to growth are preserved and their effects are correctly incorporated in k, which must not be therefore further adjusted.

(iii) Our model is general in nature. It simulates the assembly of ions or atoms into supramolecular structures and can be used to estimate the ratio between the rates of growth of several materials, including non-metallic materials, provided that the process of deposition is limited by diffusion and the process of chemical reaction at the interface is fast compared to diffusion. Under similar assumptions, a material can be examined using DLA. As a rule of thumb, if a material exhibits a dendritic structure, with some level of order and recursive patterns in the structure – like crystals – it can be analyzed using this model.

In its current form, the model makes the following assumptions and has the following limitations:

(1) In electroless deposition, motion of metal ions in a solution is much slower than the chemical reaction of reduction of those ions into metal on a patterned silicon substrate (DLA assumption).

(2) The diffusion lengths of different metals is the same – assuming that the thermodynamics conditions of growth are held constant, using the Stokes–Einstein relationship¹⁰ this is equivalent to state that the sizes of different metals is the same.

(3) The differences between different species are lumped in the sole chemical reaction of deposition, i.e. in p.

(4) It is a comparative method of analysis. This depends on the fact that, in this current form, the model shows a very high sensitivity to the geometry of the system. Since it may be impractical and computationally intractable to reproduce the entire physical system and simulate a real electroless process of growth, we use ratios between variables rather than the absolute values of variables to describe particle growth. In doing so, we can examine exclusively and focus on the parameter of interest, i.e. k, with the remaining conditions held constant. The effects of an inaccurate representation of the systems would cancel each other out when we compare two different metals. This in turn reduces uncertainty and increases the precision of the model. The cost of an augmented precision is the use of the model as a comparative rather than absolute method of analysis.

More sophisticated formulations of the model that will be developed over time may potentially relax the above constraints and be used for a more complete and precise description of metal growth for applications in materials science.

4 Conclusions

We used DLA simulations to develop models that correlate the ratio (γ) between the rates of particle growth of different metals with the ratio (ξ) between the rates of deposition of those metals in the bulk form. The best fit of numerical simulations yields a linear relationship between γ and ξ that reads ξ = 14(γ − 1). Particle growth kinetics is more than a factor of 10 lower than material growth kinetics. We validated the model using electroless deposition of gold and silver nanoparticles on patterned silicon substrates. Different from other existing methods, including CV, LSV or QCM, the method enables determination of the rates of growth (k) from the geometry of clusters of metal nanoparticles and a known reference. SEM or AFM micrographs, routinely taken in a laboratory during sample inspection and analysis, can be used as a source for the model, and the yield and estimation of k even without dedicated equipment.

5 Materials and methods

5.1 Simulating metal nanoparticle growth

We used a diffusion limited aggregation (DLA) model to simulate metal nanoparticle growth. In DLA, randomly displaced particles stick together to form the intended structures.^{4,5,12,14–16} The model assumes that the main mechanism of metal ion transportation in a solution is Brownian diffusion. This assumption implies that the metal reduction and deposition on the silicon surface is instantaneous, which is true when diffusion times are much larger than the time associated with surface chemical reaction on the silicon surface.^4,5,18 The assumption also implies that convective flows are vanishingly small everywhere in the domain. Brownian motion is a continuous-time probabilistic process, it can be discretized in a grid where each position of the grid can take 0 or 1 values, 0 (1) indicating the absence (presence) of a particle or ion in that position (Fig. 1a and b).

Brownian motion is a mechanism of ion transportation implemented at the atomic level – it describes the process of particle deposition using a discrete sequence of events in time, and guarantees maximum accuracy in reproducing the structure of an aggregate. Fick's laws of diffusion are derived from Brownian motion and molecular diffusion: they represent a generalization of Brownian motion at the continuum limit. Thus, direct numerical simulation using DLA is more appropriate for micro- and nano-scale systems, for which the continuity assumption breaks down, while continuous Fick diffusion is more efficient when one considers macro-scale systems.

In what follows, the word particle is used interchangeably with ion, they both indicate the smallest building blocks that assemble together to form the final aggregates. Ions drift in the domain (grid) due to pure diffusion. Therefore, displacement of ions in the grid obeys a random uniform distribution; at each time of the process particles move from a position of the grid to one of its nearest neighbors, each of those positions being equally probable. We call Δx the size of the cells (pixels) of the grid. It is the size of the smallest features in the grid and corresponds to the resolution of the system. Since the motion of a particle in a grid takes discrete steps, in a macroscopic interpretation Δx is also the mean path length of the particles, i.e. the average distance traveled by a particle between successive collisions. Cases in which the mean path length is different from the resolution of the systems are discussed below. During motion, particles have a fixed velocity v. The mean time interval between collisions is thus . The nature of Δx, v and τ is determined by the energy of the system and its temperature.¹⁰ Consider the scheme in Fig. 1c.

5.1.1 Geometry. The portion of the silicon substrate (well) exposed to electroless deposition is represented at the bottom of the diagram as a line of nucleation sites Λ. Its length is w. The height of resist walls that delimit the well is h. The ratio determines the aspect ratio of the well. If the aspect ratio is lower than one, then we have a short geometry. In contrast, tall geometries have aspect ratios greater than one. Short geometry reproduces more accurately the real physical prototype and the experimental set up, where the dimension of lithographed patterns (10 μm) is larger than the thickness of the resist (∼1 μm). Here, we chose w and h as w = 100 pixels, h = 50 pixels, such that . Dimensions of the entire domain are 400 pixels (length) and 500 pixels (height).

5.1.2 Initial and boundary conditions. At the left and right boundaries of the domain, periodic boundary conditions are imposed. At the upper boundary of the domain, we enforce a bouncing boundary condition. Moreover, we assume that particle deposition occurs in an excess of solute, i.e. the number n of metal ions in the domain is generally high and it is maintained constant as n = 1000. At the initial time of growth, metal ions are positioned in a region of the domain at a distance l from the well.

5.1.3 Aggregation rules. Then, the system is left free to evolve. At any iteration, particles move in the grid by one pixel. After an iteration, a particle may find itself in two separate, mutually exclusive states: (i) isolated and (ii) contact state. In the isolated state (s = I) all sites surrounding the particle are 0 (the particle misses the aggregate). In the contact state (s = C) the minimum distance between the particle and the elements of the seed Λ is 0 (the particle hits the aggregate). The algorithm starts with the condition s = I and proceeds until the test s = C for at least one particle in the domain yields true. When a particle hits the aggregate, it is incorporated into the seed to form an aggregate of particles β if a randomly generated number q is such that q < p. In the simulations, the sticking probability p can be varied between 0 and 1. If p = 0, particles would indefinitely migrate in the domain never accumulating to the aggregate. If p = 1, particle deposition is deterministic, and particles would be captured by the aggregate anytime that they hit it. After a certain number of repetitions aggregates shall have the aspect depicted in Fig. 1d. The multi-branched arrangement of particles recalls the dendrite, fractal nature that electroless systems reveal under certain growth conditions. On changing p, one can modulate the size, shape and density of the aggregates.

5.1.4 Stop condition. Simulations are halted after ∼1.2 × 10⁷ iterations.

5.1.5 Considering the thermodynamic state of a system. Thermodynamic variables are implicitly contained in the model. Super-saturation, thermodynamic potentials and the temperature itself depend on the sole temperature field T of the system, assuming that pressure gradients are vanishingly small everywhere in the domain and that the number of ions n in the system is maintained constant over time (i.e. particle deposition proceeds in an excess of solute). Under these conditions, we can express the spatial variables of the simulation in terms of the sole T. Specifically: the root mean square distance (i.e. the variance) of a bolus of diffusing particles is

〈r²〉 = 〈x²〉 + 〈y²〉 = 4Dt (7)

where x, y, r are the coordinates of the bolus, and D is the molecular coefficient of diffusionwith k_b the Boltzmann constant, T the absolute temperature of the system, μ the viscosity, and a the diameter of the particles. Moreover, using results from the kinetic theory of gases, we can write

D = δx²/2τ (9)

where δx is the mean path length of an ion in solution and τ is the time between collisions. In molecular-scale systems, the kinetic energy of a particle with mass m and velocity v is

Combining these equations, one obtains:

δx = 2D/v = f(T) (11)

which is a function of time. Eqn (11) correlates the mean path length of an ion with the thermodynamic variables of a system and the temperature T. Since the mean path length is a variable of the model, one can change the temperature of the system by tuning δx: the thermodynamic variables of the system are lumped in the term δx.

5.1.6 Additional remarks. The model is a numerical diffusion limited aggregation simulation model – it reproduces aggregation of smaller ions into larger particles under the assumption of a process that is limited by diffusion, i.e., movement of ions is a random walk, and there are no preferential directions of motion in the domain. The model neglects external driving forces or fields including convective fields, electro-magnetic fields, the generation of mass in the domain through internal chemical or nuclear reactions. Nevertheless, assuming linearity, in a more sophisticated evolution of the model these effects can be easily included by adding a constant term to the analysis and considering superposition of individual stimuli. Since ions add up to the aggregate immediately upon contact, another assumption is that motion of ions is much slower than the chemical reaction of reduction of ions into metals at the liquid–solid interface, which is instantaneous. However, since different materials have different values of reaction kinetics – that in turn influences the rate k of growth of the metal – we introduced in the model an additional variable, p, comprised between 0 and 1, to account for this diversity. The assumption that k depends on the sole p is heuristically supported by the numerical diagram in Fig. 3b, where the mass of the system (i.e. the size of the aggregate N) varies linearly with p. Moreover, the model is 2D. The model accepts as inputs (i) the geometry of the system, (ii) the resolution of the system, i.e. discretization, (iii) the mean path length of ions (that depends on the thermodynamic variables of the system), (iv) the concentration of ions in solution, (v) the total length of the simulation. The output of the simulation is an aggregate of pixels, where the geometrical characteristics of the aggregate, measured through mathematical variables including fractal dimension, are indicative of the characteristics of the true particle assemblies.

5.2 Deriving the fractal dimension of numerical aggregates

Fractals are irregular, broken geometries for which conventional Euclidean metrics cannot provide a complete description.³⁰ Fractal dimension quantifies the complexity of a fractal as a change of detail to a change of scale. A possible method to calculate the fractal dimension of numerical aggregates is to make use of a density–density correlation functionwhere ρ(r) is a local density at the point r. It is equal to 1 if the point belongs to the fractal object and 0 otherwise. N is the number of particles forming the aggregate. In a discrete space where the objects can occupy only lattice points, ρ(r) is the density at the lattice site with coordinates r. One can think of c(r) as a probability density of finding two points belonging to a fractal at distance r between each other. The density–density correlation function follows the power law dependency³⁶ on r: c(r) ∼ r^−α, from which the relationship between the usual Euclidian dimension d, fractal dimension D_f and the density–density correlation function exponent α is D_f = d − α. In a computer simulation D_f is defined as the slope of linear fit of ln(c(r)) as a function of ln(r). In the implemented method, a rectangle of height r was incremented by a quantity of one lattice unit for each step starting from the base to reach the last particle of the aggregate. The particles progressively included in rectangles with height r were counted in order to estimate the related density values. The obtained values were then correlated in a log–log plot with dimension radius r and fractal dimension is evaluated by the difference between Euclidean dimension and graph slope as shown in Fig. 2.

5.3 Electroless deposition of metals on a substrate

In electroless deposition, metal ions in solution are transported and deposited as metals on a substrate. We consider here the case in which the substrate is silicon. Metal reduction on a silicon surface is regulated by the general law:

Meⁿ⁺ + Red = Me⁰ + Ox (13)

where Meⁿ⁺ are metal ions and Red represents the reducing agent – the reaction uses n electrons necessary to produce the atomic metal Me⁰. In the case of silver ions in solution, the reaction takes the form

Si + 2H₂O → SiO₂ + 4H⁺ + 4e⁻ (14)

at the anode (silicon oxidation) and

Ag⁺ + e⁻ → Ag⁰ (15)

at the cathode (silver reduction). In the reaction, the redox potentials are E₀₁ = −0.86 V, E₀₂ = 0.8 V.³⁷ In the case of gold ions in solution, the reaction takes the form

Si + 2H₂O → SiO₂ + 4H⁺ + 4e⁻ (16)

at the anode (silicon oxidation) and

Au³⁺ + 3e⁻ → Au⁰ (17)

at the cathode (gold reduction). Here the redox potentials are E₀₃ = −0.86 V, E₀₄ = 1.52 V.³⁸ Particle growth involves an initial nucleation phase with the formation of metallic nuclei and a steady state phase of growth. The shape and size of the final structures depend on the concentration of reducing agents, temperature, time of reaction and type of metals. Typical values of solute concentrations in solution are 0.15 M for hydrofluoric acid (HF), 0.1–5 mM for auric chloride (AuCl₃), and 0.1–5 mM for silver nitrate (AgNO₃). Electroless deposition enables production of nanoparticles, nanoparticle clusters and nanoparticle arrays where the mean size of the particles may be finely adjusted in the 5–100 nm range.

5.4 Synthesis of metal nanoparticle clusters by electroless deposition

The fabrication of aggregates of metal nanoparticles occurred in two steps: (i) the realization of micropatterns by optical lithography and (ii) the electroless deposition of silver nanoparticles within the patterns. De-ionized (D.I.) water (Milli-Q Direct 3, Millipore) was used for all experiments. Silver nitrate (AgNO₃), gold chloride (AuCl₃) and hydrofluoric acid (HF) were purchased from Sigma. All chemicals, unless mentioned otherwise, were of analytical grade and were used as received. P-type (100) silicon wafers were used as substrates. Substrates were accurately cleaned with acetone and ethanol to remove possible contaminants and then etched with a 20% in weight HF solution (hydrofluoric acid 50% RPE ACS-ISO, Carlo Erba Reagents) to eliminate the superficial native oxide and create the dangling bonds necessary for metal reduction. Using optical lithography (i) we obtained ordered arrays of holes in the positive resist S1813, with a diameter d = 10 μm and a depth h ∼ 1 μm. We then exposed patterned silicon surfaces to solutions of gold chloride and hydrofluoric acid for depositing gold nanoparticles, and silver nitrate and hydrofluoric acid for depositing silver nanoparticles, for 5, 20, 60 and 120 s. In both cases, the concentration ψ of metal salts in solution was chosen to be ψ(AgNO₃) = ψ(AuCl₃) = 5 mM, the concentration of hydrofluoric acid was set as ψ(HF) = 150 mM, and the temperature of the system was maintained constant as T = 20 °C throughout the whole duration of the experiments. Resulting nanostructures were imaged using scanning electron microscopy (SEM) and atomic force microscopy (AFM).

5.5 SEM sample characterization

SEM images of metal nanoparticle clusters were captured using a Dual Beam (SEM-FIB) – FEI Nova 600 NanoLab system. During the acquisitions beam energies of 5 and 15 keV, and corresponding electron currents of 0.98 pA and 0.14 nA, were used. In some cases the mode 2 configuration was set, through which images could be magnified over 2.5 × 10⁶ times to achieve ultrahigh resolution.

5.6 AFM sample characterization

Atomic force microscopy (integrated Raman AFM system, Alpha 300 RA, Witec) was used for sample characterization. All the measurements were performed in a dry environment in intermittent contact mode over a sampling area of 1 × 1 μm². Room temperature was held fixed for all the acquisitions. Ultra-sharp Si probes with a nominal tip radius less than 5 nm were used for achieving high resolution. Multiple measurements were done in different scan directions to avoid artefacts. At least four images in height mode (trace and retrace) were recorded for each sample. The images had a resolution of 512 × 512 points and were acquired at a scanning rate of 1 Hz. Images were processed using either flattening or plane fit according to the relief characteristics, with the minimal polynomial order needed. The characteristic average surface roughness S_a was thus deconvolved for each substrate. Fast Fourier Transform (FFT) algorithms were used for data processing and fractal extraction of the characteristic dimension of the sample surface.

5.7 Deriving the fractal dimension of metal nanoparticles

AFM profiles of metal nanoparticle clusters were processed using the algorithms developed and described in ref. 34. We derived the characteristic power density function for gold and silver nanoparticle clusters at 60 s after immersion in the electroless solution. In a log–log plot, the power spectrum density appears as a line with a slope β. The slope β is related to the Hurst parameters as β = 2(H + 1). The fractal dimension D_f of the surface can be equivalently derived as or D_f = 3 − H. The fractal dimension D_f of a surface ranges from 2 (Euclidean dimension of a flat surface) to 3 (representing an extremely rough surface).

5.8 UV spectrophotometric measurement of metal nanoparticle growth kinetics

We used indirect measurements to determine the growth dynamics of Ag and Au nanoparticles. We put in an aqueous solution nanoporous silicon submicrometric particles with (i) silver (AgNO₃) and (ii) gold (AuCl₃) salts. The resulting electroless reactions were examined using UV/Vis spectrophotometry.^39–42 Metal ions in solution react with the silicon sub-micrometric particles forming metal nanoparticles in the porous matrix. We prepared two samples of HF (0.15 M) plus either (i) AgNO₃ 1 mM or (ii) AuCl₃ 1 mM in one ml of DI water. Addition of silicon particles in the solutions caused the start of the electroless reaction of reduction. We added nearly ∼10¹⁷ silicon sub-micrometric particles in each milliliter sample. Then, we used a LAMBDA 25 UV/Vis spectrophotometer to determine the growth dynamics of Au and Ag nanoparticles as the variation of absorbance of solutions at different times. Absorbance was measured at wavelengths λ₁ = 560 nm for gold and λ₂ = 460 nm for silver, which are the values of wavelength for which spectral intensities are maximum in the considered ranges of frequencies. Rate constants of metal nanoparticle growth were determined as the slope of the absorbance versus time diagrams at the initial time of metal deposition. Nanoporous silicon sub-micrometric particles were fabricated through electrochemical porosofication of silicon as described in ref. 22.

Conflicts of interest

There are no conflicts to declare.

References

1. A. R. Tao, S. Habas and P. Yang, Small, 2008, 4, 310–325.
2. S. R. Schofield, P. Studer, C. F. Hirjibehedin, N. J. Curson, G. Aeppli and D. R. Bowler, Nat. Commun., 2013, 4, 1649.
3. T. Qiu and P. Chu, Mater. Sci. Eng., R, 2008, 61, 59–77.
4. F. Gentile, M. Coluccio, A. Toma, E. Rondanina, M. Leoncini, F. De Angelis, G. Das, C. Dorigoni, P. Candeloro and E. Di Fabrizio, Microelectron. Eng., 2012, 98, 359–362.
5. M. L. Coluccio, F. Gentile, M. Francardi, G. Perozziello, N. Malara, P. Candeloro and E. Di Fabrizio, Sensors, 2014, 14, 6056–6083.
6. H. Han, Z. Huang and W. Lee, Nano Today, 2014, 9, 271–304.
7. F. Gentile, M. L. Coluccio, R. P. Zaccaria, M. Francardi, G. Cojoc, G. Perozziello, R. Raimondo, P. Candeloroa and E. D. Fabrizio, Nanoscale, 2014, 6, 8208–8225.
8. D. Fine, A. Grattoni, R. Goodall, S. S. Bansal, C. Chiappini, S. Hosali, A. L. van de Ven, S. Srinivasan, X. Liu, B. Godin, L. Brousseau III, I. K. Yazdi, J. F. Moure, E. Tasciotti, H.-J. Wu, Y. Hu, S. Klemm and M. Ferrari, Adv. Healthcare Mater., 2013, 2, 625.
9. Y.-W. Fan, F.-Z. Cui, S.-P. Hou, Q.-Y. Xu, L.-N. Chen and I.-S. Lee, J. Neurosci. Methods, 2002, 120, 17–23.
10. M. Saltzmann, Drug Delivery, Oxford University Press, Oxford, 2001.
11. A. Einstein, Ann. Phys., 1905, 17, 549–560.
12. Z. Zhang and M. Lagally, Science, 1997, 276, 377–383.
13. D. C. Wyld, J. Zizka and D. Nagamalai, Advances in Computer Science, Engineering and Applications, New Delhi, India, 2012.
14. T. Witten and L. Sander, Phys. Rev. Lett., 1981, 47, 1400–1403.
15. T. Witten and L. Sander, Phys. Rev. B: Condens. Matter Mater. Phys., 1983, 27, 5686–5697.
16. Z. Racz and T. Vicsek, Phys. Rev. Lett., 1983, 51, 2382–2385.
17. P. Meakin, Phys. Rev. B: Condens. Matter Mater. Phys., 1984, 30, 4207–4214.
18. F. Gentile, M. L. Coluccio, P. Candeloro, M. Barberio, G. Perozziello, M. Francardi and E. Di Fabrizio, J. Vac. Sci. Technol., B: Nanotechnol. Microelectron.: Mater., Process., Meas., Phenom., 2014, 32, 031804.
19. R. Pastor-Satorras and J. Wagensberg, Physica A, 1996, 224, 463–479.
20. G. Marinaro, R. La Rocca, A. Toma, M. Barberio, L. Cancedda, E. Di Fabrizio, P. Decuzzi and F. Gentile, Integr. Biol., 2015, 7, 184–197.
21. H. Foll, M. Christophersen, J. Carstensen and G. Hasse, Mater. Sci. Eng., 2002, 39, 93–141.
22. F. De Angelis, A. Pujia, C. Falcone, E. Iaccino, C. Palmieri, C. Liberale, F. Mecarini, P. Candeloro, L. Luberto, A. de Laurentiis, G. Das, G. Scala and E. Di Fabrizio, Nanoscale, 2010, 2, 2230–2236.
23. B. Zheng, L. P. Wong, L. Y. L. Wu and Z. Chen, Int. J. Electrochem., 2016, 2016, 1–11.
24. T. Yang and Y. Han, Cryst. Growth Des., 2016, 16, 2850–2859.
25. J. G. Velasco, Electroanalysis, 1997, 9, 880–882.
26. S. Zhang, H. Bai, J. Pi, P. Yang and J. Cai, Anal. Chem., 2015, 87, 4797–4805.
27. M. E. J. Newman, SIAM Rev., 2003, 45, 167–256.
28. I. Mannelli, M. Minunni, S. Tombelli and M. Mascini, Biosens. Bioelectron., 2003, 18, 129–140.
29. S. Glotzer, Nature, 2012, 481, 450–452.
30. F. Gentile, E. Battista, A. Accardo, M. Coluccio, M. Asande, G. Perozziello, G. Das, C. Liberale, F. De Angelis, P. Candeloro, P. Decuzzi and E. Di Fabrizio, Microelectron. Eng., 2011, 88, 2537–2540.
31. P. Decuzzi, R. Pasqualini, W. Arap and M. Ferrari, Pharm. Res., 2009, 26, 235–243.
32. A. Carpinteri and N. Pugno, Nat. Mater., 2005, 4, 421–423.
33. V. Amendola, O. M. Bakr and F. Stellacci, Plasmonics, 2010, 5, 85–97.
34. F. Gentile, L. Tirinato, E. Battista, F. Causa, C. Liberale, E. Di Fabrizio and P. Decuzzi, Biomaterials, 2010, 31, 7205–7212.
35. F. Gentile, R. Medda, L. Cheng, E. Battista, P. Scopelliti, P. Milani, E. Cavalcanti-Adam and P. Decuzzi, Sci. Rep., 2013, 3, 1461.
36. F. Gentile, M. Monteferrante, L. Chiodo, A. Toma, M. L. Coluccio, G. Ciccotti and E. Di Fabrizio, Molecular Physics, 2014, 112, 1375–1388.
37. D. Goia and E. Matijevic, New J. Chem., 1998, 22, 1203–1215.
38. S. Yae, N. Nasu, K. Matsumoto, T. Hagihara, N. Fukumuro and H. Matsuda, Electrochim. Acta, 2007, 53, 35–41.
39. S. A. Wadhwani, U. U. Shedbalkar, R. Singh, P. Vashisth, V. Pruthi and B. A. Chopad, Indian J. Microbiol., 2016, 56, 439–444.
40. M. Luty-Blocho, K. Pacławski, W. Jaworski, B. Streszewski and K. Fitzner, Prog. Colloid Polym. Sci., 2011, 138, 39–44.
41. K. Esumi, T. Hosoya, A. Suzuki and K. Torigoe, J. Colloid Interface Sci., 2000, 226, 346–352.
42. S. A. AL-Thabaiti, F. M. Al-Nowaiser, A. Y. Obaid, A. O. Al-Youbi and Z. Khan, Colloids Surf., B, 2008, 67, 230–237.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c8na00040a

---