We theoretically investigate homogeneous crystal nucleation in a solution containing a solute and a volatile solvent. The solvent evaporates from the solution, thereby continuously increasing the concentration of the solute. We view it as an idealized model for the far-out-of-equilibrium conditions present during the liquid-state manufacturing of organic electronic devices. Our model is based on classical nucleation theory, taking the solvent to be a source of the transient conditions in which the solute drops out of the solution. Other than that, the solvent is not directly involved in the nucleation process itself. We approximately solve the kinetic master equations using a combination of Laplace transforms and singular perturbation theory, providing an analytical expression for the nucleation flux. Our results predict that (i) the nucleation flux lags slightly behind a commonly used quasi-steady-state approximation. This effect is governed by two counteracting effects originating from solvent evaporation: while a faster evaporation rate results in an increasingly larger influence of the lag time on the nucleation flux, this lag time itself is found to decrease with increasing evaporation rate. Moreover, we find that (ii) the nucleation flux and the quasi-steady-state nucleation flux are never identical, except trivially in the stationary limit, and (iii) the initial induction period of the nucleation flux, which we characterize as a generalized induction time, decreases weakly with the evaporation rate. This indicates that the relevant time scale for nucleation also decreases with an increasing evaporation rate. Our analytical theory compares favorably with results from a numerical evaluation of the governing kinetic equations.

Organic electronic devices have received a significant amount of interest over the past decades, in large part due to their solution processability, which enables low-cost and large-scale manufacturing.1 These devices are promising candidates for a wide range of applications, including (but not limited to) batteries, transistors, and organic photovoltaics (OPV).2–4 The (dry) active layer in OPV devices, for instance, typically consists of a blend of a polymeric electron donor and a small molecular acceptor with a very complex morphology.5,6 The blend film is prepared by depositing a solution containing these species onto a substrate, which consequently dries due to solvent evaporation that, e.g., in spin-coating processing, is enhanced.6 The complex (phase-separated) morphology emerges spontaneously during this drying process. The optimal design of organic electronic devices is hampered by an incomplete understanding of the impact of the processing conditions on the emergent morphologies.7,8

Different morphologies emerge depending on the type of demixing and the processing conditions, resulting in different functionalities.1 The interplay between the processing conditions and the emergent morphology has been studied extensively in the case of liquid–liquid phase separation, not only experimentally9–11 but also by means of molecular simulations9,12 and phase-field simulations.13–15 If the demixing involves the formation of a crystalline solid phase in a liquid host, the understanding of the effect of the processing conditions on the morphology is much less complete.7,16 Since the level of crystallinity and the grain morphology of a semiconducting thin film govern charge carrier mobility and modulate the trap density, the lack of insight arguably slows down material development and process optimization.17–19 This holds irrespective of the material type, that is, irrespective of whether it is organic, inorganic, or hybrid, and whether the film comprises a single component or a blend.17,20 Hence, there is a need for predictive models that provide not only generic7,8 but also quantitative descriptions of the influence of evaporation on crystal nucleation and growth.

Accurate modeling of crystallization requires a proper understanding of the nucleation of crystallites under the out-of-equilibrium conditions associated with the casting of thin films from an evaporating solvent. Our understanding of nucleation phenomena is largely based on the phenomenological classical nucleation theory (CNT), presuming invariant ambient conditions and a constant concentration of free solute particles. CNT posits that nucleation is a balance between the gain in free energy for transitioning from a metastable to a stable phase and the free energy cost for the required formation of an interface between these two phases.21 The key result of CNT is that the stationary nucleation flux or rate, i.e., the number of stably growing nuclei that emerge per unit volume per unit time, is proportional to the Boltzmann factor associated with the maximum of the nucleation free energy.21 

Clearly, a stationary CNT description cannot be valid during the far-out-of-equilibrium solution deposition of organic thin films, as the processing precludes the emergence of stationary conditions. To shed light on how such processing conditions might influence the nucleation of crystallites in organic electronic devices, we study the effect of solvent evaporation in the context of CNT, making use of analytical theory and numerical methods. For clarity and conciseness, we focus on a liquid solution that contains a single solute and a single solvent that evaporates out of the system, where we treat the evaporation rate as a free parameter and only the solute is explicitly involved in the (unary) nucleation process.

The current understanding of nucleation under non-stationary conditions is almost exclusively based on the so-called quasi-steady-state approximation, which presumes that the stationary CNT description remains valid if we replace all relevant physical quantities by their current-time values.21 However, this kind of ad-hoc approach represents an uncontrolled approximation because the effect of solvent evaporation is not explicitly taken into account but dealt with a posteriori. Notable exceptions are, among others, the works of Shneidman22,23 and those of Trinkaus and Yoo,24 which use a variety of techniques to approximately solve the nucleation master equation while explicitly including the non-stationarity of the problem at hand. These studies have focused on a variety of experimentally relevant origins of non-stationarity, such as the depletion of monomers,25 finite-rate temperature quenches,23 and periodically modulated external pressures in the context of the nucleation of cavities in a liquid.22 It transpires that non-stationary conditions may have a significant effect on the nucleation flux because their impact is not captured by the quasi-steady-state approximation, even if the non-stationarity is not very pronounced.

In the present contribution, we study the time-dependent nucleation flux, defined as the number of supercritical nuclei that appear per unit time per unit volume in a solution that becomes increasingly concentrated due to solvent evaporation. The mathematical description for our model can be mapped onto that of Shneidman,26 for which an approximate analytical solution applies if the nucleation barrier changes slowly or (nearly) linearly in time.26 Our approximate expression for the nucleation flux compares favorably with a numerical solution of the governing kinetic equations that we also pursued. We find that the nucleation flux can be expressed in terms of a (renormalized) quasi-steady-state nucleation flux, evaluated at some shifted time, indicating that the quasi-steady-state flux remains one of the central quantities, in agreement with the description by Shneidman and others.26 The nucleation flux is initially vanishingly small and reaches a (renormalized) quasi-steady state only after some time. We quantify the associated lag time by introducing a generalized induction time, which is a generalization of the time lag to arbitrary cluster sizes that acts as a measure for the time required for a nucleus to grow from a single monomer to a size beyond that of the critical nucleus. The lag time decreases non-trivially with an increasing evaporation rate, in agreement with our numerical calculations. Our analytical results highlight the complexity of non-stationary nucleation in light of the large parameter space and the multitude of time scales that emerge. This makes it far from trivial to extract the relevant physics from numerical calculations alone.

The remainder of this paper is structured as follows. In Sec. II, we summarize classical nucleation theory for simple solutions and extend it to the case where the solvent evaporates out of the solution in Sec. III. Here, we discuss the set of equations that govern the nucleation kinetics and introduce the simplifications we require to solve them. In Sec. IV, we discuss the method of solution, the resulting nucleation flux, and its implications. Next, in Sec. V, we introduce a generalized induction time and study how the relevant time scale of the nucleation process is influenced by the non-stationary conditions. We compare our analytical calculations with a numerical integration of the relevant kinetic master equation in Sec. VI, showing a good overall agreement. Finally, in Sec. VII, we discuss our results and conclude our work. We include a list of symbols in the Nomenclature.

We consider isothermal homogeneous nucleation in a thin film consisting of a fluid containing a single solute in a volatile solvent that evaporates out of the solution. See Fig. 1 for a schematic image of the model. Experimentally, maintaining isothermal conditions requires some form of temperature control since the temperature of the thin film decreases due to solvent evaporation and the associated (latent) heat of vaporization that is required. This can be achieved, e.g., by using a substrate with a large heat capacity and a solute that conducts heat sufficiently fast or by actively maintaining the temperature of the substrate onto which the solution is cast.27 We presume that the solvent acts as a so-called carrier fluid, a background fluid in which only the solute component is involved in the (unary) nucleation process.21 In our model description, solvent evaporation drives nucleation as it continuously increases the local solute concentration and, with that, the degree of supersaturation. We also presume that the evaporation is sufficiently slow so that we can ignore any spatial inhomogeneities along the normal of the fluid film originating from the accumulation of solute at the evaporating surface.28 If this homogeneity approximation is not warranted, we could, in principle, take the spatial dependence of the supersaturation into account.29,30 We stress, however, that in this case, homogeneous nucleation is likely to be preempted by heterogeneous nucleation at the solution–air interface.21 We return to this issue at the end of this paper.

FIG. 1.

Schematic representation of our model. The fluid solution has a thin-film geometry, and the solvent evaporates at the top liquid–gas interface, as indicated by the wavey arrows. Solute molecules are schematically represented by spherical particles, which can aggregate into larger clusters consisting of multiple solute molecules.

FIG. 1.

Schematic representation of our model. The fluid solution has a thin-film geometry, and the solvent evaporates at the top liquid–gas interface, as indicated by the wavey arrows. Solute molecules are schematically represented by spherical particles, which can aggregate into larger clusters consisting of multiple solute molecules.

Close modal

Following the standard description of CNT, we identify the nuclei of the emerging solid phase using the number of solute particles n in a cluster as the sole reaction coordinate. Clusters are treated as non-interacting entities, and we neglect any precipitation-induced fluid flow, which is justified during the early stages of nucleation.21 In our analysis, we treat n as a continuous variable and insist that the clusters evolve via the attachment and detachment of single solute particles (monomers) only. As usual, we presume these stochastic processes to be Markovian, even though recent work suggests that in reality this might not necessarily be true.31 We choose to ignore this for two reasons. First, a non-Markovian description requires explicit knowledge of the so-called memory kernels,31 which, as far as we are aware, are presently unknown for our model system, even in the absence of evaporation. Second, if solvent evaporation has an effect on nucleation in solution, we expect it to already present itself in CNT. Hence, for us, a CNT approach provides a convenient way to evaluate the role of externally imposed non-stationarity in the cluster size distribution. As we shall see, even within the framework of CNT, our calculations turn out to not be straightforward.

We expect the time evolution of the dimensionless cluster number densities ρ(n, t) for n > 1 to be described by the continuity equation in n-space,21,32
ρ(n,t)t=J(n,t)n+K(n,t),
(1)
where J(n, t) is the cluster flux density in n-space, i.e., the number of nuclei of size n that emerge per unit of time at time t, and K(n, t) is a source term expressed per unit of time that enters due to solvent evaporation. We render the cluster densities ρ(n, t) dimensionless via implicit multiplication with the molecular volume of a solute molecule v0. (A list of symbols can be found in Nomenclature.) The time evolution of the monomer density ρ(n = 1, t) is included as a boundary condition, which we introduce at a later stage.
As is well known, the cluster flux density J(n, t) obeys21,32
J(n,t)=k+(n,t)Feq(n,t)nρ(n,t)Feq(n,t),
(2)
with
Feq(n,t)=ρ1(t)eβΔW(n,t),
(3)
the instantaneous “equilibrium” cluster density, k+(n, t) is the (microscopic) attachment frequency for a monomer to an n-mer, which depends on the (time-dependent) monomer concentration and may depend on the cluster size via, e.g., the surface area for reaction- or interfacial-limited growth or the linear dimensions of a cluster for diffusion-limited growth.21 In Eq. (3), β denotes, as usual, the reciprocal thermal energy 1/kBT and ΔW(n, t) and the thermodynamic work required to form an n-mer, that is, a cluster consisting of n solute particles. For the remainder of this work, we set β = 1. If K(n, t) = 0, the set of equations given by Eqs. (1) and (2) is commonly known as the Zeldovich equation.33,34
Within the framework of CNT, the dimensionless work ΔW(n, t) to form an n-mer from n monomers reads35,
ΔW(n,t)=(n1)Δμ(t)+γA0(n2/31),
(4)
with Δμ(t) a difference in chemical potential between the condensed and dissolved states, γ the interfacial tension of a flat interface, and A0 a shape factor with units of area. The first term encodes that in an over-saturated solution, the lowest free energy state corresponds to the condensed, solid state, not the metastable dissolved state. The free energy difference per monomer between these states is given by the difference of their chemical potentials, Δμ.21 As is usually performed, for this, we take as a reference a dilute solution,
Δμ(t)=lnρ(1,t)ρ=lnS(t),
(5)
which is generally reasonable for the production processes for the thin-film organic electronics we study.21 Here, ρ is the solubility limit and S(t) ≡ ρ(1, t)/ρ is the degree of supersaturation. The condensed phase is more stable than the dissolved phase if S(t) > 1. We neglect the effect of the Laplace pressure in Eq. (5) on account of the (near) incompressibility of the solid, condensed phase.21 

The second term in Eq. (4) describes the cost of forming an interface between the solid and fluid phases. We presume that a well-defined and sharp interface exists between the two phases, also known as the capillarity approximation. Note that we tacitly assume the interfacial tension does not depend strongly on the crystal face exposed to the solvent, so we need not take the crystal structure into account. The subtraction by minus one in both contributions in Eq. (4) ensures that the free energy to form a monomer is zero, providing a thermodynamically consistent description. Although this monomer correction has been subject to criticism (see, e.g., Refs. 21 and 35), we argue that we require it under conditions of solvent evaporation to ensure that the number of solute monomers remains conserved, as we show below.

The most relevant properties of the free energy barrier ΔW are35 (i) the maximum height
ΔW*(t)=3ε213n*2/3+2n*1,
(6)
(ii) the critical nucleus size associated with this maximum height
n*(t)=2σ3lnS(t)3,
(7)
and (iii) the so-called Zeldovich factor
Z(t)12π2ΔWn2n=n*=1π1εn*.
(8)
The Zeldovich factor is a measure of the width of the nucleation barrier.21 All three properties depend on time on account of the ongoing process of the evaporation of the solvent. Here, we define ε3n*1/3σ1/2 for reasons of convenience, with σγA0 a dimensionless interfacial tension, recalling that A0 is a geometric constant. For future reference, we notice that the parameter ɛ can be seen as a measure for the (square root of) the reciprocal barrier height.
Our final model ingredient is the evaporation source term K(n, t). In our idealized model, we assume that solvent evaporation and nucleation are independent processes. The former indirectly affects the latter, as solvent evaporation steadily increases the concentrations of the clusters. We can obtain the evaporation source term by a thought experiment in which the number of clusters of a given size N(n) is fixed in time. Since solvent evaporation only affects the total volume V(t) of the solution, we must have ρ(n, t) = v0N(n)/V(t), with v0 the volume of a solute molecule. This is identical to neglecting nucleation, so setting J = 0. The evaporation source term K(n, t) can now be expressed in terms of the volume V(t), where at time t = 0, the evaporation starts and the film has the volume V(0), because then ∂ρ(n, t)/∂t = K(n, t). This results in
K(n,t)=ρ(n,t)χ(t)tχ1(t)=ρ(n,t)tlnχ(t),
(9)
where, for convenience, we define the reciprocal dimensionless volume of the film χ(t) as
χ(t)=V(0)V(t).
(10)
To make the functional dependence of χ(t) on time t explicit, we would, in principle, need to set up a model for solvent evaporation. We leave it unspecified for now and introduce our solvent evaporation model at a later stage. We shall presume that Eq. (9) holds even if nucleation takes place and the number of clusters of a certain size is time dependent or, equivalently, that nucleation and evaporation are independent processes. This is reasonable during the early stages of nucleation, where the number density of clusters remains low and solute depletion is negligible.

