Nanocrystal growth via the precipitation method

TL;DR

This paper presents a mathematical model for nanocrystal growth via precipitation, including single and multiple particle systems, validated against experimental data, and offers insights for process optimization.

Contribution

It introduces a comprehensive mathematical framework for nanocrystal growth, extending from single particles to multiple particles, including Ostwald ripening, with analytical solutions and experimental validation.

Findings

Excellent agreement with experimental data

Single controlling group influences growth process

Model can describe bimodal distributions

Abstract

A mathematical model to describe the growth of an arbitrarily large number of nanocrystals from solution is presented. First, the model for a single particle is developed. By non-dimensionalising the system we are able to determine the dominant terms and reduce it to the standard pseudo-steady approximation. The range of applicability and further reductions are discussed. An approximate analytical solution is also presented. The one particle model is then generalised to $N$ well dispersed particles. By setting $N=2$ we are able to investigate in detail the process of Ostwald ripening. The various models, the $N$ particle, single particle and the analytical solution are compared against experimental data, all showing excellent agreement. By allowing $N$ to increase we show that the single particle model may be considered as representing the average radius of a system with a large number of particles. Following a similar argument the $N=2$ model could describe an initially bimodal distribution. The mathematical solution clearly shows the effect of problem parameters on the growth process and, significantly, that there is a single controlling group. The model provides a simple way to understand nanocrystal growth and hence to guide and optimise the process.

Figures

Tables

Table 1. Table 1: Physical parameters for the cadmium selenide (CdSe) nanoparticle synthesis method used by Peng et al .  [ 23 ] . The parameters in italics are not given in explicitly and are thus obtained via a fitting approach.

Quantity	Symbol	Value	Units
Universal gas constant	$R_G$	8.31	$J{\mathrm{mol}}^{-1}{K}^{-1}$
Density	${\rho }_p$	5816	$\mathrm{kg}{m}^{-3}$
Molar volume	$V_M$	$3.29\times {10}^{-5}$	$m^3{\mathrm{mol}}^{-1}$
Molar mass	$M_p$	0.19	$\mathrm{kg}{\mathrm{mol}}^{-1}$
Solution temperature	$T$	573.15	$K$
Surface energy	$\sigma$	0.44	$J{m}^{-2}$
Capillary length	$\alpha$	$6.00\times {10}^{-9}$	$m$
Initial bulk concentration	$c_0^{\ast }$	55.33	$\mathrm{mol}{m}^{-3}$
Volume of the liquid	$V$	$7.21\times {10}^{-6}$	$m^3$
Maximum particle radius	$r_m$	$3.78\times {10}^{-9}$	$m$
Diffusion coefficient	$D$	$3.01\times {10}^{-18}$	$m^2{s}^{-1}$
Reaction rate	$k$	$7.97\times {10}^{-10}$	$m{s}^{-1}$
Solubility of bulk material	$s_{\infty }^{\ast }$	$5.53\times {10}^{-2}$	$\mathrm{mol}{m}^{-3}$
Population density	$N_0$	$8.04\times {10}^{21}$	$\mathrm{No}{m}^{-3}$

Equations

$$s^{\ast }=s_{\infty }^{\ast }\exp \left(\frac{2\sigma V_M}{r_p^{\ast }{R}_{G}T}\right)\equiv s_{\infty }^{\ast }\exp \left(\frac{\alpha }{r_p^{\ast }}\right),$$

$$\frac{\partial c^{\ast }}{\partial t^{\ast }}=\frac{D}{r^{\ast 2}}\frac{{\partial }^{\ast }}{\partial r^{\ast }}\left(r^{\ast 2}\frac{\partial c^{\ast }}{\partial r^{\ast }}\right).$$

$$c^{\ast }(r_p^{\ast },t^{\ast })=c_i^{\ast }(t^{\ast }),c^{\ast }(r_p^{\ast }+{\delta }^{\ast },t^{\ast })=c_b^{\ast }(t^{\ast }),c^{\ast }(r^{\ast },0)=c_{b,0}^{\ast }\text{for r∗>rp∗},$$

$$J=4\pi r^{\ast 2}D\frac{\partial c^{\ast }}{\partial r^{\ast }}.$$

$$J=4\pi r^{\ast 2}k(c_i^{\ast }-s^{\ast }),$$

$$c_{i}^{*}=s^{*}+\frac{D}{k}\frac{\partial c^{*}}{\partial r^{*}}\bigg{|}_{r^{*}=r_{p}^{*}}\,,$$

$$\frac{1}{N_0}{M}_p{c}_0^{\ast }=M_p{c}_b^{\ast }(t)\left[\frac{1}{N_0}-\frac{4\pi }{3}{r_p^{\ast }}^3\right]+\frac{4\pi }{3}{\rho }_p{r_p^{\ast }}^3$$

$$c_b^{\ast }(t)\approx c_0^{\ast }-\frac{4\pi N_0}{3V_M}{r_p^{\ast }}^3.$$

$$V_{M}J=\frac{d}{d{t}^{\ast }}\left(\frac{4}{3}\pi {r_p^{\ast }}^3\right)=4\pi {r_p^{\ast }}^2\frac{d{r}_p^{\ast }}{d{t}^{\ast }}.$$

$$\frac{\mathop{}\!\mathrm{d}r_{p}^{*}}{\mathop{}\!\mathrm{d}t^{*}}=V_{M}D\frac{\partial c^{*}}{\partial r^{*}}\bigg{|}_{r^{*}=r_{p}^{*}}\,.$$

$$r=\frac{r^{\ast }}{r_{p,0}^{\ast }},r_p=\frac{r_p^{\ast }}{r_{p,0}^{\ast }},t=\frac{t^{\ast }}{{\tau }^{\ast }},c=\frac{c^{\ast }-s_0^{\ast }}{\Delta c},s=\frac{s^{\ast }-s_0^{\ast }}{\Delta c},$$

$$\epsilon \frac{\partial c}{\partial t}=\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial c}{\partial r}\right),$$

$$\displaystyle\frac{\mathop{}\!\mathrm{d}r_{p}}{\mathop{}\!\mathrm{d}t}=\frac{\partial c}{\partial r}\bigg{|}_{r=r_{p}}\,,$$

$$\displaystyle c(r_{p},t)=s+\operatorname{Da}\frac{\partial c}{\partial r}\bigg{|}_{r=r_{p}}\,,\qquad c(r_{p}+\delta,t)=\bar{c}_{0}-\beta r_{p}^{3},\,,\qquad$$

$$c(r,0)=1,r_p(0)=1,$$

$$s=s_{\infty }\exp (\omega /r_p)-s_0,$$

$$\epsilon =V_M\Delta c,\delta =\frac{{\delta }^{\ast }}{r_{p,0}^{\ast }},\mathrm{Da}=\frac{D}{k{r}_{p,0}^{\ast }},\omega =\frac{\alpha }{r_{p,0}^{\ast }},\beta =\frac{4\pi N_0{r}_{p,0}^{\ast 3}}{3V_M},{\bar{c}}_0=\frac{c_0^{\ast }-s_0^{\ast }}{\Delta c}.$$

$$\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial c}{\partial r}\right)=0.$$

$$c=-\frac{A}{r}+B,$$

$$A=\frac{r_p^2(r_p+\delta )(c_b-s)}{r_p\delta +\mathrm{Da}(r_p+\delta )},B=s+A\left(\frac{1}{r_p}+\frac{\mathrm{Da}}{r_p^2}\right).$$

$$A=\frac{r_p^2(c_b-s)}{r_p+\mathrm{Da}},B=c_b.$$

$$\frac{\mathrm{d}{r}_p}{\mathrm{d}t}=\frac{c_b-s}{\mathrm{Da}+r_p}=\frac{{\bar{c}}_0-\beta r_p^3+s_0-s_{\infty }\exp (\omega /r_p)}{\mathrm{Da}+r_p}.$$

$$r_p=\frac{r_m}{2}\frac{\left[1+2f(r_{p0})\exp \left(\frac{t-t_0}{G}\right)-\sqrt{-3+12f(r_{p0})\exp \left(\frac{t-t_0}{G}\right)}\right]}{\left[-1+f(r_{p0})\exp \left(\frac{t-t_0}{G}\right)\right]},$$

$$G=\frac{1}{6ab}\frac{ak+bD}{akbD},a^3=V_m(c_0^{\ast }-c_{eq}^{\ast }),b^3=\frac{4}{3}\pi N_0.$$

$$\frac{1}{N_0}{M}_p{c}_0^{\ast }=M_p{c}_b^{\ast }(t^{\ast })\left[\frac{1}{N_0}-\frac{4\pi }{3}\sum _{i=1}^N{r_i^{\ast }}^3\right]+\frac{4\pi {\rho }_p}{3}\sum _{i=1}^N{r_i^{\ast }}^3$$

$$c_b^{\ast }(t^{\ast })\approx c_0^{\ast }-\frac{4\pi N_0}{3V_M}\sum _{i=1}^N{r}_i^{\ast 3}.$$

$$\frac{\mathrm{d}{r}_i}{\mathrm{d}t}=\frac{c_b-s_i}{\mathrm{Da}+r_i}=\frac{{\bar{c}}_0-\bar{\beta }{\sum }_{i=1}^N{r}_i^3+{\bar{s}}_0-s_{\infty }\exp (\omega /r_i)}{\mathrm{Da}+r_i}$$

Full text

Keywords: Nanocrystal growth; Size focussing; Ostwald ripening; Mathematical model

1 Introduction

Nanoparticles (NPs) are small units of matter with dimensions in the range 1-$100\text{\,}\mathrm{n}\mathrm{m}$. They exhibit many advantageous, size-dependent properties such as magnetic, electrical, chemical and optical, which are not observed at the microscale or larger [17, 9, 4, 30]. Consequently the ability to produce monodisperse particles that lie within a controlled size distribution is critical.