In an experimental setting, the solution arguably starts out in an under-saturated state with Δμ(0) < 0, which is then gradually quenched via solvent evaporation into a supersaturated state where Δμ(0) > 0. In this work, we consider a solution that is supersaturated at time zero when the evaporation starts. This is justified because nucleation barriers in excess of a few hundred times the thermal energy kBT, relevant just beyond the binodal, cause nucleation to be an exceedingly improbable and slow process. This is even more pertinent if we consider the relatively small volumes relevant for thin films, on the order of milli- to microliters,36,37 resulting in negligible nucleation.21 A critical supersaturation exists below which the expected number of nucleation events in a thin film on the time scale of the experiments is (much) lower than unity.38,39 Hence, we can justifiably set the initial time t = 0 at or slightly below this critical supersaturation. This critical supersaturation can be estimated as the supersaturation for which a single nucleus emerges on the experimental time scale, i.e., J(t) × V(t) × texp = 1. Here, J is the nucleation flux, implicitly dependent on the supersaturation S(t), V(t) is the volume of the solution, and texp is the experimental time scale. The critical supersaturation can, therefore, only be estimated if the kinetic and thermodynamic properties of the nucleating solute are known, which, as far as we are aware, are not known for the compounds relevant to us. Therefore, for the sake of simplicity, we shall presume that the initial degree of supersaturation is a free parameter.

If the degree of supersaturation is stationary, S(t) = S, which is true in the absence of monomer depletion and solvent evaporation, implying that ρ(1, t) = ρ > ρ, K(n, t) = 0, and ∂ρ(n, t)/∂t = 0, the steady-state solution of Eqs. (1) and (2) reads21 
JSSρZk*+eΔW*=Zk*+Feq(n*),
(11)
at least if the nucleation barrier is sufficiently large, ΔW* ≫ 1. Here, k*+ denotes the attachment frequency of a monomer to the critical nucleus, Z is the earlier introduced Zeldovich factor that obeys Eq. (8), and Feq(n) is the hypothetical equilibrium distribution given by Eq. (3).

Under transient conditions and provided solvent evaporation is sufficiently slow on the time scale of nucleation, the so-called Quasi-Steady-State (QSS) approximation based on Eq. (11) yields a reasonably accurate prediction for the nucleation flux. This corresponds to letting all quantities in Eq. (11) depend on time. Since this does not actually account for how precisely solvent evaporation influences the nucleation rate, we focus next on solving Eqs. (1) and (2) explicitly while allowing for non-stationary effects. This means that we will go beyond the QSS approximation. Nevertheless, we shall see that the nucleation flux obtained from the QSS approximation remains a central quantity in our results.

In this section, we first rewrite our governing equations into a form that is better suited for analytical treatment and introduce the “non-stationarity index,” which is proportional to the ratio of the nucleation and evaporation time scales, as the key control parameter throughout our work. Second, we discuss the approximations we require to analytically, albeit approximately, solve the nucleation equation, and finally introduce the evaporation and aggregation models we use throughout the remainder of this article.

In order to investigate how solvent evaporation affects nucleation, we rewrite our governing equations into a form better suited for analytical treatment. To this end, we introduce the normalized density of clusters of size n at time t as
ν(n,t)ρ(n,t)Feq(n,t),
(12)
where Feq(n, t) is given by Eq. (3). As such, ν(n, t) expresses deviations from the instantaneous “equilibrium” cluster size distribution. We obtain the time-evolution equation for the normalized cluster density by combining Eqs. (1)(4), (9), and (12), reading
ν(n,t)t+ν(n,t)tlnχ1(t)Feq(n,t)=Feq1(n,t)nk+(n,t)Feq(n,t)ν(n,t)n,
(13)
with χ(t) defined by Eq. (10). Invoking the no-depletion boundary condition that conserves the number of solute molecules (“monomers”) results in a boundary condition for the monomer density ρ(1, t) = ρ(1, 0)χ(t) that depends on time on account of solvent evaporation, marking the deviation from standard CNT. This in turn allows us to simplify the second term on the left-hand side of Eq. (13) by defining the quantity
ζ(n,t)tlnχ1(t)Feq(n,t)=tΔW(n,t)=n1tΔμ(t),
(14)
where we make use of the fact that the nucleation-free energy Eq. (4) depends on time only through the chemical potential difference of the monomers in the two phases, Δμ(t).
For the sake of clarity, we rewrite Eq. (13) in the form of a diffusion-advection equation or a Fokker–Planck equation if we interpret ν(n, t) as a scaled probability distribution function,
ν(n,t)t+ν(n,t)ζ(n,t)=nk+(n,t)ν(n,t)n+g(n,t)ν(n,t)n,
(15)
where
g(n,t)k+(n,t)ΔW(n,t)n.
(16)
Equation (15) shows that cluster growth has a diffusive contribution, which according to the standard interpretation of Fokker–Planck equations represents how at the microscopic level thermal fluctuations affect cluster growth, and an advective contribution that is driven by the generalized force −ΔW(n, t)/∂n.40 We interpret this advective contribution as the deterministic part of the cluster growth, as is usual in CNT.32,34 The shape of the free energy landscape in Eq. (4) enforces g(n, t) to be negative for n < n*, implying shrinking clusters, and positive for n > n*, implying growing clusters, consistent with our interpretation of the deterministic part of Eq. (15).
As it turns out, the (n, t) coordinate system is not well-suited to study non-stationary nucleation. The reason is that the non-stationary conditions induced by solvent evaporation must result in a time-dependent nucleation time scale tcrit(t), which we introduce below, and a characteristic “length” scale, given by the size of the critical nucleus n*(t). Since both drift with time, we switch to a more natural reaction coordinate space.26 First, we introduce the reduced size coordinate x = n/n*(t), with n*(t) the critical cluster size that depends on the instantaneous supersaturation S(t) according to Eq. (7). The derivatives in Eq. (15) transform to /nn*1/x and to /tn/tx+(x/t)n/x, where the subscript indicates which variable is to be held constant while taking the time derivative. This implies that we introduce a co-moving coordinate system with x and t as independent variables. Second, in order to find the appropriate transformation for the time coordinate, we must first obtain the appropriate nucleation time scale, tcrit(t). To do so, we introduce in Eq. (15) the reduced attachment frequencies k̄+(x,t)=k+(x,t)/k+(1,t) with k+(1,t)k+*(t), the critical attachment frequency evaluated at the critical cluster size n = n* (or x = 1), and make use of Eq. (8) to find
tcrit(t)ν(x,t)t+tcrit(t)ν(x,t)ζ(x,t)=x12k̄+ε2ν(x,t)x+g(x,t)xν(x,t).
(17)
Here, we redefine the deterministic growth rate as
g(x,t)12k̄+ε2ΔWx+xtcrit(t)lnn*(t)t,
(18)
with the nucleation free energy barrier in the (x, t)-coordinate system
ΔW(x,t)=3ε23x2/32x3n*2/3(t)2n*1(t).
(19)
The second square-bracketed term on the right-hand side in Eq. (19) is due to the monomer correction introduced earlier. Recall that the parameter ɛ is a measure for the reciprocal of the nucleation barrier height [see Eq. (6)], which corresponds to x = 1 in Eq. (19).
In Eq. (17), we can identify the nucleation time scale tcrit(t) as21,26,41,42
tcrit(t)=2π1k+(n*,t)1Z(t)2,
(20)
where Z(t) is the Zeldovich factor expressed in Eq. (8). We recall that the Zeldovich factor is a measure of the reciprocal width of the critical region. The nucleation time scale represents the lifetime of a critical nucleus. It expresses the notion that the growth of near-critical clusters is governed by a random walk in n space over a distance Z−1(t) due to thermal motion. Under stationary conditions, the nucleation time scale can also be derived straightforwardly from Eq. (15), noting that g(x) is negligibly small in the critical region. The factor of two accounts for the fact that only half of the clusters grow to a supercritical size, and the factor π−1 is a common, although non-standard, factor that we include for convenience.41,43
Next, we introduce a reduced time τ, making use of the lifetime of a critical nucleus tcrit(t) in Eq. (20). From Eq. (17), we observe that the derivatives with respect to time, one of which is implicit in ζ(x, t), are multiplied by the lifetime of a critical nucleus, tcrit(t). Hence, we conclude that the convenient time variable would be one that yields dτ = dt/tcrit(t) because it accounts for the difference between the rate at which time passes for the observer and the time experienced by a cluster.44 The reduced or integrated time τ reads44 
τ=0tdttcrit(t),
(21)
expresses the time that passes in units of the lifetime of the critical nucleus tcrit and simply takes into account that the unit of time drifts with time, emphasizing the highly non-linear character of the problem at hand.
As the final step in our effort to rewrite the nucleation equation in a useful form, we focus on the dominant contribution of the non-stationary solution conditions in Eq. (17) due to solvent evaporation, ζ(x, τ), and recall that this quantity is defined as the (negative) time derivative of the nucleation barrier. Using the definition of this quantity in Eq. (14), we may now write it in terms of the (x, τ)-coordinate system as
ζ(x,τ)=τΔW=n*τlnχ(τ)xn*1ω(τ)ψ(x,τ),
(22)
with ω(τ) = n*(τ) ln χ(τ)/∂τ a “non-stationarity coefficient” or “non-stationarity index”26 that measures how strongly the non-stationary conditions affect the nucleation of a critical nucleus of the solid phase of the solute. This parameter we find to play a key role in our work. The non-stationarity index (or coefficient) depends on time but, as we shall argue below, only weakly if the evaporation is sufficiently slow. The function ψ(x,τ)=(xn*1(τ)) is the cluster size-dependent part of ζ(x, τ) and follows straightforwardly from Eq. (14).
The physical interpretation of the non-stationarity coefficient is straightforward if we express it in terms of the evaporation time scale. This time scale is quantified by the relative rate of change of the monomer concentration due to solvent evaporation, ρ(1/n*,τ)1ρ(1/n*,τ)/τ=τlnχ(τ). See also Eqs. (1) and (9). The resulting time scale of solvent evaporation is
tevap(τ)=lnχ(τ)/τ1.
(23)
Hence, using Eqs. (21)(23), we can now express the non-stationarity coefficient in units of the actual time t as
ω(t)n*(t)tcrit(t)/tevap(t).
(24)
Clearly, the quantity ω(t) measures the relative magnitude of the time scales for evaporation and nucleation. Multiplication with the critical nucleus size n*(t) accounts for the fact that the evaporation timescale enters via the time derivative of the chemical potential difference, which is measured per monomer. For the remainder of this paper, unless explicitly stated otherwise, we consider all quantities to be in the (x, τ)-coordinate system.

We have now rewritten the nucleation master Eq. (17), the growth rate g(x) (18), and the rate of change of the nucleation barrier due to solvent evaporation (22) in the (x, τ)-coordinate system. We have performed so from the (n, t)-coordinate space equations without approximations, except for the introduction of the (very customary) no-depletion boundary condition. Since the equations are in the form of a Fokker–Planck equation with size- and time-dependent coefficients, exact analytical solutions are generally not known.45 Before we introduce the approximations that allow us to approximately solve (17), we first point out that our model in the (x, τ)-coordinate frame maps onto the single component model by Shneidman,26 wherein the supersaturation changes due to, e.g., an externally applied time-dependent pressure or temperature, rather than a source. The equivalence of the models stems from the cancellation of the evaporation source term by the division of the cluster densities ρ by the hypothetical equilibrium cluster densities Feq. Hence, our theory is a generalization to mixtures in which only one component is conserved. Because the remainder of this work uses the nucleation master Eq. (17) in the (x, τ)-coordinate frame, we can apply the methods and results of Shneidman to our case.26 

Following Ref. 26, very accurate approximate solutions to Eq. (17) can still be obtained in experimentally relevant limits by presuming that (i) the nucleation barrier is much larger than the thermal energy, so ΔW* ≫ 1, by itself a prerequisite for applying CNT, and (ii) solvent evaporation is sufficiently slow so that the following conditions hold:
131n*n*τ=1εετε2|ω|1,1ωωτ1.
(25)
Under these constraints, the critical cluster size, n*, the reciprocal of the square root of the nucleation barrier, ɛ, and the non-stationarity index, ω, become so weakly time-dependent in the co-moving coordinate system that they are essentially constant. Furthermore, the time derivative of both the free energy barrier and the reduced attachment frequencies is, for all intents and purposes, independent of time so that ζ(x, τ) ≡ ζ(x) and k̄+(x,τ)k̄+(x).46 Considering that a high nucleation barrier implies ɛ≪ 1 [see Eq. (6)], we find from Eq. (25) that such an approach is warranted for the variables n* and ɛ if ω<O(ε2). We justify our treatment of ζ as a constant a posteriori.

With these simplifications, we may neglect for the scaled cluster sizes x(ε2ω)1 the contribution from the second term in Eq. (18), reading x( ln n*/∂τ) ∼ −2ω in reduced reaction coordinates, which would otherwise make analytical treatment exceedingly more difficult. Since the clusters grow essentially unidirectionally for x ≫ 1, going downhill in the free energy landscape, our results for x<(ε2ω)1 are not affected by the omission of this term. We stress that because ɛ2ω is assumed to be very small, this approach is justified even if the scaled cluster size x is not particularly small. However, the omission would result in an overprediction of the growth rate of larger clusters, which are outside the scope of this work as they do not contribute to nucleation but only to subsequent domain growth. In this late stage regime, our nucleation model, which treats clusters as non-interacting entities, becomes less accurate, if only because the likelihood that clusters interact increases with the cluster size. For stationary nucleation, this situation is typically dealt with by applying the Johnson–Mehl–Avrami–Kolmogorov model, where one treats nucleation and growth as separate processes and explicitly corrects for the excluded volume of clusters.47 