There exist a number of NP synthesis methods, including gas phase and solution based synthesis techniques. Although the first method can produce large quantities of nanoparticles, it produces undesired agglomeration and nonuniformity in particle size and shape. Precipitation of NPs from solution avoids these problems and is one of the most widely used synthesis methods [16]. The typical strategy is to cause a short nucleation burst in order to create a large number of nuclei in a short space of time, and the seeds generated are used for the latter particle growth stage. Thus the temporal separation of nucleation and growth occurs, as proposed by La Mer and Dinegar [11, 12], is applied. The resulting system consists of varying sized particles. Small NPs are more unstable than larger ones and tend to grow or dissolve faster. Thus at relatively high monomer concentrations size focussing occurs (leading to monodispersity). When the monomer concentration is depleted by the growth some smaller NPs shrink and eventually disappear while larger particles continue to grow, thus leading to a broadening of the size distribution (Ostwald ripening).

The particle size distribution (PSD) can be refocused by changing the reaction kinetics. For example, Peng et al. [24] observed size focusing during Cadmium Selenide growth following the injection of additional solute. Bastús et al. [2, 3] were also able to induce size focusing of gold and silver nanoparticles by the addition of extra solute and adjusting the temperature and pH. This type of technique for size focussing is still rather ad hoc in that the precise relationships between particle growth, system conditions and the final PSD are not fully understood [26]. Hence, in practice, the optimal reaction conditions are usually ascertained empirically or intuitively.

In the 1960’s Lifshitz and Slyozov [13] and, independently, Wagner [35] were amongst the first to provide theoretical descriptions of Ostwald ripening. Their classical theory, hereafter referred to as LSW theory, consisted of a system of three coupled equations: a growth equation for a single particle, a continuity equation for the PSD and a mass conservation expression for the concentration. They solved the model to obtain pseudo-steady-state asymptotic solutions for the average particle radius and PSD. Lifshitz and Slyozov [13] focused on diffusion-limited growth, where growth is limited by the diffusion of reactants to the particle surface, while Wagner [35] considered growth limited by the reactions at the particle surface (i.e. reaction-limited growth). In fact recent work described by Myers and Fanelli [19] have shown that, within the restrictions of the steady-state assumption, the model cannot distinguish between diffusion or reaction driven growth, so both approaches are equally valid. For this reasons authors using either mechanism, or both, have been equally successful in approximating experimental data.

Experimental studies on NP growth [6, 22, 27] show that LSW theory may provide good predictions for the particle size but the observed PSDs were typically broader and more symmetric. Possible explanations for this disparity is that LSW theory does not account for the finite volume of the coarsening phase $\phi$, and that it assumes a particle’s growth rate is independent of its surroundings. In addition, LSW theory does not indicate how long it takes to reach the final state. A further issue is that it purports to describe the dynamics in the initial stages of the growth process. In [19] it is proven that the pseudo-steady solution does not hold for small times.

Many studies have modified and built on the pioneering analysis of LSW theory. Ardell [1] and Sarian and Weart [25] extended LSW theory to systems for non-zero $\phi$ (i.e. the mean distance between particles is finite). Several authors [5, 33, 34] have addressed the shortcomings of LSW theory by statistically averaging the diffusional interaction of a particle of a given size with its surroundings to demonstrate that the resulting PSD becomes broader and more symmetric with increasing $\phi$. The inclusion of stochastic effects, due to temperature and changes in concentration, in the modified population balance model of Ludwig et al. [15] led to broader PSDs in line with experimental data. The population balance approach of Iggland and Mazzotti [10] was used to examine the evolution of non-spherical particles at the beginning of growth.

Most of the above studies were in relation to micron or larger-sized particles. As measurement techniques have advanced many researchers have applied LSW theory and the related modifications to the study of nanoparticle growth. Talapin et al. [29] used a Monte Carlo approach to simulate the evolution of a nanoparticle PSD subject to diffusion-limited growth, reaction-limited growth and mixed diffusion-reaction growth. In contrast to other treatments, their simulations gave PSDs narrower than those predicted by LSW theory. This was explained by the fact that they considered much smaller particles. Their main conclusion was that Ostwald ripening occurs much more rapidly for nanoparticles while PSDs are narrower than in their microscale counterparts. Similarly, Mantzaris [16] used a population balance formulation and a moving boundary algorithm to study the diffusion and reaction-limited growth regimes.

Another issue which is particularly relevant in the context of nanoparticles is the applicability of the Ostwald-Freundlich condition (OFC) which relates the radius of the particle, $r_{p}^{*}$, to its solubility, $s^{*}$. This condition can be written as

[TABLE]