Under the assumptions expressed in Eq. (25), that is, taking ɛ, n*, and ω to be constant and ζ(x) and k̄+(x) to be independent of time, Eq. (17) reduces to
ν(x,τ)τ+ν(x,τ)ζ(x)=12ε2xk̄+(x)ν(x,τ)x+g(x)xν(x,τ).
(26)
The deterministic growth rate of the clusters, g(x), is time-independent in the current coordinate system and given by
g(x)=12ε2k̄+(x)ΔWx=3(1x1/3)k̄+(x),
(27)
where we inserted Eq. (19). Due to the cancellation of the factor ɛ2 in Eq. (27) by multiplication with ɛ−2 in Eq. (19), g(x) represents the dominant mechanism for cluster growth over the whole size domain, except for sizes near that of the critical cluster. The extension of our central equation Eq. (26) to heterogeneous nucleation or more complex nucleation barriers is relatively straightforward. In the  Appendix, we outline how evaporation modifies heterogeneous nucleation within CNT.
We solve Eq. (26) for an instantaneously quench into a metastable state with an arbitrary degree of supersaturation at τ = 0, which also marks the point when solvent evaporation commences. Subsequently, we invoke our initial and boundary conditions that in the co-moving coordinate system become
ν(x>n*1,0)=0,ν(n*1,τ)=1,limxν(x,τ)=0.
(28)
Here, ν(x, τ) = ρ(x, τ)/Feq(x, τ), with ρ(x, τ) the cluster number densities and Feq(x, τ) the hypothetical equilibrium cluster densities in Eq. (3). We reiterate that Eq. (28) ascertains that (i) only monomers are present in the solution at τ = 0, (ii) the number of free solute molecules remains fixed, and (iii) the reduced cluster size distribution ν(x, τ) approaches zero for large clusters: x. The second condition follows from Eqs. (3) and (12), when expressed in x, τ coordinates, and only holds if the “monomer correction” is applied. The third condition enforces a finite number of solute molecules per volume, or cluster density, noting that the (hypothetical) “equilibrium” system described by Eq. (3) predicts a diverging cluster density in the limit of very large clusters. This, obviously, is unphysical.
Although we can (approximately) solve Eq. (26) without specifying an evaporation model or how monomers attach to a cluster, we introduce both to make our results more tangible. Molecular aggregation may be diffusion-limited if the attachment of a monomer to a cluster is much faster than the diffusion toward it, or reaction-limited if the opposite applies. For the growth of clusters of solution-processed organic molecules, either of these limits or an intermediate regime may be relevant, depending on the processing conditions and material properties.48–50 We opt to invoke a reaction-limited growth model,21 as it allows us to calculate all quantities of interest analytically. This, as it turns out, appears not to be possible for the diffusion-limited growth model.26 In the reaction-limited case, the attachment frequencies are proportional to the surface area of a cluster and can be written as21,51
k+(x,t)=k0n*2/3ρ1(t)x2/3,
(29)
where the dimensionless quantity x2/3=(n/n*)2/3 is proportional to the surface area of a cluster, ρ1(t) is the dimensionless monomer density, and k0 is a constant with units of reciprocal seconds and proportional to the so-called attempt frequency. Note that n* enters Eq. (29) to express k+(x, t) in terms of the reduced cluster size x rather than the number of monomers per cluster n. For simplicity, we neglect that the attempt frequency may depend on the cluster size if the latter affects the energetic cost of forming an activation complex.21 
Next, we introduce a model for solvent evaporation, assuming that the evaporation rate is proportional to the concentration of solvent molecules.27 If the initial solute concentration is low and we limit ourselves to the initial stages of the drying process, this then predicts the volume V(t) of our solution to decrease linearly with time as
V(0)/V(t)χ(t)=11αt.
(30)
We treat the reciprocal of the drying time α, a measure for the rate of evaporation, as a freely adjustable parameter and maintain our assumption that monomer depletion due to aggregation is negligible. This description of the evaporation kinetics applies, e.g., in the context of spin-coating of solutions in the initial stages of drying,10,27 where the drying time can be expressed in terms of the height h of the film α1=hdh/dt1. Experimentally, drying times for micrometer-thick films can be mere seconds.16 
Using the models for monomer attachment Eq. (29) and solvent evaporation Eq. (30), we can also find an analytical expression for the non-stationary index ω, which makes explicit under what conditions our constraints (25) apply. Deriving an expression for ω requires an explicit expression for the scaled time τ(t) as a function of the actual time t. Upon combining Eqs. (20), (21), and (29), we find
τ(t)=1tcrit(t=0)0tdtχ(t)(1+Δμ01lnχ(t))2,
(31)
which, combined with Eq. (30), yields
τ(t)=τs(1+Δμ01lnχ(t))31,
(32)
expressed in terms of the reciprocal volume χ(t) defined by Eq. (30). Here, Δμ0 = ln S(t = 0) is the difference in the chemical potential between the two phases at the initial time, and τs = Δμ0/3αtcrit(0) is a time scale that we interpret below. Unfortunately, the inversion of Eq. (32) does not produce a simple expression of the actual time t in terms of the scaled time τ. For early times, a Taylor expansion yields t/tcrit(0) = τ(1 − /τs), with a positively valued B = (2 + Δμ0)/6, showing that the actual time t increases sub-linearly with the scaled time τ. This is expected because solvent evaporation speeds up nucleation, and the time τ measures the time that passes from the point of view of the nucleating solute particles.
From the above-mentioned analysis, combined with Eq. (22), we find
ω(t)=ω(0)1+Δμ01lnχ5=ω(0)τ(t)τs+15/3,
(33)
where ω(0) = αtcrit(0)n*(0) is evaluated at t = 0, which is consistent with Eq. (24). The timescale τs represents the characteristic time for which ω drifts with time. For high nucleation barriers ΔW* and sufficiently large critical cluster sizes n*, we can express the time scale τs as the ratio of the height of the dimensionless nucleation barrier and the non-stationarity index ω, both evaluated at τ = 0, as
τs=2ΔW*(0)/3ω(0)2/3ε2ω(0)Oε3,
(34)
where we neglect a correction term of order n*2/3, use the expression for the nucleation barrier Eq. (6), and make explicit that ω(0) implicitly depends on the initial barrier height via n* and tcrit as ω(0) ∼ ɛ−5, see Eqs. (7), (20), and (22). Although Eqs. (33) and (34) are specific for the reaction-limited aggregation model, we do not expect that different attachment models will qualitatively change our results, as the relevant physics remains the same.

From Eq. (34), we conclude that the two constraints presented in Eq. (25) can actually be written in terms of a single constraint: τs ≫ 1. Hence, this time scale must play a central role in our work. It also highlights that the non-stationarity index ω, the square root of the barrier height ɛ, and the critical cluster size n* cannot actually be constant in our model system where the solvent evaporates out of our film. Nevertheless, if τs is sufficiently large, treating them as constants is justified, provided ττs. For now, we assume that this condition always holds and return to this constraint and how to account for a weak drift in time in Sec. VI.

In this section, we derive an analytical albeit approximate expression for the nucleation flux J(x, τ). As we shall see, this quantity exhibits two distinct time regimes. The early time regime, which can be seen as a lag or induction time, describes the response of the nucleation flux to the initial quench at τ = 0, after which the late-time regime commences. We shall discuss how the nucleation flux is affected by solvent evaporation at both early and late times. In order to derive an approximate analytical expression for the nucleation flux J(x, τ), we need to solve Eq. (26) for ν(x, τ) first. This we do by first applying a Laplace transform with respect to time and subsequently using the framework of singular perturbation theory in order to find a perturbative solution using ɛ, the square root of the reciprocal barrier height, as the perturbation parameter. We point out that the region around the critical cluster |x − 1| ≪ 1 represents a so-called transition layer of width ɛ, characterized by a rapidly changing ν(x, τ).52–55 Based on standard singular perturbation techniques,55 this allows us to subdivide the size-domain into subcritical cluster x < 1, critical cluster |x − 1| ≪ 1, and supercritical cluster x > 1, and we can solve Eq. (26) separately for each region. We make the solutions in the three separate regions consistent with each other via a procedure known as matching.55 We present these calculations in supplementary material I.

In Fig. 2, we show the solution for the scaled cluster densities ν(x, τ) in late time under stationary conditions in blue and under non-stationary conditions for ω*=(1n*1)ω=101 and ɛ = 10−2. The solution for subcritical clusters is shown in green and for critical clusters for x − 1 ≪ 1 in orange. The supercritical solution for x > 1 is not shown, as we find it to be zero everywhere (see supplementary material I). The vertical dashed lines represent the edges of the critical region. In the right panel, we use the stretched variable X = (x − 1)/ɛ to better visualize the scaled cluster densities ν(x, τ) near x = 1. For stationary nucleation, the scaled cluster density is nearly unity for subcritical clusters and decreases rapidly at the critical cluster size. This highlights that the subcritical cluster densities approach their “equilibrium” value given by Eq. (3). For non-stationary nucleation, on the other hand, we observe that for subcritical cluster densities, ν(x, τ) deviates significantly from unity. The origin of the smaller value for the reduced cluster size is that the aggregation process cannot immediately respond to changes in the ambient conditions. Hence, the cluster densities “lag” behind the equilibrium value. This effect increases in magnitude with an increasing rate of evaporation.

FIG. 2.

Left: The cluster size distribution ν(x, τ) under stationary conditions (ω* = 0, blue),21 and under non-stationary conditions for ω* = 10−1 and ɛ = 10−2, with the solution for x < 1 (green) given in Eq. (S.6) and the inner solution around x = 1 (orange), given by the Laplace inverse of Eq. (S.8). The solution for x > 1 is not shown, as it is zero everywhere. The vertical dashed lines indicate the edges of the critical region. Right: data and colors same as left, shown as a function of the stretched variable X = (x − 1)/ɛ and zoomed in on the critical region. Here, we assume that the attachment frequencies vary very slowly with x in the subcritical (x < 1) region for simplicity, so set k̄+(x)=1.

FIG. 2.

Left: The cluster size distribution ν(x, τ) under stationary conditions (ω* = 0, blue),21 and under non-stationary conditions for ω* = 10−1 and ɛ = 10−2, with the solution for x < 1 (green) given in Eq. (S.6) and the inner solution around x = 1 (orange), given by the Laplace inverse of Eq. (S.8). The solution for x > 1 is not shown, as it is zero everywhere. The vertical dashed lines indicate the edges of the critical region. Right: data and colors same as left, shown as a function of the stretched variable X = (x − 1)/ɛ and zoomed in on the critical region. Here, we assume that the attachment frequencies vary very slowly with x in the subcritical (x < 1) region for simplicity, so set k̄+(x)=1.

Close modal

In Fig. 3, we show the scaled cluster densities ν(x, τ) for ɛ = 10−2 and non-stationarity indexes ω* = 0.01 (top) and ω* = 0.1 (bottom) at different moments in time. In the left panels, we show the subcritical cluster distributions. After the initial quench into the meta-stable region at τ = 0, the cluster densities populate in order. For subcritical clusters, the “growth front” is infinitely sharp, which is not true in reality. This is a consequence of the singular perturbation method we employ to find the approximate solution for the scaled cluster densities. In the right panels, we show the scaled cluster densities for the critical clusters. Using the stretched variable X = (x − 1)/ɛ to magnify the critical region, we observe that the wavefront broadens upon entering the critical region. After a given amount of time, a (quasi)-stationary state is reached, and the scaled cluster densities no longer change. Note that this quasi-stationary state is not identical to the well-known quasi-steady-state approximation of nucleation.21 

FIG. 3.

The cluster size distribution as a function of the cluster size for various sizes for the asymptotic inverse barrier height ɛ = 10−2, the reduced cluster frequencies k̄+(x)=1, and non-stationarity indexes ω* = 0.01 (top) and ω* = 0.1 (bottom). Top left and bottom left: The subcritical cluster size distribution at various moments in time after the initial quench. As the wavefront reaches the critical cluster size, the subcritical size distribution ν(x, τ) no longer changes. After this time, the wavefront enters the inner region, as shown in the right figures.

FIG. 3.

The cluster size distribution as a function of the cluster size for various sizes for the asymptotic inverse barrier height ɛ = 10−2, the reduced cluster frequencies k̄+(x)=1, and non-stationarity indexes ω* = 0.01 (top) and ω* = 0.1 (bottom). Top left and bottom left: The subcritical cluster size distribution at various moments in time after the initial quench. As the wavefront reaches the critical cluster size, the subcritical size distribution ν(x, τ) no longer changes. After this time, the wavefront enters the inner region, as shown in the right figures.

Close modal

Using the approximate solutions for ν(x, τ) that we show in Figs. 2 and 3, we can now find an approximate expression for the nucleation flux J(x, τ) by combining them with Eq. (2). Here, we focus attention on the supercritical clusters, as these are the ones that continue to grow. This is actually far from trivial to derive because our expression for the supercritical cluster densities ν(x > 1, τ) is always zero due to the singular nature of our perturbation scheme, which is not true in reality. To address this problem, we follow Shneidman26 and Shi et al. 54 and introduce a mass-conservation argument to propagate critical clusters to supercritical size. As we show in supplementary material II, the nucleation flux J(x, τ) under non-stationary conditions can actually be expressed as an integral of the nucleation flux under stationary conditions. This shows that stationary and non-stationary nucleation are intrinsically connected. Upon calculating this integral expression, we obtain an analytical expression for the non-stationary nucleation flux. In the remainder of this section, we study the non-stationary nucleation flux J(x, τ) that we derive in supplementary material II.

As shown in supplementary material II, the non-stationary nucleation flux of supercritical clusters with x > 1, in the context of homogeneous unary nucleation in a solution that becomes increasingly concentrated as a result of solvent evaporation, reads
J(x,τ)=JtQSS(x,τ)Qω*+1,e(ττi(x))J0(x),
(35)
with
JtQSS(x,τ)=Γ(ω*+1)JQSS(ττi(x)+τcor),
(36)
a renormalized, time-shifted quasi-steady-state, or tQSS for short, nucleation. Furthermore, Γ(⋯) denotes the usual gamma function, Q(n, z) = Γ(n, z)/Γ(n), a regularized gamma function with Γ(n, z) the upper incomplete gamma function of the variables n and z,56 and
ω*=ω(1n*1),
(37)
with ω the non-stationary index defined in Eq. (24). The quasi-steady state flux JQSS is defined as
JQSS(τ)=2Z1π1Feq(x=1,τ)=JSSeω*τ,
(38)
with JSS the steady-state nucleation flux given in Eq. (11). Equation (35) is, by construction, only valid for x > 1, but this fits, as already pointed out, our purpose.
The timescales that emerge in Eqs. (35) and (36) are the so-called incubation time τi(x) that is a measure for the mean growth time for a single cluster from monomer size to a size x, and the time scale τcor that emerges only in the context of non-stationary nucleation. The former is defined as
τi(x)=P1/n*xdxg(x)+ln6ΔW*2τstat,=x1/3+ln(x1/31)2+ln18ε2,
(39)
where P denotes that we must take the principal value for reasons outlined in supplementary material II, g(x′) is the growth rate defined in Eq. (27), and we define
τstat=1/n*1dx1g(x)1x1ln(1n*1).
(40)
The incubation time also emerges in nucleation under stationary conditions. We use the reaction-limited model Eq. (29) to evaluate the incubation time Eq. (39). It increases with increasing x but does so very weakly (sublinearly). Note that the incubation time diverges for x = 1. This nonphysical behavior is irrelevant, as our expression for the nucleation rate is only valid for x > 1.
The second time scale τcor that emerges is given by
τcor=τstatτnon-stat=1312+On*2/3,
(41)
where we use the definition
τnon-stat=1/n*1dxψ(x)ψ(1)1g(x)1x1ln(1n*1),
(42)
with ψ(x)=xn*1. Again, for the evaluation of Eq. (41), we apply the reaction-limited aggregation model Eq. (29). We ignore the corrections in n*, as they are very small for critical cluster sizes n* ≫ 1. The quantity Eq. (41) is a contribution to the time-shift in Eq. (35) and emerges only under conditions of non-stationary nucleation.

The gamma function that acts as a prefactor in Eq. (36) can be absorbed in the time shift, resulting in the time shift ττi(x)+τcor+ω*1lnω*Γ(ω*), making explicit that the evaporation rate actually influences the delay or lag time. We return to this in Sec. V and show that the time shift due to solvent evaporation also emerges naturally from a (generalized) induction time. The quantity J0(x) in Eq. (35) is equal to the first term in Eq. (35), evaluated at τ = 0, and arises from the condition that the nucleation flux at τ = 0 must be zero. These expressions for the nucleation flux also apply to the case of heterogeneous nucleation, introducing minor modifications only as shown in the  Appendix.

Equation (36) is the first main result of this paper. It applies under the presumptions of (i) a high nucleation barrier, implying that its reciprocal value ɛ2 must be very small, and (ii) that the non-stationarity index ω* is constant in the co-moving reference frame and satisfies ω*ɛ2 ≪ 1. For ω* = 0, JtQSS(x, τ) = JQSS(0) = JSS. The shift in time present in Eq. (36) arguably represents some form of memory effect, which we argue makes sense because the nucleation flux describes the rate of emergence of clusters of size x. Obviously, such nuclei do not form instantaneously at this size but grow from a single monomer to a supercritical nucleus under conditions that change continuously due to the ongoing solvent evaporation.

For the remainder of this article, we shall, for simplicity, set J0(x) = 0. This is justified for sufficiently small values of ω* for the regularized gamma function Q entering Eq. (36); in that case, it is very many orders of magnitude smaller than unity and can, for all intents and purposes, be considered to be identical to zero. As we show in Sec. V, for ω* ≳ 5.3, this is no longer true albeit that the exact value depends on the incubation time τi(x) and, therefore, the cluster size x. This we find to be associated with a vanishing and ultimately negative (“virtual”) induction time. As a result, for ω* > 5.3, J0(x) can no longer be justifiably set equal to zero. Remarkably, accounting for a non-zero J0(x) does not remedy this issue because the resulting expression still lacks a lag time.

The origin of this problem lies in an inconsistency in the expression for the stationary nucleation flux that we require in the derivation of Eq. (35). This expression itself, as it turns out, does not accurately account for the (very) early times, as it incorrectly predicts a non-zero nucleation flux at the initial time τ = 0. As far as we are aware, this inconsistency is present in all results on the kinetics of classical nucleation theory obtained using similar singular perturbation techniques under both stationary and non-stationary ambient conditions.26,52,54,57,58 It has been suggested that it is a direct consequence of the use of a singular perturbation approach.52 See supplementary material II for additional details. We postpone an in-depth discussion for what values of the non-stationary index ω* our results remain accurate to Secs. V–VII and first discuss the implications of Eq. (35).

Equation (35) makes explicit both the initial and late time behavior, which we illustrate in Fig. 4. Instead of evaluating the dimensionless surface tension and the initial supersaturation, we somewhat arbitrarily set the values of the critical cluster size at n* = 152 and the nucleation barrier ΔW*(τ = 0) = 70 at time zero. The surface tension corresponding with these values is consistent with a relevant model system for organic electronics and corresponds with that of fullerene dissolved in carbon disulfide for an initial degree of supersaturation of S(0) = 2.7.59 These values are chosen to ensure that the constraints expressed in Eq. (25) hold, i.e., that the critical cluster size and the nucleation barrier drift only very weakly with time. For different initial parameter values, we find qualitatively similar results if the constraints in Eq. (25) hold and if we take into account that the non-stationarity index depends implicitly on the initial conditions as ω* ∼ ΔW*(τ = 0)a, where the value for the exponent depends on the aggregation kinetics, e.g., for reaction-limited aggregation a = 5/2 and a = 3 for diffusion-limited aggregation. Hence, we only show results for a single initial condition.

FIG. 4.

Left: The nucleation rate Eq. (35) evaluated at reduced cluster size x = 1.5 and normalized by quasi-steady-state nucleation flux JQSS(τ = 0) [see Eq. (38)] as a function of the scaled time τ. The critical cluster size is set to n* = 152, and the initial nucleation barrier is set to ΔW*(τ = 0) = 70. We include data for the non-stationarity index values of ω* = 0.005 (red), 0.05 (green), 0.5 (orange), and 1 (blue), discussed in more detail in Sec. VI. Crosses are the numerical solution of the discrete master equations for identical parameter values with a time interval Δτ = 1/4 between the data points, discussed in more detail in Sec. VI. Right: The nucleation rate Eq. (35) is normalized by the renormalized, time-shifted quasi-steady-state nucleation rate Eq. (36), both evaluated at the same reduced cluster size x = 1.5. The color coding is the same as in the left panel.

FIG. 4.

Left: The nucleation rate Eq. (35) evaluated at reduced cluster size x = 1.5 and normalized by quasi-steady-state nucleation flux JQSS(τ = 0) [see Eq. (38)] as a function of the scaled time τ. The critical cluster size is set to n* = 152, and the initial nucleation barrier is set to ΔW*(τ = 0) = 70. We include data for the non-stationarity index values of ω* = 0.005 (red), 0.05 (green), 0.5 (orange), and 1 (blue), discussed in more detail in Sec. VI. Crosses are the numerical solution of the discrete master equations for identical parameter values with a time interval Δτ = 1/4 between the data points, discussed in more detail in Sec. VI. Right: The nucleation rate Eq. (35) is normalized by the renormalized, time-shifted quasi-steady-state nucleation rate Eq. (36), both evaluated at the same reduced cluster size x = 1.5. The color coding is the same as in the left panel.

Close modal

The left and right panels of Fig. 4 show the time dependence of the nucleation flux J(x, τ) for an increasing rate of evaporation, represented by increasing values of the non-stationarity index ω* = 0, 0.005, 0.05, 0.5, and 1.0, and the timescale τs over which ω* drifts with time as τs = , 1.1 × 105, 1.1 × 104, 1.1 × 103, and 5.6 × 101, respectively. We normalize this quantity (left) by the constant quasi steady-state flux JQSS(τ = 0) and (right) the time-dependent renormalized and time-shifted version of it, JtQSS(τ), given in Eq. (36). As both JtQSS(τ) and J(x, τ) are invariants of x, if we express them in terms of a function of the shifted time ττi(x), changing x only affects the time delay, it suffices to show our results for a single, arbitrarily chosen scaled cluster size. We set x = 1.5. The symbols in the left panel are (color-matched) numerical calculations of the (discrete) nucleation master equation and show excellent agreement for the values of ω* shown. We discuss the numerical calculations in more detail in Sec. VI. The agreement deteriorates somewhat for high values of ω*, which of course is not surprising given that Eq. (25) demands that τs1ε2ω*(ΔW*)1ω*1. For the largest value of ω* shown, the time scale τs = 1.1 × 101 is not particularly large, so we should expect that the theoretical expressions are less accurate. Still, we find that the deviations are about 20% or smaller for the tQSS flux for scaled times τ > 12. Hence, our results still provide reasonably accurate predictions, even if the underlying assumptions are no longer necessarily justified.

The left panel of Fig. 4, where the nucleation flux is scaled to the (constant) QSS value at time zero, shows that in the absence of evaporation and ω* = 0 (purple, indistinguishable from the red curve), a constant asymptote is reached after some induction time. The asymptote is given by the steady-state value JQSS(0) = JSS. In the presence of evaporation, for ω* > 0, the flux in this late time regime increases exponentially with time, as we in fact expect from Eq. (36). For very small values of ω*, there is a weak drift in the time window shown in the figure, resulting in what looks like a pseudo plateau. For sufficiently rapid solvent evaporation, particularly if ω* ≫ 1, the nucleation flux does not seem to exhibit this kind of pseudo plateau but increases exponentially immediately after the incubation time. Considering that ω* is a ratio of time scales, this should be expected. For ω* < 1, nucleation is more rapid and dictates the kinetics. For ω* > 1, solvent evaporation is more rapid, so the nucleation kinetics is dictated by solvent evaporation.

The right panel shows that normalized by the time-dependent and rapidly increasing tQSS flux, the nucleation flux reduces to a sigmoidal shape. The distinction between the early and late time regimes, which becomes less evident in the left panel for high ω* due to the apparent absence of the pseudo-plateau, is still evident in the right graph. From this, we conclude that the quasi steady-state nucleation flux remains the central quantity for the late stages of the nucleation process, except that it is renormalized by a factor Γ(ω* + 1) and evaluated at some earlier time [see Eq. (36)]. The magnitude of this renormalization factor grows rapidly with ω* only if ω* > 1, but remains within 12% of unity if the evaporation time scale is much larger than the nucleation time scale (0 < ω* < 1).

Let us now focus specifically on early time behavior. The right panel in Fig. 4 shows that initially the normalized nucleation flux remains vanishingly small for some time, after which it steeply rises to reach the tQSS value. That this must be so is actually evident by considering the limiting behavior of the regularized gamma function in Eq. (35), limz→0Q(n, z) = 1 for τ. This time delay may be considered a lag or induction time that evidently depends on the magnitude of ω*. With an increasing evaporation rate, the curves become steeper and the time lag shortens. Clearly, solvent evaporation influences the early time behavior, but only if it is sufficiently fast, considering the fact that the curves obtained for ω* ≤ 0.05 overlap.

For vanishing solvent evaporation ω* → 0, Eq. (35) reduces to a transient time-shifted quasi-steady-state nucleation flux, realizing that Q(ω*+1,z)exp(z)1+O(ω*),56 yielding
J(x,τ)=JQSS(ττi(x)+τcor)ee(ττi(x)),
(43)
expressed as a product of the late-time quasi-steady-state nucleation flux JQSS(ττi(x) + τcor) evaluated at some shifted time, and a so-called double exponential that describes the initial transient regime, i.e., the sigmoidal curves shown in the right panel of Fig. 4. For stationary nucleation, JQSS(ττi(x) + τcor) = JSS, in which case Eq. (43) correctly reproduces the known expression for the transient stationary flux.26,52,53

The late-time, time-shifted QSS flux is reached after some initial lag period.21,32 Equation (43) suggests that even in the regime where the QSS approximation is justifiable, ω* → 0, the QSS nucleation flux should actually be evaluated at some shifted time ττi(x) + τcor, incorporating the earlier mentioned memory effect. Such a memory effect is not included in the standard QSS theory, even though it is a natural consequence of the existence of a lag time in the kinetics of nucleation; it simply accounts for the fact that the nucleation flux cannot instantaneously respond to a change in the ambient conditions.21 We find that the “standard” QSS flux JQSS(τ) and the time-shifted QSS flux obtained from our theory JQSS(ττi(x) + τcor) are only identical in the stationary limit for ω* = 0. The relative difference between the two fluxes scales exponentially with ω*(τi(x) − τcor). Similar arguments apply if we consider some other source of non-stationarity instead of solvent evaporation.26 In other words, the QSS flux is never correct, except for the trivial case of stationary nucleation.

For sufficiently slow solvent evaporation, however, the difference between the QSS and time-shifted QSS flux does become small. Sufficiently slow here means that the height of the nucleation barrier varies over much longer times than the nucleation time scale. For cluster sizes moderately close to the critical size x = 1 and barrier heights of the order of 10–100 kBT, we conclude from Eq. (39) that the incubation time τi(x) is of the order of the lifetime of a critical cluster, and so is the timescale τcor, as can be deduced from Eq. (41), which is the contribution to the memory effect that has its origin in the evaporation. Hence, in the regime where the QSS flux is expected to be reasonably accurate, if the non-stationarity index ω* is very small, the difference between the quasi-steady nucleation flux JQSS evaluated at τ or at ττi(x) + τcor should be small. Practically, this means that the impact of the memory effect implicit in Eq. (43) should be expected to be minor in the regime where the QSS approximation is thought to be warranted, presumably explaining why the QSS flux without the time-shift is still being used in the field.21,60

In Sec. V, we make the shift in the early time behavior observed in Fig. 4 (right) explicit by introducing the generalized induction time as a measure for the duration of the early time regime and find it to depend in a complex fashion on the value of the non-stationarity coefficient ω*.