where $s_{\infty}^{*}$ is the solubility of the bulk material, $\sigma$ the interfacial energy, $R_{G}$ the universal gas constant, $T$ the absolute temperature. The capillary length $\alpha=2\sigma V_{M}/(R_{G}T)$ defines the length scale below which curvature-induced solubility is significant [29]. This equation shows that the particle solubility increases as the size decreases (which promotes Ostwald ripening). One approximation to the OFC is to assume that the exponential term in (1) can be linearised to give the two term expression $s^{*}\approx s_{\infty}^{*}(1+\alpha/r_{p}^{*})$[13, 14, 28, 35]. Obviously this expansion, which is based on $\alpha/r_{p}^{*}$, is invalid for nanoparticles where capillary length is of the same order of magnitude as the particle radius [19]. Mantzaris [16] used an expansion for the exponential term in the OFC with $n$ terms and showed that increasing $n$ led to higher average growth rates and a narrowing of the PSD. However, when comparing his simulation to experimental data for CdSe nanoparticles from [23], he applied the linearised version for the solubility. Talapin et al. [29], noting that for nanoparticles of the order 1-$5\text{\,}\mathrm{n}\mathrm{m}$ the linearised OFC may be incorrect, applied the full condition.

In the following we begin by analysing the growth of a single particle. This is the basic building block for more complex models. The treatment leads to equations similar to those of standard LSW theory, however we arrive at them following a non-dimensionalisation which highlights dominant terms and those which may be formally neglected. In this way we can ascertain which standard assumptions are appropriate and, more importantly, which are not. Under conditions which appear easily satisfied for nanocrystal growth the governing ordinary differential equation has an explicit solution, in the form $r_{p}=r_{p}(t)$ and also shows that the growth is controlled by a single parameter which may be calculated by comparison with experiment. This section closely follows the work described in [19]. The single particle model is obviously incapable of reproducing Ostwald ripening, where larger particles grow at the expense of smaller ones. Consequently we then generalise the model to deal with a large number of particles. In the results section we compare the analytical solution with that of a full numerical solution and experimental data for the growth of a single particle and show excellent agreement between all three. By setting the number of particles to two in the general model we are able to clearly demonstrate Ostwald ripening. Simulations with $N=$10 and 1000 particles demonstrates that increasing $N$ leads to increasingly good agreement between the average radius and that predicted by the single particle model. The single particle model may thus be considered as a viable method for predicting the evolution of the average radius of a group of particles.

2 Growth of a single particle

As shown in Figure 1, we initially focus on a single, spherical nanoparticle, with radius $r_{p}^{*}$ in a system of particles. The ∗ notation represents dimensional quantities. The assumption is that particles are separated at large but finite distances compared to their radius. Their morphologies remain nearly spherical and particle aggregation is neglected. Thus, the mass flow from each particle can be represented as a monopole source located at the center of the particle [32] and the problem becomes radially symmetric. We assume the standard La Mer model [11], such that there has been a short nucleation burst and the system is now in the period of growth.

The monomer concentration, $c^{*}$, is described by the classical diffusion equation in spherical coordinates

[TABLE]

This holds in the diffusion layer $[r_{p}^{*},r_{p}^{*}+\delta^{*}]$ where $r^{*}$ is distance from the centre of the particle, $t^{*}$ is time and $D$ is the constant diffusion coefficient. To conform with standard literature (see for example [16, 28, 31]), we have included a diffusion layer of length $\delta^{*}$ around the particle, where the concentration adjusts from the value at the particle surface to the value in the far-field. Equation (2) is subject to

[TABLE]

where $c^{*}_{i}$ is the concentration adjacent to the particle surface, $c_{b}^{*}$ is the concentration in the far-field and $c_{b,0}^{*}$ is a constant describing the initial concentration when the solution is well-mixed. The monomer concentration in the far-field $c_{b}^{*}(t^{*})$ will be derived via mass conservation. The value at the particle surface $c_{i}^{*}$ is very difficult to measure [28], hence it is standard to work in terms of the particle solubility.

The particle solubility $s^{*}$ (with the same dimensions as concentration) is given by the Ostwald–Freundlich condition (1). If $s^{*}<c_{b}^{*}$ then monomer molecules diffuse from the bulk towards the particle to react with the surface and the particle grows, whereas if $s^{*}>c_{b}^{*}$ the particle shrinks.

In order to determinate an expression for the concentration at the particle surface, we consider two equivalent expressions for the mass flux at the particle surface, $J$. Firstly, Fick’s first law states that the flux of monomers passing through a spherical surface of radius $r^{*}$ is

[TABLE]

At the surface of the sphere the flux must also follow a standard first order reaction equation

[TABLE]

where $k$ is the reaction rate, which is assumed to be constant for both growth and dissolution contributions. Equating (4) with (5) gives

[TABLE]

which defines the concentration $c_{i}^{*}$ for the surface condition of (3).

To complete the boundary conditions in the system, we require an expression for the time-dependent bulk concentration, $c_{b}^{*}(t^{*})$. The mass of monomer in the system is constant because the bulk material is assumed well-mixed during the entire process. Mass conservation of the monomer atoms in the particle and surrounding solution is then

[TABLE]

where $\rho_{p}$ is density, $M_{p}$ is molar mass and $N_{0}$ the population density. Since particle occupy a small region in comparison to their radius $4\pi N_{0}{r_{p}^{*}}^{3}/3\ll 1$, also the molar volume $V_{M}=M_{p}/\rho_{p}$. Equation (7) then leads to

[TABLE]

which will be used for the far-field concentration in (3).

As stated earlier, the diffusion equation must be solved on a domain $r^{*}>r_{p}^{*}$, where the particle radius is an unknown function of time. The flux of monomer to the particle is responsible for the particle growth

[TABLE]

Eliminating $J$ between (9) and (4) yields

[TABLE]

This is subject to the initial condition $r_{p}^{*}(0)=r_{p,0}^{*}$, where $r_{p,0}^{*}$ is the initial particle radius.

The governing system is now fully defined and consists of equation (2), subject to the initial and boundary conditions (3), where $c^{*}_{i}$ is defined by (6) and $c_{b}^{*}$ by (8), and the unknown particle radius satisfies (10). As there is no analytical solution to the system, we proceed to simplify the problem and use numerical approximations in order to understand the behaviour of the solution.

2.1 Nondimensionalisation

A complex mathematical model can often be reduced to a simpler form by estimating the relative magnitude of terms. A standard way to achieve this is by writing the system in non-dimensional form. Here the model is nondimensionalised via

[TABLE]

where $\Delta c=c^{*}_{b,0}-s^{*}_{0}$ represents the driving force for particle growth and $s^{*}_{0}=s_{\infty}^{*}\exp{(\alpha/r_{p,0}^{*})}$ is the initial particle solubility. The concentration and growth equations yield two possible time scales $\tau^{*}_{D}\sim r_{p,0}^{*2}/D$ and $\tau^{*}_{R}\sim r_{p,0}^{*2}/(V_{M}Dc_{0})$, respectively. To focus on particle growth we choose the growth time scale $\tau^{*}=\tau^{*}_{R}$ and the governing system is now transformed to

[TABLE]

where

[TABLE]

and

[TABLE]

The above system contains a number of nondimensional groups. The first, $\varepsilon$, is generally very small for nanoparticle growth. For example, Peng et al. [23] studied Cadmium Selenide nanoparticles, with a capillary length of 6nm and initial radii in the range $1-100\,$nm, so that $\varepsilon=\mathcal{O}(10^{-3})$. In general it should be expected that $\varepsilon\ll 1$. If we look at the time scales, we see that $\tau^{*}_{D}/\tau^{*}_{G}=V_{M}\Delta c=\varepsilon\ll 1$. Physically, this indicates that growth is orders of magnitude slower than the diffusion time scale, that is, the concentration adjusts much faster than growth occurs and so the system can be considered as pseudo-steady. In terms of the mathematical model, this means that the time derivative can be omitted from the concentration equation, but since time also enters into the problem through the definitions of $r_{p}$ and $c_{b}$ this is a pseudo-steady-state situation rather than a true steady-state.

The parameter $\operatorname{Da}$ is an inverse Damköhler number measuring the relative magnitude of diffusion to surface reactions [16]. In the past similar models have been simplified by considering diffusion-limited growth ($\operatorname{Da}\ll 1)$ or surface reaction limited growth ($\operatorname{Da}\gg 1$). In practice both mechanisms play a role. In [19] it is shown that the diffusion limited case requires either $k\to 0$, which results in zero growth, or the concentration adjacent to the particle surface matches the solubility $c_{i}\simeq s^{*}$ throughout the process. Similary the reaction driven growth requires $D\to 0$, which again indicates no growth, or $c_{i}\simeq c_{b}$ throughout the process. Therefore, we will place no restrictions on $\operatorname{Da}$.

A common simplification is to assume $\omega\ll 1$ which reduces the OFC, (1), to a constant $s^{*}=s_{\infty}^{*}$ or a linear approximation is used, see [14, 28]. This significantly simplifies the analysis. However, for particles that have just nucleated or very small nanoparticles $\omega$ is not small and the simplification is not appropriate. Despite the large errors in the prediction of $s^{*}$ caused by the small $\omega$ assumption authors obtain good matches to data. In [19] it is shown that this is because the pseudo-steady model is not valid for early times when the particle is small. By the time the model is valid so is the linearisation. Basically, the variation of $s^{*}$ plays a minor role in the study of the growth of a single nanoparticle. However, this is not the case with multiple particles where Ostwald ripening is driven by the delicate balance between the bulk concentration and the particle solubility.

The reason why the pseudo-steady model is invalid at small times is due to the thickness of the boundary layer $\delta(t)$. The model involves the assumption $\delta(t)\gg r_{p}$ yet initially, when the fluid is well-mixed $\delta(0)=0$. Only when the boundary layer is sufficiently thick is it reasonable to apply the pseudo-steady model. In [19] it is shown through comparison with experiments that the initial stage can last for the order of 100s. The shift to the pseudo-steady model can often be identified simply by looking at the trend in the data. In the following we will present the model with the full OFC and then an approximation where it is neglected. We will also neglect early data points when matching to experimental data.

2.2 Pseudo-steady state solution