This section focuses attention on the effect of solvent evaporation on early time behavior. Here, we calculate a generalized induction time, which we find to depend non-trivially on the solvent evaporation rate via the non-stationarity coefficient, ω. A common approach to quantify the duration of the early time regime under conditions of stationary nucleation is to calculate the lag or induction time that is associated with the initial transient regime before some sort of steady state is reached.21,42,61 For stationary nucleation, the induction time, θ, is defined implicitly by the integral equality,21,42,58
limtθtdtJSS=limt0tdtJ(x,t).
(44)
According to this definition, the induction time marks the point in time after which integration of a steady-state nucleation rate yields the same number of clusters per unit of volume larger than the reduced size x as for the actual nucleation rate, starting from the time t = 0 at which the quench takes place. We opt to use this definition because it is valid for clusters of any reduced cluster size x ≥ 1. From an experimental point of view, it is also more straightforward to determine, as it is not limited by the observation of clusters of the critical size. Such clusters are not necessarily directly observable, e.g., because they tend to be too small or present in a concentration that is too low for direct observation.62 Of course, the induction time defined via Eq. (44) is only valid under stationary conditions but has a straightforward generalization to our non-stationary setting,
limτθτdτJtQSS(x,τ)=limτ0τdτJ(x,τ),
(45)
where on both sides of the equal sign, the limit τ must be taken simultaneously. This expression we construct by demanding that the integrated nucleation flux (given by the right-hand side) equals the integrated late-time time-shifted quasi-steady nucleation flux (the left-hand side) if it were initiated at the induction time θ. Here, we presume that J0(x) = 0 in the nucleation flux Eq. (35). Notice that the integrals on both sides of the equal sign diverge but do so identically by construction because the nucleation flux J(x, τ′) approaches JtQSS(x, τ′) asymptotically for τ; see also the discussion near Eq. (35).
We combine our analytical expressions for the nucleation flux Eq. (35), presuming J0(x) = 0, and the tQSS flux Eq. (36) with Eq. (45), and solve for the induction time θ. To do this, we evaluate the integrals in Eq. (45) separately. We evaluate the integral on the right-hand side of Eq. (45) by first inserting the expression for the nucleation flux Eq. (35). In evaluating the resulting expression, we make use of the work by Shneidman,26 who calculated a similar integral to find the number density of clusters larger than some cluster size x, arriving at
0τdτJ(x,τ)=JQSS(0)eω*τcorE1ω*e(ττi(x)),
(46)
where JQSS(0) = JSS given by Eq. (11) and Ep is the exponential integral of order p.56 To evaluate the integral on the left-hand side of Eq. (45), we make use of the QSS nucleation flux Eq. (11) and the tQSS nucleation flux Eq. (36), which allows us to write the latter as
JtQSS(x,τ)=Γ(ω*+1)JQSS(ττi(x)+τcor)=Γ(ω*+1)JQSS(0)eω*τcoreω*ττi(x).
(47)
Integration with respect to τ′ yields
θτdτJtQSS(x,τ)=Γ(ω*)JQSS(0)×eω*τcoreω*(ττi(x))eω*(θτi(x)),
(48)
where we use the recurrence relation for the Gamma function z−1Γ(z + 1) = Γ(z). Combined with Eqs. (45) and (46), and after a reordering of terms, we find
Γ(ω*)eω*(θτi(x))=limτΓ(ω*)eω*(ττi(x))E1ω*e(ττi(x)),
(49)
with the left hand side independent of the time τ.
Both bracketed terms on the right-hand side of this expression diverge but do so at the same rate and turn out to cancel each other exactly. This follows from the fact that the argument of the exponential integral E1ω* becomes zero in the late-time limit τ. We conveniently introduce the series expansion for the exponential integral as56 
Ep(z)=zp1Γ(1p)k=0(z)kk!(1p+k),
(50)
of which the first term is identical to the first term on the right hand side in Eq. (49). This reduces Eq. (49) to
Γ(ω*)eω*(θτi(x))=limτk=0e(ττi(x))kk!(ω*+k)=1ω*,
(51)
where the last equality follows from the observation that under conditions of solvent evaporation, ω* ≥ 0, only the k = 0 term survives in the limit τ. From Eq. (51), we find an explicit expression for the induction time as a function of the non-stationarity index ω*,
θ(x,ω*)=τi(x)1ω*lnΓ(ω*)1ω*lnω*=τi(x)1ω*lnω*Γ(ω*).
(52)
The contribution to the induction time due to solvent evaporation also emerges in the time shift of the tQSS flux, Eq. (36), if we absorb its prefactor Γ(1 + ω*) into the time shift. Hence, it appears to be a fundamental, non-stationary contribution to the induction or lag time in non-stationary nucleation.
For the sake of clarity, we express Eq. (52) in terms of the induction time under stationary conditions θ0(x) ≡ θ(x, ω* = 0) = τi(x) + γE ≥ 1,58 with γE = 0.577… the Euler–Mascheroni constant,56 as
θ(x,ω*)=θ0(x)γE+1ω*lnω*Γ(ω*).
(53)
This expression indeed reduces to the stationary induction time for the limiting case ω* → 0, which can be verified using the expansion56 
lnΓ(z)=γEzlnz+k=1zkln1+zk.
(54)
Equation (53) is novel and is the second main result of this paper, the first being the nucleation flux Eqs. (35) and (36). The shift in the induction time in Eq. (53) due to solvent evaporation is independent of the reduced cluster size x since ω* is independent of this quantity. Consequently, the relative impact of solvent evaporation on the induction time must decrease with increasing x because Eq. (53) is an increasing function of x. To shed light on how Eq. (53) depends on the evaporation rate that is implicit in the quantity ω*, we focus on the limit ω* ≪ 1. For that purpose, we expand Eq. (53) near the stationary case ω* = 0, yielding the lowest order in powers of ω*,
θ(x,ω*)=θ0(x)π212ω*+Oω*2,
(55)
showing that the function is well-behaved in the limit ω* → 0.

It seems that the induction time θ decreases from its stationary value θ0 to an increasing value of ω*. This agrees with the expectation that solvent evaporation must increase the supersaturation with time, which in turn must induce an increasingly rapid nucleation and hence lead to a reduction of the associated time scale. In Fig. 5, we plot the effect of solvent evaporation on the induction time θθ0 using the non-stationarity index ω* as the control parameter. Here, the blue, solid line is the full result of Eq. (53), and the dashed lines are the linearized result of Eq. (55). Outside of the linearized regime, Eq. (55) overestimates the actual reduction of the induction time.

FIG. 5.

The difference between the induction time at non-zero solvent evaporation θ and at zero solvent evaporation θ0 expressed in reduced time units as a function of the non-stationarity coefficient ω*. The blue solid line corresponds with the full result Eq. (53); the black dashed line is the linearized result Eq. (55).

FIG. 5.

The difference between the induction time at non-zero solvent evaporation θ and at zero solvent evaporation θ0 expressed in reduced time units as a function of the non-stationarity coefficient ω*. The blue solid line corresponds with the full result Eq. (53); the black dashed line is the linearized result Eq. (55).

Close modal

From Eq. (55) and Fig. 5, we conclude that the decrease in induction time is of the order ω* for ω* ≪ 1. Hence, as θ0 ≥ 1, a discernible shift is only observable if ω* is sufficiently large compared to the induction time under stationary conditions, θ0(x). This is only the case if the time scale for solvent evaporation is of similar magnitude as the nucleation time scale, multiplied by the number of molecules in the critical cluster. See also Eq. (22). If that happens to be the case, then for some value of ω*, the induction time θ(x, ω*) may actually vanish and become negative beyond that value. The root cause of this lies in our expression for the nucleation flux Eq. (35), which should be zero for very early times but is not. Apparently, the magnitude of this inconsistency becomes magnified with an increasing value of ω*. We stress that this cannot be amended by invoking a non-zero value for J0(x) in Eq. (35) because this turns out to not affect the value of ω* for which the induction time becomes negative. In our implementation, Eq. (35) incorrectly predicts a non-vanishing and non-negligible nucleation flux for supercritical clusters at τ = 0+. It is not clear to us how to remedy this issue.

Actually, causality dictates that the time lag must always be non-zero and always be positive, for clusters cannot instantaneously grow from monomeric size to the size of a (super)critical cluster. This suggests that the smallest possible lag time should be identical to the time required for the first monomer to grow to critical sizes. As we show in supplementary material I, Eq. (S.4), this growth time is always positive, finite, and independent of the solvent evaporation rate described by the quantity ω*.54 For values of ω* that yield an induction time smaller or approximately equal to the (positively valued) subcritical growth time, Eq. (53) should already be treated with caution. For values of ω* that yield an induction time much larger than the subcritical growth time, Eq. (53) should be accurate.

Taking Eq. (53) at face value, we argue that it provides an upper bound on ω* for the validity of our results for the nucleation flux Eq. (35). Numerical inversion of Eq. (53) shows that the value of ω* for which the induction time becomes negative grows very rapidly with an increasing value of θ0(x). Considering the physically lower bound on the stationary induction time, θ0(x) = 1, we find that beyond ω* ≈ 5.3, the induction time becomes negative. On the other hand, if θ0(x) = 10, we find this to happen for ω* ≈ 6 × 104. The fact that these values are (much) larger than unity means that we expect the theory to break down under conditions where nucleation is predominated by the effect of solvent evaporation, consistent with Eq. (24) and the constraints expressed in Eq. (25).

Having discussed the implications of our result on the induction time, we next focus on a comparison of the theoretical results presented in this section with the numerical results in Fig. 4 in order to verify how accurate Eq. (53) is for high values of ω*. In Fig. 6, we plot the nucleation flux obtained from numerical calculations scaled to the theoretical tQSS flux, Eq. (36), allowing for comparison to the results shown in the right panel of Fig. 4. Note that the curves approach unity for low ω*, so the tQSS flux is accurate, but this is not true for the larger values of ω* shown, where in the late time (after the induction time), the tQSS overestimates the numerical nucleation flux by about 10%–20%. Qualitatively, the agreement is rather good, however, showing that the induction time must indeed decrease with an increasing non-stationarity index ω*. A quantitative comparison requires the explicit integration of the implicit definition of the induction time in Eq. (45) using the numerically obtained nucleation flux. This is only possible if we can extract the full nucleation flux, including the wind-up period, and the late-time nucleation flux, excluding the wind-up period. Unfortunately, based on Fig. 6, we must conclude that the theoretical tQSS flux is not sufficiently accurate for this purpose, and it is also not clear how to extract the late-time nucleation flux from the numerical nucleation flux alone. Therefore, we do not compare the results quantitatively.

FIG. 6.

The nucleation rate Jsim(x, τ) as obtained from numerical calculations at scaled cluster size x = 1.5, normalized by the tQSS nucleation rate [Eq. (36)]. The initial critical cluster size is set to n* = 152, and the initial nucleation barrier is set to ΔW*(τ = 0) = 70. The colors correspond to different values for the non-stationarity index ω*(0) = 0 (purple), ω*(0) = 5 × 10−2 (green), ω*(0) = 5 × 10−1 (orange), and ω*(0) = 1 (blue).

FIG. 6.

The nucleation rate Jsim(x, τ) as obtained from numerical calculations at scaled cluster size x = 1.5, normalized by the tQSS nucleation rate [Eq. (36)]. The initial critical cluster size is set to n* = 152, and the initial nucleation barrier is set to ΔW*(τ = 0) = 70. The colors correspond to different values for the non-stationarity index ω*(0) = 0 (purple), ω*(0) = 5 × 10−2 (green), ω*(0) = 5 × 10−1 (orange), and ω*(0) = 1 (blue).

Close modal

In Sec. VI, we investigate in detail under what conditions our results should be expected to apply. In particular, we discuss how the rate of solvent evaporation influences the accuracy of our prediction, how the initial supersaturation affects our results, and return to the comparison of our results with the predictions of the QSS flux and with numerical calculations.

Before summarizing our findings, we discuss in more detail the impact of our assumptions and approximations on the predictions of the theory. The predictions presented in Secs. IV and V are based on an extension of classical nucleation theory while presuming (i) a nucleation barrier that is considerably larger than the thermal energy, as well as (ii) a solvent evaporation rate that is sufficiently slow to give an approximately constant non-stationarity index ω* in an appropriately chosen co-moving coordinate system. This coordinate system is obtained by scaling the time and the reaction coordinate (the size of a cluster) to the time-dependent lifetime of a critical nucleus tcrit(t) and the critical cluster size n*(t), respectively. These assumptions allow us to treat the reciprocal of the nucleation barrier, ɛ2, and the critical cluster size, n*, as constants in the co-moving coordinate system. However, strictly speaking, solvent evaporation causes these quantities to be time dependent even in the co-moving coordinate system. In addition, and this would be an assumption (iii), in our calculations, nucleation and evaporation start at some initial finite supersaturation, whereas in experiments, they commence under undersaturated conditions.

The question then arises under what practical conditions are these three assumptions justified? To start addressing this question, assumption (i) is implicit in classical nucleation theory and implies that for the time scales considered, the solution should never be close to a spinodal, as the vicinity of a spinodal implies a vanishing interfacial tension and hence a vanishing nucleation barrier.63–65 Given that for a phase transition in which molecules in a dissolved state drop out of solution and form a solid, crystalline state, a spinodal is not thought to exist, we do not expect a low surface tension under any conditions.66,67 We also note that colloidal dispersions do tend to have a very low surface tension if the particles are sufficiently large,68,69 implying that for these, the regime where CNT applies may be limited. For such particles, a metastable gas–liquid spinodal might in addition be hidden under the binodal, as is known for aqueous solutions of certain proteins70 and is suspected for dispersions of fullerenes under appropriate conditions.71 This scenario is explicitly excluded from our analysis.

To investigate under what conditions assumption (ii), expressed by the set of conditions Eq. (25), remains reasonable, we need to analyze the expression for the non-stationarity index ω in Eq. (33) as expressed in terms of the volumetric evaporation rate and the initial nucleation barrier height, which is related to the initial supersaturation. As discussed in Sec. III, the conditions in Eq. (25) can actually be expressed as a single constraint on the time scale over which the non-stationarity index ω drifts, τs ≫ 1, as defined in Eq. (34). Perhaps more insightful is to combine the constraint on ω of Eq. (25) with Eq. (33), resulting in
1ωdωdτ=531τs+τ1,
(56)
indicating that the non-stationarity index decreases with time, since for solvent evaporation, ω > 0. Equation (56) shows that our predictions should indeed be valid if τs ≫ 1. Here, τs = 2/3 × ɛ2ω(0) with ω(0) = αtcrit(0)n*(0), with ɛ the reciprocal square root of the barrier height, α the reciprocal drying time of the thin film, tcrit the lifetime of a critical cluster, and n* the critical cluster size. We note that ω(0), by its very definition, depends on the barrier height. For the reaction-limited model used throughout this work, ω(0)ΔW*5/2, whereas for a diffusion-limited aggregation model, ω(0)ΔW*3.

Based on Eq. (56), we might conclude that our results should also apply for τ ≫ 1, even if τs is not much larger than unity. Because we presume a fixed value for ω*ω*(0) with ω*=(1n*1)ω for our theoretical analysis, this is not the case. Indeed, even if ω* only drifts very weakly with time, the deviation between our prediction for the nucleation flux and the actual nucleation flux, which can, e.g., be obtained by numerically solving Eqs. (2) and (15), must accumulate with time. Equation (56) actually suggests a natural way to correct for this: we could a posteriori make the non-stationarity index ω*, the critical cluster size n*, and ɛ time-dependent, similar in spirit to how the quasi-steady-state nucleation flux is derived from the steady-state nucleation flux. For reasons of self-consistency, these quantities should depend on the shifted time ττi(x), with τi(x) the incubation time, as was in fact suggested by Shneidman.26 Upon introducing these two ad hoc corrections into our work, we deduce from Eq. (56) that our results are expected to become accurate if τ ≫ 1 also.

The question now arises how strongly does solvent evaporation actually affects nucleation under experimentally relevant conditions? Solvent evaporation need not be fast in an absolute sense for it to affect nucleation, only fast relative to nucleation. This is certainly the case very close to the binodal, where the non-stationarity index ω* can become very large as both the number of molecules in a critical cluster and the lifetime of a critical cluster diverge, causing even slow solvent evaporation to be rapid relative to nucleation. This situation presents itself during a gradual quench from the stable to the metastable solution, relevant to experiments for the casting of thin film organic electronics. Very close to the binodal; however, nucleation is so slow that it cannot be observed on experimental time scales. Hence, as also discussed in Sec. II, the magnitude of the effect of solvent evaporation on nucleation depends on the conditions present at the moment that the first (supercritical) nuclei form and is, therefore, related to our final assumption (iii): how the choice of the initial supersaturation affects our results.

Studying how the choice of initial conditions affects our results requires a reasonable estimate of the initial nucleation barrier height and a reasonably good understanding of the nucleation kinetics so that we can estimate the attachment frequencies. These attachment frequencies can be determined by combining Eq. (29) with parameter estimates from either experimental or atomistic numerical studies.72 As far as we are aware, these are not known for the organic compounds of interest to us. Hence, at present, we cannot accurately estimate the instant in time we can set to t = 0, nor, consequently, the magnitude of the evaporation-related effects on nucleation under experimentally relevant conditions.

In order to make headway in spite of this, we discuss how changing the initial supersaturation S(0) changes our results qualitatively. Within our theory, this is equivalent to changing the initial barrier height ΔW*(0) ∼ ln−2S(0). Nucleation barriers for homogeneous nucleation observable on experimental time scales are generally estimated between 10 and 100 kBT.21 For higher nucleation barriers, nucleation almost exclusively commences via heterogeneous nucleation.21 Hence, let us take 100 kBT as a reasonable estimate for the maximum barrier that results in nucleation on experimental timescales.

As long as the conditions corresponding to the initial supersaturation or barrier height result in the constraint τs ≫ 1 presented in Eq. (34) to be valid and the induction time Eq. (53) to be positively valued, the expressions for the nucleation rate Eqs. (35) and (36) and the result for the induction time in Eq. (53) are expected to be accurate. Upon changing the initial supersaturation for fixed values of the other kinetic and thermodynamic parameters, the results shown in Fig. 4 for the nucleation flux and in Fig. 5 for the induction time remain correct if we account for implicit changes in two model parameters. First, we note that the lifetime of a critical cluster tcrit scales linearly with the barrier height, tcrit ∼ ΔW*. Since all results are presented in terms of the scaled time τ instead of the actual time t, this change is actually implicit in the results presented. Second, the non-stationarity index ω* depends on the initial nucleation barrier and evaporation rate α, as ω*αΔW*5/2. Clearly, the impact of solvent evaporation on nucleation decreases rapidly (super linearly) with an increasing initial degree of supersaturation. For example, if the initial barrier height changes from 100 to 10 kBT at a fixed evaporation rate α, ω* decreases by a factor of about 3 × 102.

This evidently presumes that homogeneous nucleation is the dominant nucleation mechanism for the experimental systems relevant to us. In a solution confined to a thin film geometry, this might not always be the case, and homogeneous nucleation could be preempted by heterogeneous nucleation. This might occur if solvent evaporation is very rapid, resulting in the accumulation of solute near the solution–air interface.28 Alternatively, homogeneous nucleation becomes much less probable simply because the available space (volume) for homogeneous nucleation relative to the available space for heterogeneous nucleation (the substrate and solution–air interfaces) decreases with decreasing thickness of the wet film. As we discuss in the  Appendix, the extension of our work to heterogeneous nucleation is straightforward to carry out, and all results in this paper apply with only minor modifications.

To further investigate the limitations of our theory, we return to numerics, reminding the reader that we compared our theory with numerical calculations of Eq. (57) in Fig. 4 without discussing them in detail. We found in Sec. V that the induction time vanishes at very small supersaturation, so we expect that we only need to consider the tQSS nucleation flux of Eq. (36) for experimental situations. Hence, we now focus attention on the nucleation flux at times much later than the induction time and discuss the differences between the tQSS and the QSS fluxes.

For the numerical calculations, we use discrete kinetic equations, which are numerically simpler to solve. The evolution equations for the cluster densities then read21,73
ρn(t)t=Jn1,nJn,n+1+Kn(t)
(57)
where the subscripts indicate the (discrete) cluster size n, Jn−1,n the cluster flux density in n-space, and Kn(t) again the evaporation source term defined by Eq. (9). The cluster flux density is given by21 
Jn1,n=kn1,n+(t)ρn1(t)Feq(n1,t)Feq(n,t)ρn(t).
(58)
Here, kn,n+1+(t) is the attachment frequency for a monomer to an n-mer, and Feq is again the “equilibrium” cluster densities defined in (3). We solve Eq. (57) using a standard backward Euler method with time step Δτ = 0.05, assuming nmax = 10 000 ≫ n* as the maximum cluster size. This value for nmax is sufficiently large not to affect our results for the nucleation flux.

The numerical results presented here use the same arbitrarily chosen nucleation barrier of ΔW*(0) = 70 and critical cluster size of n*(0) = 152 we use in Sec. IV and compare different values for the non-stationarity index ω*. These values for ΔW*(0) and n*(0) ensure that the nucleation barrier is sufficiently high such that the high-barrier assumption holds and that our analytical results should principally be accurate. As long as the constraint presented in Eq. (56) is satisfied, other values for the initial conditions yield results that are qualitatively the same if we account for the change in the non-stationarity index ω*ΔW*5/2 at a fixed evaporation rate, as discussed earlier.

Before we can discuss the comparison of our theoretical and numerical results, we need to point out possible differences between the solutions to the continuous and discrete nucleation equations.58,74 Specifically, they may result in slightly different induction times. The magnitude of these differences turns out to depend strongly on the aggregation model used. For the reaction-limited growth model we use, these differences are known to be very minor or even negligible, justifying a comparison between results from the discrete and continuum equations.26,58

We expect that only the tQSS flux and not the initial wind-up period can be observed under experimental conditions, where at time zero, we cross the binodal during evaporation. The reason is that the induction time vanishes very close to the binodal, referring to our discussions earlier in this section and Sec. V. The question arises how much more accurate is the tQSS flux compared to the QSS flux? This is what we focus on in Fig. 7, where we plot the numerical results from the same calculations as those shown in the left panel of Fig. 4 (crosses), so n* = 152 and ΔW(0) = 70, and compare them with the QSS flux Eq. (11) (dotted, color-matched), as well as the tQSS flux of Eq. (36) (dashed, color matched). As discussed earlier, we let ω*, n*, and ɛ depend on the shifted time ττi(x) to account for the weak drift in these quantities. In Fig. 7, we do not include the initial time for τ < 8 but focus on the theoretical expressions that are accurate at later times. For these later times, we find that for low values of ω* (red and green curves), the QSS and tQSS fluxes are, for all intents and purposes, identical. This we also expect based on the small-ω* limit of Eq. (36). See also Eq. (43). At τ = 12, we are (far) beyond the induction time and find that at this instant in time, the QSS flux overpredicts the nucleation flux by about 70% for ω* = 0.5 (orange) and by about 90% for ω* = 1.0 (blue), at least for larger values of ω*. The tQSS flux also overpredicts the nucleation flux but does so to a significantly smaller extent, about 10% for ω* = 0.5 and about 20% for ω* = 1.0.

FIG. 7.

The nucleation flux J(x, τ) evaluated at the scaled cluster size x = 1.5, normalized by the JQSS(τ = 0) [see Eq. (38)] as a function of the scaled time τ. Crosses are the numerical solutions with a time interval between the data points of Δτ = 1/4; dotted lines are the QSS flux Eq. (38); dashed lines are the tQSS flux given by Eq. (36). We set the critical cluster size n* = 152 and the initial nucleation barrier ΔW*(τ = 0) = 70. We show results for ω*(0) = 5 × 10−3, τs = 1.1 × 105 (red), ω*(0) = 5 × 10−2, τs = 1.1 × 104 (green), ω*(0) = 5 × 10−1, τs = 1.1 × 103 (orange), and ω*(0) = 1.12, τs = 5.6 × 101 (blue).

FIG. 7.

The nucleation flux J(x, τ) evaluated at the scaled cluster size x = 1.5, normalized by the JQSS(τ = 0) [see Eq. (38)] as a function of the scaled time τ. Crosses are the numerical solutions with a time interval between the data points of Δτ = 1/4; dotted lines are the QSS flux Eq. (38); dashed lines are the tQSS flux given by Eq. (36). We set the critical cluster size n* = 152 and the initial nucleation barrier ΔW*(τ = 0) = 70. We show results for ω*(0) = 5 × 10−3, τs = 1.1 × 105 (red), ω*(0) = 5 × 10−2, τs = 1.1 × 104 (green), ω*(0) = 5 × 10−1, τs = 1.1 × 103 (orange), and ω*(0) = 1.12, τs = 5.6 × 101 (blue).

Close modal

Clearly, the tQSS nucleation flux (taking into account the time-drifting ω*, n*, and ɛ) provides a significant improvement over the QSS approximation, even at relatively fast solvent evaporation where the assumptions required for its derivation are not necessarily justified. We attribute the eventual decrease in accuracy of the tQSS flux to the (relatively) small value of τs [see Eq. (56)]. The temporal drift in the non-stationarity index ω*, the critical cluster size n*, and in ɛ can in this case not be accounted for a posteriori. For τs small, the temporal drift should actually be taken into account in the calculations presented in Sec. III and in supplementary material I and II. This significantly complicates finding analytical solutions to the nucleation master equation, and we are at present not aware of any method available that allows one to achieve this.

In this work, we study the effect of solvent evaporation on the rate of homogeneous nucleation of a solute in a two-component solution to form a solid precipitate. Solvent evaporation causes the nucleation conditions to change continuously, implying that a description based on stationary nucleation theory cannot apply. In order to explicitly account for the effect of solvent evaporation on the rate at which supercritical clusters emerge, we extend classical nucleation theory to an instantaneous quench in the metastable region upon which nucleation commences under the conditions of ongoing evaporation. Our findings demonstrate that additional timescales emerge that are not present in either stationary or quasi-steady-state approaches that have been put forward in the literature.21 We find that the ratio of two of those time scales dictates how strongly evaporation affects the nucleation process.

The three most important time scales we extract from our theory are (1) the timescale associated with solvent evaporation, (2) the time scale related to the lifetime of a critical cluster, and (3) a generalized induction time, which is a measure for how quickly nucleation responds to an initial quench. The latter two we find to also be influenced by the rate of solvent evaporation. It turns out that the magnitude of the effect of solvent evaporation can be measured by a single quantity that we dubbed the non-stationarity index, which is proportional to the ratio of the lifetime of a critical cluster to the evaporation timescale and proportional to the number of molecules in a critical cluster. If the value of this quantity is very much smaller than unity, evaporation is slow, and the relevant evaporation timescale is large. Hence, solvent evaporation only weakly affects nucleation, and we retrieve a time-shifted, quasi-steady nucleation.21 Surprisingly, this time-shifted quasi-steady-state is never equivalent to the “standard” (not time-shifted) quasi-steady-state limit of nucleation, except in the trivial case of stationary nucleation.21 In the actual quasi-steady limit, we must always account for a delayed response to ambient conditions, which, of course, is a direct consequence of the time lag of nucleation. The customary quasi-steady-state nucleation flux and time-shifted quasi-static flux we obtain are related via a time-shift and a factor of proportionality that scales exponentially with the non-stationarity index. This is not only true for solvent evaporation but can actually be generalized to arbitrary sources of non-stationarity26 and should, therefore, have consequences in other fields where nucleation is non-stationary.

If the non-stationarity index is sufficiently large, then solvent evaporation strongly affects nucleation. The nucleation cascade can, in that case, not sufficiently rapidly accommodate changes in the ambient conditions, causing its rate to increase less rapidly than is to be expected, presuming an instantaneous response. Faster evaporation causes this delay to increasingly strongly affect the nucleation rate. The delay in response to the ambient conditions we find to be related to a generalized induction time for which we have derived an analytical expression. Our prediction for the generalized induction time shows that it must decrease with an increasing evaporation rate. This agrees with our expectation that an increasing supersaturation due to solvent evaporation must result in increasingly faster nucleation. The product of the evaporation rate and the delay time sets the magnitude of the effect of solvent evaporation on nucleation. Hence, the rate of nucleation is governed by two counteracting effects originating from solvent evaporation: while a faster evaporation rate results in an increasingly larger influence of the delay time on the nucleation flux, this delay time itself effectively decreases with increasing evaporation rate. For sufficiently fast evaporation, which we argue to always occur sufficiently close to the binodal of the component dropping out of solution, our predictions indicate that the generalized induction time becomes negative and, in a sense, becomes virtual. This we do not expect to actually happen because our theory should break down for a sufficiently large non-stationarity index. Under those circumstances, the key assumption underlying our theory—that the timescale over which the nucleation barrier drifts with time is much larger than the nucleation time scale—no longer applies.