Since $\varepsilon=\mathcal{O}(10^{-3})$ and all variables have been scaled to be $\mathcal{O}(1)$ then neglecting terms of order $\varepsilon$ should result in errors of the order 0.1%. Consequently, we neglect the time derivative in the diffusion equation and obtain the pseudo-steady state form

[TABLE]

After integrating and applying the boundary conditions we obtain

[TABLE]

where

[TABLE]

There is no way to calculate $\delta(t)$ in the pseudo-steady approach. A time-dependent treatment, such as that described in [18] is required. Hence the standard method is to assume $r_{p}\ll\delta$, which reduces (20) to

[TABLE]

Substituting (19) into the Stefan condition (13) leads to

[TABLE]

Hence, the problem has been reduced to the solution of a single first-order ordinary differential equation for $r_{p}$. It is a highly nonlinear equation which must be solved numerically. The assumption that $r_{p}\ll\delta$ means it only holds for relatively large times. Approximate solutions, in various limits, may be found in the literature. For example if we take $c_{b}$ constant and $\omega$ sufficiently small for the linear approximation to the exponential to hold then equation (22) may be integrated in the limits of large and small $\operatorname{Da}$. In [19] it is shown that for sufficiently large times, for a single particle, the variation of $e^{\omega/r_{p}}$ does not affect the solution in which case equation (22) may be integrated analytically to find an implicit solution of the form $t=t(r)$. By identifying negligible terms they are able to invert this to find an explicit solution, $r=r(t)$ which depends on a single parameter,

[TABLE]

where $r_{m}$ is the experimental maximum radius, $t_{0}$ is the time at which the second growth stage is judged to have begun $r_{p0}$ the radius at this time and $f(r_{p0})=(r_{m}^{2}+r_{m}r_{p0}+r_{p0}^{2})/(r_{m}-r_{p0})^{2}$. If $t_{0}$ is greater than the true value, this should not affect the results. This is discussed in further detail in [19]. The unknown parameter $G$ is defined as

[TABLE]

Its value is obtained by comparison with experimental data. Once $G$ is determined, then the diffusion coefficient ($D$), the reaction rate ($k$), the solubility of the bulk material ($s_{\infty}$) and population density ($N_{0}$) may be systematically retrieved. In [19] it is stated that $ak\approx bD$, hence $G\approx 1/(3a^{2}bk)=1/(3ab^{2}D)$. Further, since $c_{0}^{*}\gg c_{eq}^{*}$ a reasonable approximation is $a^{3}=V_{m}c_{0}^{*}$. Growth stops when the maximum radius $r_{m}^{*}=a/b$ is achieved.

3 Evolution of a system of $N$ particles

We now extend the single particle model to an arbitrarily large system of particles. The particle radii, initial radii and solubilities are denoted $r^{*}_{i}$, $r^{*}_{i,0}$ and $s_{i}^{*}$, respectively, where $i$ represents the $i^{\text{th}}$ particle and $i=1\ldots N$. We nondimensionalise via (11) with the only difference being that the mean value of the radii $\bar{r}_{p,0}^{*}$ replaces the length scale $r_{p,0}^{*}$. It has to be noted that this also affects the concentration scale, $\Delta c$, through the initial solubility. Hence in what follows, all dimensionless parameters are the same as those defined in (17), except with $r_{p,0}^{*}$ and $s_{0}^{*}$ replaced by $\bar{r}_{p,0}^{*}$ and $\bar{s}_{0}^{*}$, respectively.

Under the pseudo-steady approximation and assuming that there are no interparticle diffusional interactions, the growth of each particle is now described by an equation of the form (22) with $r_{i}$ and $s_{i,0}$ substituting for $r_{p}$ and $s_{0}$, respectively. The bulk concentration equation must account for all particles, that is

[TABLE]

where $N$ may decrease with time due to Ostwald ripening. Assuming that the solution is sufficiently dilute that $4\pi N_{0}/3\sum\limits_{i=1}^{N}r_{i}^{*3}\ll 1$ and also $V_{M}=M_{p}/\rho_{p}$, we obtain

[TABLE]

In dimensionless form the problem is then governed by the system of differential equations

[TABLE]

for each $i=1,\ldots,N$. Equation (27) represents a system of $N$ non-linear ODEs which must be solved numerically.

4 Comparison of model with experiment

The accuracy of the various forms of the mathematical model will now be ascertained through comparison with the experiments on CdSe nanocrystal synthesis reported by Peng et. al. [23]. Certain parameter values concerning the experiment and CdSe are provided in that paper, others, such as $D,k,s_{\infty}^{*},N_{0}$ must be inferred. Here they will be determined through fitting to equation (23). Since this only contains one free parameter the fitting is a very simple process. We then show that the pseudo-steady state model (PSS model) gives virtually identical results to equation (23). Once it has been established that the analytical solution and the PSS model give such good correspondence and match to experimental data we move on to comparing with the numerical results. The PSS is an approximation to the full system defined by equations (12-16), the analytical solution equation (23) is an approximation to the PSS. Using the full system in the $N$ particle model would be extremely computationally expensive, for this reason the PSS model is the basic component of the $N$ particle model. Consequently we demonstrate that the PSS closely matches the numerical solution of the full model (and consequently so does equation (23)). The $N$ particle model is then examined. First, by setting $N=2$ we are able to demonstrate Ostwald ripening. We go on to show that as $N$ increases the prediction for the average radius tends to the analytical solution for a single particle as $N$ becomes large.