Experimentally, producing an instantaneous, non-zero supersaturation upon which the solvent evaporation commences seems difficult to achieve. Indeed, the inspiration for this work is the wet deposition of a spin-coated or slot-die cast thin film of a crystalline organic semiconductor used in, e.g., the production of organic transistors,2 where the solution quickly but gradually changes from undersaturated to supersaturated.2 In early times, just after crossing the binodal, nuclei arguably do not emerge because nucleation is exceedingly slow and very close to the binodal. Hence, the distinction between an initially undersaturated solution and an initially supersaturated solution becomes, in a sense, somewhat arbitrary. To connect with spin-coating or slot-dying experiments, we argue that in the theory, the initial supersaturation should be set to the conditions corresponding to the time before the instant at which the first nuclei emerge. The associated effective supersaturation at this kind of “time zero” can even be calculated if a microscopic understanding of the underlying nucleation kinetics is available, informed, for instance, by experiments. Here, we steer away from attempting to precisely pinpoint this quantity and suffice by arguing that any resulting large value of the nucleation barrier will do. Our results depend only weakly on the nucleation barrier at zero time, provided it is sufficiently large compared to the thermal energy.

Theoretically, the analysis of nucleation under time-dependent supersaturation has always been hindered by the presence of multiple time and reaction coordinate (the cluster size) scales that generally both also depend on time and do so differently.24 This causes the kinetic nucleation equations to be exceedingly hard to solve analytically.21,24,52,53 This is also true for the relevant time and reaction coordinate scales within our model. We solved this problem by making use of a coordinate system that is co-moving with respect to both the reaction coordinate and time, following the earlier work of Shneidman.26 To obtain analytical predictions within the prescription of a co-moving coordinate frame, we presume that the rate of change of the nucleation barrier is constant in the co-moving time coordinate. We find that this is reasonable as long as the timescale over which the nucleation barrier drifts in time is much longer than the nucleation time scale.

While this is strictly not true in reality, we can post hoc account for such a drift in time as long as this drift in time is sufficiently slow. This procedure is, of course, reminiscent of what is performed to obtain the quasi-steady-state nucleation flux.21,26 The rate of change of the rate of change of the nucleation barrier that emerges sets a timescale that depends on the ratio of the nucleation barrier height to the non-stationarity index. If this timescale is sufficiently large (on the scale of the nucleation time), then the drift of the nucleation barrier in terms of the co-moving time frame can be ignored. Our numerical calculations show that even if this timescale is not large, our theory turns out to be much more accurate than the usual quasi-steady-state theory.

To conclude, we believe that our results increase the current understanding of the crystallization or solidification that occurs during the casting of solutions for thin film materials. The rate of solvent evaporation has been used experimentally to exert control over crystallization processes during the solution casting of functional thin films.75,76 How solvent evaporation affects the final crystal structure is not only determined by nucleation itself but also by the subsequent growth processes, which are material specific.77 For example, for perovskite thin film solar cells, a higher evaporation rate results in more and smaller crystalline domains,75,76 whereas for experimental work on a fullerene-based solar cell, the rate of evaporation appears to have little-to-no effect on the mean size of the fiber-like crystal structure.16 To be able to make tangible predictions on crystallized morphology based on the present contribution, we need to account for the relevant growth process as well. This can be achieved by, for example, combining our results with phase-field simulations8 or incorporating them into a nucleation-and-growth model such as the Johnson–Mehl–Avrami–Kolmogorov theory,47 which is often applied in combination with in situ diagnostics of the formation of functional polycrystalline thin films.78 This would allow us to study not only the initial nucleation but also the subsequent crystallite growth as well as the grain statistics of the final polycrystalline structure.

A possible extension of our work would be to account for the additional solutes and solvents that are commonly used to control the morphology and the performance of the devices.11,16 Such additional components are known to affect the nucleation-and-growth process,16 which is believed to originate from composition-dependent evaporation rates and surface tensions. These might also induce novel sources of non-stationarity and provide a possible route to gain better control over the nucleation process and the subsequent growth of the crystallites.

The supplementary material contains the derivation of the scaled cluster densities ν(x, τ) and the nucleation flux J(x, τ) and a discussion on the properties of both quantities.

R.d.B., J.J.M., and P.v.d.S. acknowledge funding by the Institute for Complex Molecular Systems at Eindhoven University of Technology. We are grateful to D. Reguera (University of Barcelona, Spain) for helpful discussions.

The authors have no conflicts to disclose.

René de Bruijn: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Jasper J. Michels: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). Paul van der Schoot: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

A

area (m2)

Feq

dimensionless “equilibrium” cluster distribution

g

(deterministic) cluster growth rate

J

nucleation flux (s−1)

JQSS

quasi-steady-state nucleation flux

JSS

steady-state nucleation flux

JtQSS

renormalized, time-shifted quasi-steady-state nucleation flux

K

evaporation source nucleation master equation (s−1)

k+

attachment frequencies (s−1)

k̄+

dimensionless attachment frequencies

n

cluster size reaction coordinate (number of molecules in cluster)

N

number of solute molecules in the solution

n*

critical cluster size

β = 1/(KBT)

reciprocal thermal energy (J−1)

S

supersaturation

tcrit

lifetime of a cluster (s)

V

volume in which the solution is contained (m3)

v0

molecular volume of solute molecules (m3)

x

reduced cluster size reaction coordinate

X

stretched reduced cluster size

Z

Zeldovich factor

γ

surface tension (J)

ΔW

nucleation free energy (J)

Δμ

chemical potential difference between phases (J)

ζ

non-stationarity function

ν = ρ/Feq

reduced cluster densities

ρ

dimensionless cluster densities

ρ

solubility limit

σ

dimensionless surface tension

τ

reduced time coordinate

τcor

time shift due to no-stationary effects

τi

incubation time

τnon-stat

non-stationary contribution to the lag time

τs

time scale over which ω drifts with time

τstat

part of the incubation time

χ

dimensionless reciprocal volume

ψ

cluster size dependence of ζ

ω

unmodified non-stationarity index

ω*

non-stationarity index, defined as ω*=ω(1n*1)

A

area (m2)

Feq

dimensionless “equilibrium” cluster distribution

g

(deterministic) cluster growth rate

J

nucleation flux (s−1)

JQSS

quasi-steady-state nucleation flux

JSS

steady-state nucleation flux

JtQSS

renormalized, time-shifted quasi-steady-state nucleation flux

K

evaporation source nucleation master equation (s−1)

k+

attachment frequencies (s−1)

k̄+

dimensionless attachment frequencies

n

cluster size reaction coordinate (number of molecules in cluster)

N

number of solute molecules in the solution

n*

critical cluster size

β = 1/(KBT)

reciprocal thermal energy (J−1)

S

supersaturation

tcrit

lifetime of a cluster (s)

V

volume in which the solution is contained (m3)

v0

molecular volume of solute molecules (m3)

x

reduced cluster size reaction coordinate

X

stretched reduced cluster size

Z

Zeldovich factor

γ

surface tension (J)

ΔW

nucleation free energy (J)

Δμ

chemical potential difference between phases (J)

ζ

non-stationarity function

ν = ρ/Feq

reduced cluster densities

ρ

dimensionless cluster densities

ρ

solubility limit

σ

dimensionless surface tension

τ

reduced time coordinate

τcor

time shift due to no-stationary effects

τi

incubation time

τnon-stat

non-stationary contribution to the lag time

τs

time scale over which ω drifts with time

τstat

part of the incubation time

χ

dimensionless reciprocal volume

ψ

cluster size dependence of ζ

ω

unmodified non-stationarity index

ω*

non-stationarity index, defined as ω*=ω(1n*1)

The extension of our work to heterogeneous nucleation is straightforward, as we show next. Let us consider the generic case where the nuclei that form partially wet the solution–air interface or the substrate–solution interface, which we both presume to be flat. Here, for the sake of simplicity, we presume that the nuclei are cap-shaped (see Fig. 8) for a schematic representation. While such cluster shapes seem unphysical for heterogeneous crystal nucleation due to the underlying crystal symmetry and associated anisotropic interfacial tension, it turns out that the phenomenology of this model is still applicable to crystal nucleation.79–81 Hence, we use it for illustrative purposes.

FIG. 8.

Schematic representation of cluster geometries during heterogeneous nucleation. The substrate is located at the bottom and is dashed. The solution–air interface is located at the top. Nuclei attach to either the solution–air interface (top-right) or the solution substrate interface (bottom-left). The clusters are cap-shaped with contact angles ϕsub and ϕair for the solution–substrate and solution–air interfaces, respectively.

FIG. 8.

Schematic representation of cluster geometries during heterogeneous nucleation. The substrate is located at the bottom and is dashed. The solution–air interface is located at the top. Nuclei attach to either the solution–air interface (top-right) or the solution substrate interface (bottom-left). The clusters are cap-shaped with contact angles ϕsub and ϕair for the solution–substrate and solution–air interfaces, respectively.

Close modal
Focusing first on the thermodynamics of heterogeneous nucleation, where we neglect any contribution from a line tension, the nucleation free energy is given by21,
ΔWHet=nΔμ+f(ϕ)σn2/3,
(A1)
where σ is the dimensionless surface tension between the solution and the nucleus [see also Eq. (4) in the main text], and the scaling function f(ϕ) is a function that describes (i) the change in volume-to-surface ratio present in heterogeneous nucleation and (ii) the change in surface energetics as the different types of interfaces have different surface tensions.21,80 In general, additional physics, such as that arising from the presence of a line tension, can be incorporated into the scaling function f(ϕ), in which case the scaling function might also depend on the cluster size n. Equation (A1) shows that the difference between homogeneous and heterogeneous nucleation can generally be accounted for by a renormalized surface tension. This can lower79,80 or increase81 the nucleation barrier with respect to homogeneous nucleation. We refer to Refs. 21, 79, and 80 for details on the scaling function f(ϕ).
Note that we do not need to subtract the self-consistent correction from the free energy to set the free energy to form a monomer to zero, as we did in the case of homogeneous nucleation (see Sec. II in the main text). For heterogeneous nucleation, we must distinguish between a monomer in free solution and a monomeric cluster, that is, a monomer that is adsorbed onto a surface, adding adsorption-free energy to the free energy barrier. The change in chemical potential Δμ we presume is similar to that in homogeneous nucleation; that is, we presume the chemical potential to be given by Eq. (5), where the density ρ(t) acts as the local solute density at the substrate or solution–gas interface and ρ is the solubility limit. All thermodynamic quantities such as the critical nucleation barrier, the critical cluster size, and the Zeldovich factor presented in the main text for homogeneous nucleation apply for heterogeneous nucleation if we account for the renormalized (effective) surface tension by the factor f(ϕ) and if we exclude the self-consistent monomer correction. These quantities are given by
ΔWHet*=13n*2/3f(ϕ)σ=3ε2,
(A2)
nHet*=2f(ϕ)σ3lnS(t)3,
(A3)
and
ZHet=12π2ΔWn2n=n*=f(ϕ)σπ13n*2/3.
(A4)
The expression for the equilibrium distribution does change somewhat and is now given by
Feq(n,t)=ρ0eΔWHet,
(A5)
where ρ0 is the surface density of nucleation sites, which for a flat interface is constant in time as evaporation does not affect this quantity.
The kinetics of heterogeneous nucleation is slightly different since the continuity equation does not need to be supplemented by a source term. The reason is, of course, that the number of clusters per surface area does not change over time due to solvent evaporation. Hence, we find
ρ(n,t)t=J(n,t)n,
(A6)
where J(n, t) is still given by the Zeldovich form Eq. (2). Although solvent evaporation is not explicitly present in the continuity equation, it remains present as a supersaturation in the nucleation barrier and, hence, a model for solvent evaporation is still required.
Combining the thermodynamic and kinetic descriptions of heterogeneous nucleation in the same way as we did in the main body of the article, we find
ν(x,τ)τ+ν(x,τ)ζ(x,τ)=ε2k̄+(x)ΔWHet(x,τ)x+xlnn*τν(x,τ)x+ε22k̄+(x)xν(x,τ)x+ε22k̄+(x)2ν(x,τ)x2,
(A7)
expressed in the co-moving {x, τ} coordinate system, see Sec. III, and must be supplemented by the same boundary conditions as in the homogeneous nucleation case, Eq. (28). Here, we define the function ζ(x, τ) that describes the non-stationarity
ζ(x,τ)=τΔWHet(x,τ)=ωψ(x),
(A8)
see Eq. (22) in the main text. For the non-stationarity coefficient ω, we must account not only for the changing supersaturation but also consider that the function f(ϕ) might, in principle, depend on time. For instance, for crystal nucleation, the surface tensions of the different facets might depend on the (time-dependent) solute concentration. This causes the effective surface tension in Eq. (A1) to become time-dependent. If we, for simplicity, use the illustrative example of hemispherical nuclei shown in Fig. 8, a time-dependent surface tension could be related to a time-dependent surface-nuclei contact angle. This results in a non-stationarity coefficient that can be shown to be given by
ω=n*lnvχ(τ)τσn*2/3fϕϕτ.
(A9)
Here, the first term is due to the non-stationarity originating from the time-dependent supersaturation and is also present in homogeneous nucleation [see (22)]. The second term accounts for the change in the contact angles. Similar contributions can be included if the scaling function f(ϕ) turns out to depend on time.

Presuming that ω remains approximately constant within the co-moving coordinate system described in the main text, the solution to Eq. (A7) for the nucleation flux is equivalent to the expression found for homogeneous nucleation in Eqs. (35) and (36). We must, however, replace the critical nucleus size, the critical nucleation free energy, the Zeldovich factor, and the non-stationarity factor ω with their respective expressions presented in this  Appendix. Moreover, the nucleation attachment mode might be different, so another description of the attachment frequencies might be required. The same holds for the induction time, Eq. (53).

Hence, we conclude that the general theory presented in the main text of this work applies also to heterogeneous nucleation, mutatis mutandis.