4.1 Parameter estimation via the analytical solution

In Figure 2 we show the first eleven data points from [23]. As discussed earlier not all data points correspond to the pseudo-steady regime, here it is clear that the first three points follow a linear trend so these will be neglected. In the experiment extra monomer was added after three hours, so we have ignored all data beyond the eleventh point. Using the remaining eight data points in the nonlinear least-squares Matlab solver lsqcurvefit to fit to equation (23) we obtain $G\approx 958$. The results of equation (23), with $G=958$, is shown as the solid line in Figure 2.

To determine the necessary parameters for the other models we first note that the maximum radius attained during this part of the experiment is $r_{m}^{*}\approx 3.8$nm$=a/b=D/k$. The experimental concentration at the end of the growth process is known, this defines $c_{eq}^{*}=s_{\infty}^{*}e^{\alpha/r_{m}^{*}}$. This is enough information to determine $D,k,s_{\infty}^{*},N_{0}$. The values taken from [23] are shown as the first ten rows of Table 1, the final four (in italics) are the ones calculated after $G$ has been determined. These values may be used for the pseudo-steady state (PSS) model, or the arbitrary $N$ model. The dashed line in Figure 2 represents the result predicted by the PSS model. Clearly there is excellent agreement between this and the analytical solution, thus verifying the claim that the solubility may be set to a constant without greatly affecting the solution (provided the early time data is neglected).

4.2 Validating the pseudo-steady state approximation

The PSS model is described by equations (19)-(22). Since this forms the basis for the $N$ particle model it is important to verify its accuracy. We do this by comparison with the numerical solution of the full system (12)-(15) (referred to as the full model). Although we have already shown that the PSS is very well approximated by the analytical solution and that for a single crystal the solubility variation may be neglected, we must employ the PSS in the $N$ particle model. This is because when an individual particle’s solubility drops below the bulk concentration Ostwald ripening occurs. The analytical solution neglects variation in solubility, so cannot capture this behaviour.

Problems similar to the full model frequently occur in studies of phase change where it is termed the one-phase Stefan problem (one-phase because the temperature is neglected in one of the phases, this is analogous to neglecting the concentration in the crystal). Examples of one-phase problems occur in laser melting and ablation, Leidenfrost evaporation of a droplet and in supercooled materials. At the nanoscale there are many studies on nanoparticle melting and growth, see [8, 7, 20, 21]. The nanoparticle studies are particularly relevant, since they deal with a spherical geometry and at the nanoscale the melt temperature varies in a manner similar to the variation of the solubility in the current problem. For this reason we follow the numerical scheme outlined in these studies. It requires a standard boundary immobilization transformation and then a semi-implicit finite difference scheme is applied to the resulting equations. For further details see [8, 7]. The PSS model requires the solution of a single nonlinear ordinary differential equation, (22). To do this we simply use the Matlab ODE solver ode15s. Once $r_{p}$ is determined the concentration is given by equation (21).

In Figure 3 we compare the numerical solution of the full and PSS models using the parameters of Table 1, where panel (a) shows the evolution of the particle radius and panel (b) the concentration profile at five different times. In both cases the agreement between the full and PSS models is excellent, thus justifying the use of the simpler PSS model in the $N$ particle system. Panel (a) shows how the particle grows rapidly until around $t\approx 1$ hr when the growth rate decreases, subsequently the radius slowly approaches the maximum value of $r_{p}\simeq 3.8$ nm. This behaviour can be understood by analysing the concentration profiles presented in panel (b). The growth rate is proportional to the concentration gradient, see equation (10). From Figure 3 (b) it is clear that the concentration gradient between the particle surface and the far field is relatively large at small times, leading to rapid growth. After $t\approx 1$ hr the concentration profile is practically flat, leading to a slow growth rate.

4.3 Ostwald ripening with $N=2$

Ostwald ripening (OR) occurs when the bulk concentration falls below a given particle’s solubility. With a single particle the growth rate tends to zero as the solubility and bulk concentration become similar, hence OR will never occur. However, with a group of particles of different sizes theoretically OR must always occur. In practice this could take a very long time and so be difficult to observe. To demonstrate that the current model can predict OR we now investigate the simplest possible case, with two particles.