1.
L. J.
Richter
,
D. M.
DeLongchamp
, and
A.
Amassian
,
Chem. Rev.
117
,
6332
(
2017
).
2.
J.
Mei
,
Y.
Diao
,
A. L.
Appleton
,
L.
Fang
, and
Z.
Bao
,
J. Am. Chem. Soc.
135
,
6724
(
2013
).
3.
H.
Ling
,
S.
Liu
,
Z.
Zheng
, and
F.
Yan
,
Small Methods
2
,
1800070
(
2018
).
4.
D.
Di Carlo Rasi
and
R. A. J.
Janssen
,
Adv. Mater.
31
,
1806499
(
2019
).
5.
H.
Gaspar
,
F.
Figueira
,
L.
Pereira
,
A.
Mendes
,
J.
Viana
, and
G.
Bernardo
,
Materials
11
,
2560
(
2018
).
6.
C.
McDowell
,
M.
Abdelsamie
,
M. F.
Toney
, and
G. C.
Bazan
,
Adv. Mater.
30
,
1707114
(
2018
).
7.
K.
Zhang
,
Z.
Wang
,
T.
Marszalek
,
M.
Borkowski
,
G.
Fytas
,
P. W. M.
Blom
, and
W.
Pisula
,
Mater. Horiz.
7
,
1631
(
2020
).
8.
J. J.
Michels
,
K.
Zhang
,
P.
Wucher
,
P. M.
Beaujuge
,
W.
Pisula
, and
T.
Marszalek
,
Nat. Mater.
20
,
68
(
2021
).
9.
V.
Negi
,
A.
Lyulin
, and
P.
Bobbert
,
Macromol. Theory Simul.
25
,
550
(
2016
).
10.
J. J.
van Franeker
,
D.
Westhoff
,
M.
Turbiez
,
M. M.
Wienk
,
V.
Schmidt
, and
R. A. J.
Janssen
,
Adv. Funct. Mater.
25
,
855
(
2015
).
11.
J. J.
van Franeker
,
M.
Turbiez
,
W.
Li
,
M. M.
Wienk
, and
R. A.
Janssen
,
Nat. Commun.
6
,
6229
(
2015
).
12.
T.
Lee
,
A. V.
Sanzogni
,
P. L.
Burn
, and
A. E.
Mark
,
ACS Appl. Mater. Interfaces
12
,
40548
(
2020
).
13.
C.
Schaefer
,
J. J.
Michels
, and
P.
van der Schoot
,
Macromolecules
49
,
6858
(
2016
).
14.
C.
Schaefer
,
S.
Paquay
, and
T. C. B.
McLeish
,
Soft Matter
15
,
8450
(
2019
).
15.
V.
Negi
,
O.
Wodo
,
J. J.
van Franeker
,
R. A. J.
Janssen
, and
P. A.
Bobbert
,
ACS Appl. Energy Mater.
1
,
725
(
2018
).
16.
J. J.
van Franeker
,
G. H. L.
Heintges
,
C.
Schaefer
,
G.
Portale
,
W.
Li
,
M. M.
Wienk
,
P.
van der Schoot
, and
R. A. J.
Janssen
,
J. Am. Chem. Soc.
137
,
11783
(
2015
).
17.
L.
Zhu
,
M.
Zhang
,
G.
Zhou
,
T.
Hao
,
J.
Xu
,
J.
Wang
,
C.
Qiu
,
N.
Prine
,
J.
Ali
,
W.
Feng
,
X.
Gu
,
Z.
Ma
,
Z.
Tang
,
H.
Zhu
,
L.
Ying
,
Y.
Zhang
, and
F.
Liu
,
Adv. Energy Mater.
10
,
1904234
(
2020
).
18.
A. A.
Virkar
,
S.
Mannsfeld
,
Z.
Bao
, and
N.
Stingelin
,
Adv. Mater.
22
,
3857
(
2010
).
19.
Y.
Wu
,
S.
Schneider
,
C.
Walter
,
A. H.
Chowdhury
,
B.
Bahrami
,
H.-C.
Wu
,
Q.
Qiao
,
M. F.
Toney
, and
Z.
Bao
,
J. Am. Chem. Soc.
142
,
392
(
2020
).
20.
E. M.
Tennyson
,
T. A. S.
Doherty
, and
S. D.
Stranks
,
Nat. Rev. Mater.
4
,
573
(
2019
).
21.
D.
Kashchiev
,
Nucleation
, 1st ed. (
Elsevier
,
Oxford
,
2000
).
22.
V. A.
Shneidman
and
P.
Hänggi
,
Phys. Rev. E
49
,
641
(
1994
).
23.
V. A.
Shneidman
,
J. Chem. Phys.
127
,
041102
(
2007
).
24.
H.
Trinkaus
and
M. H.
Yoo
,
Philos. Mag. A
55
,
269
(
1987
).
25.
V. A.
Shneidman
,
Phys. Rev. E
84
,
031602
(
2011
).
26.
V. A.
Shneidman
,
Phys. Rev. E
82
,
031603
(
2010
).
27.
D. E.
Bornside
,
C. W.
Macosko
, and
L. E.
Scriven
,
J. Appl. Phys.
66
,
5185
(
1989
).
28.
C.
Schaefer
,
J. J.
Michels
, and
P.
van der Schoot
,
Macromolecules
50
,
5914
(
2017
).
29.
D.
Reguera
and
J. M.
Rubí
,
J. Chem. Phys.
119
,
9877
(
2003
).
30.
D.
Reguera
and
J. M.
Rubí
,
J. Chem. Phys.
119
,
9888
(
2003
).
31.
A.
Kuhnhold
,
H.
Meyer
,
G.
Amati
,
P.
Pelagejcev
, and
T.
Schilling
,
Phys. Rev. E
100
,
052140
(
2019
).
32.
V.
Kalikmanov
,
Nucleation Theory, Lecture Notes in Physics
(
Springer
,
Netherlands, Dordrecht
,
2013
), Vol.
860
.
33.

The validity and accuracy of this equation, specifically related to treating n as a continuous variable, are discussed in Ref. 74.

34.
Y. B.
Zeldovich
,
Acta Physicochim. URSS
18
,
1
(
1943
).
35.
D. T.
Wu
, in
Solid State Physics
, edited by
H.
Ehrenreich
and
F.
Spaepen
(
Academic Press
,
San Diego
,
1996
), Vol.
50
, pp.
37
197
.
36.
J.
Danglad-Flores
,
S.
Eickelmann
, and
H.
Riegler
,
Chem. Eng. Sci.
179
,
257
(
2018
).
37.
F.
Gumpert
,
A.
Janßen
,
C. J.
Brabec
,
H.-J.
Egelhaaf
,
J.
Lohbreier
, and
A.
Distler
,
Eng. Appl. Comput. Fluid Mech.
17
,
2242455
(
2023
).
38.
Z.
Kožíšek
,
K.
Sato
,
S.
Ueno
, and
P.
Demo
,
J. Chem. Phys.
134
,
094508
(
2011
).
39.
N. S.
Tiwari
and
P.
van der Schoot
,
J. Chem. Phys.
144
,
235101
(
2016
).
40.
J. K.
Dhont
,
An Introduction to Dynamics of Colloids
(
Elsevier
,
Amsterdam
,
2003
).
41.
K. F.
Kelton
,
A. L.
Greer
, and
C. V.
Thompson
,
J. Chem. Phys.
79
,
6261
(
1983
).
42.
D. T.
Wu
,
J. Chem. Phys.
97
,
2644
(
1992
).
43.

The factor 2/π is non-standard as the lifetime of a critical cluster (20) is usually found after approximately solving Eq. (15) for all n. The precise form of the factor then depends on what approximations where introduced to solve this equation.

44.
J.
Crank
,
The Mathematics of Diffusion
(
Clarendon Press
,
Oxford
,
1975
).
45.
H.
Risken
,
Fokker–Planck Equation
(
Springer
,
Berlin
,
1996
), pp.
63
95
.
46.

Since the attachment of a monomer to a cluster is a second-order reaction, we expect the attachment frequencies to be given by k+(x, τ) = ρ(1, τ)κ(x, τ), with κ(x, τ) a reaction coefficient associated with the attachment of a monomer to a cluster of size x. As a result, the reduced attachment frequencies are ¯kk+(x, τ) = κ(x, τ)/κ(1, τ). If n* and ɛ can be treated as essentially independent of time, then, for all relevant models for the aggregation kinetics known to us, this holds for κ(x, τ)/κ(1, τ) as well.21,32,51

48.
L.
Tumbek
and
A.
Winkler
,
Surf. Sci.
606
,
L55
(
2012
).
49.
L.
Yu
,
D.
Qian
,
S.
Marina
,
F. A. A.
Nugroho
,
A.
Sharma
,
S.
Hultmark
,
A. I.
Hofmann
,
R.
Kroon
,
J.
Benduhn
,
D.-M.
Smilgies
,
K.
Vandewal
,
M. R.
Andersson
,
C.
Langhammer
,
J.
Martín
,
F.
Gao
, and
C.
Müller
,
ACS Appl. Mater. Interfaces
11
,
21766
(
2019
).
50.
J. S.
Bangsund
,
T. R.
Fielitz
,
T. J.
Steiner
,
K.
Shi
,
J. R.
Van Sambeek
,
C. P.
Clark
, and
R. J.
Holmes
,
Nat. Mater.
18
,
725
(
2019
).
51.
D.
Turnbull
and
J. C.
Fisher
,
J. Chem. Phys.
17
,
71
(
1949
).
52.
V.
Shneidman
,
Sov. Phys. Tech. Phys.
37
,
76
(
1987
).
53.
V.
Shneidman
,
Sov. Phys. Tech. Phys.
33
,
1338
(
1988
).
54.
G.
Shi
,
J. H.
Seinfeld
, and
K.
Okuyama
,
Phys. Rev. A
41
,
2101
(
1990
).
55.
A. H.
Nayfeh
,
Introduction to Perturbation Techniques
(
John Wiley & Sons
,
2011
).
56.
NIST Handbook of Mathematical Functions
, edited by
F. W.
Olver
,
D. W.
Lozier
,
R. F.
Boisverrt
, and
C. W.
Clark
(
Cambridge University Press
,
New York
,
2010
).
57.
J.
Hoyt
and
G.
Sundar
,
Scr. Metall. Mater.
29
,
1535
(
1993
).
58.
V. A.
Shneidman
and
M. C.
Weinberg
,
J. Chem. Phys.
97
,
3621
(
1992
).
59.
V.
Aksenov
,
T.
Tropin
,
M. V.
Avdeev
,
V.
Priezzhev
, and
J.
Schmelzer
,
Phys. Part. Nucl.
36
,
S52
(
2005
).
60.
S. B. P. E.
Timmermans
,
A.
Ramezani
,
T.
Montalvo
,
M.
Nguyen
,
P.
van der Schoot
,
J. C. M.
van Hest
, and
R.
Zandi
,
J. Am. Chem. Soc.
144
,
12608
(
2022
).
61.

The induction time is not to be confused with the mean-first-passage-time, which is sometimes also referred to as the induction time. It describes the time delay for the first cluster to emerge, and is inversely proportional to the product of the volume of the system under investigation and the steady-state nucleation rate.82 

62.
C.
Spinella
,
S.
Lombardo
, and
F.
Priolo
,
J. Appl. Phys.
84
,
5383
(
1998
).
63.
64.
J.
Wedekind
,
G.
Chkonia
,
J.
Wölk
,
R.
Strey
, and
D.
Reguera
,
J. Chem. Phys.
131
,
114506
(
2009
).
66.
L. S.
Bartell
and
D. T.
Wu
,
J. Chem. Phys.
127
,
174507
(
2007
).
67.
G. C.
Sosso
,
J.
Chen
,
S. J.
Cox
,
M.
Fitzner
,
P.
Pedevilla
,
A.
Zen
, and
A.
Michaelides
,
Chem. Rev.
116
,
7078
(
2016
).
68.
E. H. A.
de Hoog
and
H. N. W.
Lekkerkerker
,
J. Phys. Chem. B
103
,
5274
(
1999
).
69.
P.
van der Schoot
,
J. Phys. Chem. B
103
,
8804
(
1999
).
70.
N.
Asherie
,
A.
Lomakin
, and
G. B.
Benedek
,
Phys. Rev. Lett.
77
,
4832
(
1996
).
71.
M. H. J.
Hagen
,
E. J.
Meijer
,
G. C. A. M.
Mooij
,
D.
Frenkel
, and
H. N. W.
Lekkerkerker
,
Nature
365
,
425
(
1993
).
73.
R.
Becker
and
W.
Döring
,
Ann. Phys.
416
,
719
(
1935
).
74.
D. T.
Wu
,
J. Chem. Phys.
97
,
1922
(
1992
).
75.
H.
Hu
,
M.
Singh
,
X.
Wan
,
J.
Tang
,
C.-W.
Chu
, and
G.
Li
,
J. Mater. Chem. A
8
,
1578
(
2020
).
76.
L.
Zeng
,
S.
Chen
,
K.
Forberich
,
C. J.
Brabec
,
Y.
Mai
, and
F.
Guo
,
Energy Environ. Sci.
13
,
4666
(
2020
).
77.
W. J. E. M.
Habraken
,
J.
Tao
,
L. J.
Brylka
,
H.
Friedrich
,
L.
Bertinetti
,
A. S.
Schenk
,
A.
Verch
,
V.
Dmitrovic
,
P. H. H.
Bomans
,
P. M.
Frederik
,
J.
Laven
,
P.
van der Schoot
,
B.
Aichmayer
,
G.
de With
,
J. J.
DeYoreo
, and
N. A. J. M.
Sommerdijk
,
Nat. Commun.
4
,
1507
(
2013
).
78.
A.
Levitsky
,
S. A.
Schneider
,
E.
Rabkin
,
M. F.
Toney
, and
G. L.
Frey
,
Mater. Horiz.
8
,
1272
(
2021
).
79.
D.
Turnbull
and
B.
Vonnegut
,
Ind. Eng. Chem.
44
,
1292
(
1952
).
80.
K.
Kelton
and
A.
Greer
,
Nucleation in Condensed Matter: Application in Materials and Biology
, 1st ed. (
Pergamon
,
Amsterdam
,
2010
).
81.
Y. S.
Djikaev
,
A.
Tabazadeh
, and
H.
Reiss
,
J. Chem. Phys.
118
,
6572
(
2003
).
82.
J.
Wedekind
,
R.
Strey
, and
D.
Reguera
,
J. Chem. Phys.
126
,
134103
(
2007
).

Supplementary Material