The system is defined by equation (27) with $N=2$. We take parameter values from Table 1 and choose initial radii $2\text{\,}\mathrm{n}\mathrm{m}$ and $2.5\text{\,}\mathrm{n}\mathrm{m}$. The governing equations may be solved again using the Matlab ODE solver ode15s. Results are presented in Figure 4. The first figure shows the evolution of the radii for more than 25 hours. The solid line represents the evolution of the $2.5\text{\,}\mathrm{n}\mathrm{m}$ particle, the dashed line is the $2\text{\,}\mathrm{n}\mathrm{m}$ one. As can be seen, for small times both particles grow rapidly however, after around 1.7 hours the smaller particle starts to shrink, while the larger one grows linearly. In Figure 4(b) the variation of the particle solubility and concentration is shown, solid and dashed lines correspond to the 2.5, 2nm particle’s solubility respectively, while the dotted line is the bulk concentration. With reference to the variation of the radius it is clear that the rapid growth phase corresponds to a sharp decrease in the bulk concentration. Initially the solubility of each particle is below the bulk concentration and decreases as $r_{p}$ increases. Ostwald ripening begins when the solubility of the smaller particle crosses the $c_{b}$ curve, $c_{b}=s_{1}$ at $t\simeq 1.7$ hr, subsequently its size decreases. The solubility of the larger particle keeps slowly decreasing, in keeping with its slow growth, and remains below the bulk concentration until the end of the simulation. If we continued the simulation the smaller particle would eventually disappear.

4.4 $N$ particles system

To simulate the experiments of [23] we consider a distribution of $N$ nanoparticles where the initial distribution is generated by random numbers, with an initial mean radius $\bar{r}^{*}_{p,0}$ of $2.92\text{\,}\mathrm{n}\mathrm{m}$ and a standard deviation of $\sigma_{o}=8.9$%. In the numerical solution if a particle decreases below 2nm it is assumed to break up and all the monomer enters back into the bulk concentration.

In Figure 5 we compare the prediction for the average radius of 10 and 1000 CdSe particles (dashed lines) with the corresponding data of Peng et al. [23]. The single particle analytical solution for $r(t)$, equation (23), is shown as a solid line. The inset shows the difference between the $N$ particle and analytical solution. In Figure 5(a) the maximum difference between the two solutions is of the order 1.5%, which decreases rapidly with time. The solution with $N=1000$, shown in Figure 5(b), has a maximum difference of the order 0.15% from the analytical solution.

In practice $N$ would be much higher. In Table 1 the population density is given as $N_{0}=8.04\times 10^{21}$ crystals/m3, so in a volume $V\approx 7\times 10^{-6}$m3 we would expect around $10^{16}$ crystals. The Figure demonstrates that as $N$ increases the solution tends to the analytical solution. Given that $N$ is typically very high it is then clearly not necessary to solve the large system. Obviously the analytical solution is easier to understand and implement than a $10^{16}$ particle model. However, it is important to note that in the present example there is no significant Ostwald ripening. From Figure LABEL:twoparticlesradius_colors.eps we observe defocussing starts around 1.7 hours and after nearly 30 hours the radius of the smaller particle has only decreased by 7%. In the experimental data used here extra monomer is added to the solution after three hours and we stop our calculations then. So, in the absence of significant OR we may assume that the analytical solution may be used to predict the average evolution of nanocrystal growth. If OR is to be modelled an $N$ particle model should be used, since this accounts for the solubility of each particle.

5 Conclusions

We have developed a model for the growth of a system of $N$ particles, where $N$ may be arbitrarily large. The model involves a system of first order nonlinear ordinary differential equations, which are easily solved using standard methods. The basis of the $N$ particle model is the pseudo-steady approximation presented in [19]. This incorporates the particle solubility variation which then permits the model to capture Ostwald ripening.

It has been shown that the pseudo-steady model for a single particle has an accurate approximate explicit solution. This was verified by comparison with the full pseudo-steady model. The explicit solution shows that there is a single main parameter controlling nanocrystal growth.

The main drawback to the single particle model is that it cannot capture Ostwald ripening, whereby smaller particles disappear at the expense of larger ones. By studying the system with $N=2$ we were able to emulate Ostwald ripening on a very simple system. By allowing $N$ to become large and calculating the average particle radius we showed that the results approached the single particle explicit solution, which may thus be considered to represent the average growth of a large distribution of particles. A consequence of this is that the $N=2$ model can equally well represent the average radii for an initially bimodal distribution of nanocrystals. An $N>2$ model can represent a much larger distribution of particles.

The main advantage of the current method is that since the single particle model may be solved analytically, and this accurately describes the average radius of a distribution, then the controlling parameters are apparent. This allows us to adjust them and so optimise the growth process, paving the way for efficient large scale production. Since we only have to deal with a single particle the numerical solution is rapid (almost instantaneous), as opposed to previous large scale, time-consuming calculations.

6 Acknowledgements

Claudia Fanelli acknowledges financial support from the Spanish Ministry of Economy and Competitiveness, through the “Maria de Maeztu” Programme for Units of Excellence in R&D (MDM-2014-0445). Tim G. Myers acknowledges the support of Ministerio de Ciencia e Innovacion Grant No. MTM2017-82317-P. Francesc Font and Vincent Cregan acknowledge financial support from the Obra Social “la Caixa” through the programme Recerca en Matemàtica Col$\cdot$laborativa. The authors have been partially funded by the CERCA Programme of the Generalitat de Catalunya.