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.

## I. INTRODUCTION

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 experimentally^{9–11} but also by means of molecular simulations^{9,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 generic^{7,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 Shneidman^{22,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.

## II. CLASSICAL NUCLEATION THEORY

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.

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.

*ρ*(

*n*,

*t*) for

*n*> 1 to be described by the continuity equation in

*n*-space,

^{21,32}

*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

*v*

_{0}. (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.

*J*(

*n*,

*t*) obeys

^{21,32}

*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/

*k*

_{B}

*T*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}

*W*(

*n*,

*t*) to form an

*n*-mer from

*n*monomers reads

^{35}

^{,}

*μ*(

*t*) a difference in chemical potential between the condensed and dissolved states,

*γ*

_{∞}the interfacial tension of a flat interface, and

*A*

_{0}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,

^{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.

*W*are

^{35}(i) the maximum height

^{21}All three properties depend on time on account of the ongoing process of the evaporation of the solvent. Here, we define $\epsilon \u22613n*\u22121/3\sigma \u22121/2$ for reasons of convenience, with

*σ*≡

*γ*

_{∞}

*A*

_{0}a dimensionless interfacial tension, recalling that

*A*

_{0}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.

*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*) =

*v*

_{0}

*N*(

*n*)/

*V*(

*t*), with

*v*

_{0}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

*χ*(

*t*) as

*χ*(

*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 *k*_{B}*T*, 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*) × *t*_{exp} = 1. Here, *J* is the nucleation flux, implicitly dependent on the supersaturation *S*(*t*), *V*(*t*) is the volume of the solution, and *t*_{exp} 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.

*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) reads

^{21}

*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

*F*

_{eq}(

*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.

## III. NON-STATIONARY NUCLEATION THEORY

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.

*n*at time

*t*as

*F*

_{eq}(

*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

*χ*(

*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

*monomers*in the two phases, Δ

*μ*(

*t*).

*ν*(

*n*,

*t*) as a scaled probability distribution function,

*∂*Δ

*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).

*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

*t*

_{crit}(

*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 $\u2202/\u2202n\u2192n*\u22121\u2202/\u2202x$ and to $\u2202/\u2202tn\u2192\u2202/\u2202tx+(\u2202x/\u2202t)n\u2202/\u2202x$, 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,

*t*

_{crit}(

*t*). To do so, we introduce in Eq. (15) the reduced attachment frequencies $k\u0304+(x,t)=k+(x,t)/k+(1,t)$ with $k+(1,t)\u2261k+*(t)$, the critical attachment frequency evaluated at the critical cluster size

*n*=

*n*

_{*}(or

*x*= 1), and make use of Eq. (8) to find

*x*,

*t*)-coordinate system

*ɛ*is a measure for the reciprocal of the nucleation barrier height [see Eq. (6)], which corresponds to

*x*= 1 in Eq. (19).

*t*

_{crit}(

*t*) as

^{21,26,41,42}

*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}

*τ*, making use of the lifetime of a critical nucleus

*t*

_{crit}(

*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,

*t*

_{crit}(

*t*). Hence, we conclude that the convenient time variable would be one that yields d

*τ*= d

*t*/

*t*

_{crit}(

*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

*τ*reads

^{44}

*t*

_{crit}and simply takes into account that the unit of time drifts with time, emphasizing the highly non-linear character of the problem at hand.

*ζ*(

*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

*ω*(

*τ*) =

*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 $\psi (x,\tau )=(x\u2212n*\u22121(\tau ))$ is the cluster size-dependent part of

*ζ*(

*x*,

*τ*) and follows straightforwardly from Eq. (14).

*t*as

*ω*(

*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 *F*_{eq}. 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}

*W*

_{*}≫ 1, by itself a prerequisite for applying CNT, and (ii) solvent evaporation is sufficiently slow so that the following conditions hold:

*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\u0304+(x,\tau )\u2261k\u0304+(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 $\omega <O(\epsilon \u22122)$. We justify our treatment of

*ζ*as a constant

*a*

*posteriori*.

With these simplifications, we may neglect for the scaled cluster sizes $x\u226a(\epsilon 2\omega )\u22121$ the contribution from the second term in Eq. (18), reading *x*(*∂* ln *n*_{*}/*∂τ*) ∼ −*xɛ*^{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<(\epsilon 2\omega )\u22121$ 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}

*ɛ*,

*n*

_{*}, and

*ω*to be constant and

*ζ*(

*x*) and $k\u0304+(x)$ to be independent of time, Eq. (17) reduces to

*g*(

*x*), is time-independent in the current coordinate system and given by

*ɛ*

^{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.

*τ*= 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*,

*τ*) =

*ρ*(

*x*,

*τ*)/

*F*

_{eq}(

*x*,

*τ*), with

*ρ*(

*x*,

*τ*) the cluster number densities and

*F*

_{eq}(

*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.

^{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 as

^{21,51}

*ρ*

_{1}(

*t*) is the dimensionless monomer density, and

*k*

_{0}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}

^{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

*α*, 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 $\alpha \u22121=hdh/dt\u22121$. Experimentally, drying times for micrometer-thick films can be mere seconds.

^{16}

*ω*, 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*) 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

*αt*

_{crit}(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*/

*t*

_{crit}(0) =

*τ*(1 −

*Bτ*/

*τ*

_{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.

*ω*(0) =

*αt*

_{crit}(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

*ω*(0) implicitly depends on the initial barrier height via

*n*

_{*}and

*t*

_{crit}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.

## IV. NUCLEATION FLUX

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 $\omega *=(1\u2212n*\u22121)\omega =10\u22121$ 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.

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}

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 Shneidman^{26} 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.

*x*> 1, in the context of homogeneous unary nucleation in a solution that becomes increasingly concentrated as a result of solvent evaporation, reads

*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

*ω*the non-stationary index defined in Eq. (24). The quasi-steady state flux

*J*

_{QSS}is defined as

*J*

_{SS}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.

*τ*

_{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

*g*(

*x*′) is the growth rate defined in Eq. (27), and we define

*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.

*τ*

_{cor}that emerges is given by

*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 $\tau \u2212\tau i(x)+\tau cor+\omega *\u22121ln\omega *\Gamma (\omega *)$, 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 *J*_{0}(*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, *J*_{tQSS}(*x*, *τ*) = *J*_{QSS}(0) = *J*_{SS}. 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 *J*_{0}(*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, *J*_{0}(*x*) can no longer be justifiably set equal to zero. Remarkably, accounting for a non-zero *J*_{0}(*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.

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 × 10^{5}, 1.1 × 10^{4}, 1.1 × 10^{3}, and 5.6 × 10^{1}, respectively. We normalize this quantity (left) by the constant quasi steady-state flux *J*_{QSS}(*τ* = 0) and (right) the time-dependent renormalized and time-shifted version of it, *J*_{tQSS}(*τ*), given in Eq. (36). As both *J*_{tQSS}(*τ*) 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 $\tau s\u22121\u2248\epsilon 2\omega *\u2248(\Delta W*)\u22121\omega *\u226a1$. For the largest value of *ω*_{*} shown, the time scale *τ*_{s} = 1.1 × 10^{1} 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 *J*_{QSS}(0) = *J*_{SS}. 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), lim_{z→0} *Q*(*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.

*ω*

_{*}→ 0, Eq. (35) reduces to a transient time-shifted quasi-steady-state nucleation flux, realizing that $Q(\omega *+1,z)\u223cexp(\u2212z)1+O(\omega *)$,

^{56}yielding

*J*

_{QSS}(

*τ*−

*τ*

_{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,

*J*

_{QSS}(

*τ*−

*τ*

_{i}(

*x*) +

*τ*

_{cor}) =

*J*

_{SS}, 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 *J*_{QSS}(*τ*) and the time-shifted QSS flux obtained from our theory *J*_{QSS}(*τ* − *τ*_{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 *k*_{B}*T*, 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 *J*_{QSS} 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}

## V. INDUCTION TIME

*ω*. 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}

*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,

*τ*→

*∞*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

*J*

_{0}(

*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

*J*

_{tQSS}(

*x*,

*τ*′) asymptotically for

*τ*→

*∞*; see also the discussion near Eq. (35).

*J*

_{0}(

*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

*J*

_{QSS}(0) =

*J*

_{SS}given by Eq. (11) and E

_{p}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

*τ*′ yields

*z*

^{−1}Γ(

*z*+ 1) = Γ(

*z*). Combined with Eqs. (45) and (46), and after a reordering of terms, we find

*τ*.

*τ*→

*∞*. We conveniently introduce the series expansion for the exponential integral as

^{56}

*ω*

_{*}≥ 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

*ω*

_{*},

*ω*

_{*}) into the time shift. Hence, it appears to be a fundamental, non-stationary contribution to the induction or lag time in non-stationary nucleation.

*θ*

_{0}(

*x*) ≡

*θ*(

*x*,

*ω*

_{*}= 0) =

*τ*

_{i}(

*x*) +

*γ*

_{E}≥ 1,

^{58}with

*γ*

_{E}= 0.577… the Euler–Mascheroni constant,

^{56}as

*ω*

_{*}→ 0, which can be verified using the expansion

^{56}

*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

*ω*

_{*},

*ω*

_{*}→ 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.

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 *J*_{0}(*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 × 10^{4}. 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.

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.

## VI. LIMITATIONS OF THE THEORY

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 *t*_{crit}(*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 proteins^{70} and is suspected for dispersions of fullerenes under appropriate conditions.^{71} This scenario is explicitly excluded from our analysis.

*ω*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

*ω*> 0. Equation (56) shows that our predictions should indeed be valid if

*τ*

_{s}≫ 1. Here,

*τ*

_{s}= 2/3 ×

*ɛ*

^{2}

*ω*(0) with

*ω*(0) =

*αt*

_{crit}(0)

*n*

_{*}(0), with

*ɛ*the reciprocal square root of the barrier height,

*α*the reciprocal drying time of the thin film,

*t*

_{crit}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, $\omega (0)\u223c\Delta W*5/2$, whereas for a diffusion-limited aggregation model, $\omega (0)\u223c\Delta 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 $\omega *=(1\u2212n*\u22121)\omega $ 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^{−2} *S*(0). Nucleation barriers for homogeneous nucleation observable on experimental time scales are generally estimated between 10 and 100 *k*_{B}*T*.^{21} For higher nucleation barriers, nucleation almost exclusively commences via heterogeneous nucleation.^{21} Hence, let us take 100 *k*_{B}*T* 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 *t*_{crit} scales linearly with the barrier height, *t*_{crit} ∼ Δ*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 $\omega *\u223c\alpha \Delta 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 *k*_{B}*T* at a fixed evaporation rate *α*, *ω*_{*} decreases by a factor of about 3 × 10^{2}.

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.

^{21,73}

*n*,

*J*

_{n−1,n}the cluster flux density in

*n*-space, and

*K*

_{n}(

*t*) again the evaporation source term defined by Eq. (9). The cluster flux density is given by

^{21}

*n*-mer, and

*F*

_{eq}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

*n*

_{max}= 10 000 ≫

*n*

_{*}as the maximum cluster size. This value for

*n*

_{max}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 $\omega *\u223c\Delta 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.

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.

## VII. DISCUSSION AND CONCLUSION

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-stationarity^{26} 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 simulations^{8} 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.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

## NOMENCLATURE

*A*area (m

^{2})*F*_{eq}dimensionless “equilibrium” cluster distribution

*g*(deterministic) cluster growth rate

*J*nucleation flux (s

^{−1})*J*_{QSS}quasi-steady-state nucleation flux

*J*_{SS}steady-state nucleation flux

*J*_{tQSS}renormalized, time-shifted quasi-steady-state nucleation flux

*K*evaporation source nucleation master equation (s

^{−1})*k*^{+}attachment frequencies (s

^{−1})- $k\u0304+$
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/(*K*_{B}*T*)reciprocal thermal energy (

*J*^{−1})*S*supersaturation

*t*_{crit}lifetime of a cluster (s)

*V*volume in which the solution is contained (m

^{3})*v*_{0}molecular volume of solute molecules (m

^{3})*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

*ν*=*ρ*/*F*_{eq}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 $\omega *=\omega (1\u2212n*\u22121)$

## NOMENCLATURE

*A*area (m

^{2})*F*_{eq}dimensionless “equilibrium” cluster distribution

*g*(deterministic) cluster growth rate

*J*nucleation flux (s

^{−1})*J*_{QSS}quasi-steady-state nucleation flux

*J*_{SS}steady-state nucleation flux

*J*_{tQSS}renormalized, time-shifted quasi-steady-state nucleation flux

*K*evaporation source nucleation master equation (s

^{−1})*k*^{+}attachment frequencies (s

^{−1})- $k\u0304+$
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/(*K*_{B}*T*)reciprocal thermal energy (

*J*^{−1})*S*supersaturation

*t*_{crit}lifetime of a cluster (s)

*V*volume in which the solution is contained (m

^{3})*v*_{0}molecular volume of solute molecules (m

^{3})*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

*ν*=*ρ*/*F*_{eq}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 $\omega *=\omega (1\u2212n*\u22121)$

### APPENDIX: HETEROGENEOUS NUCLEATION

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.

^{21}

^{,}

*σ*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 lower

^{79,80}or increase

^{81}the nucleation barrier with respect to homogeneous nucleation. We refer to Refs. 21, 79, and 80 for details on the scaling function

*f*(

*ϕ*).

*μ*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

*ρ*

_{0}is the

*surface*density of nucleation sites, which for a flat interface is constant in time as evaporation does not affect this quantity.

*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.

*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

*ω*, 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

*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*.

## REFERENCES

*Nucleation Theory, Lecture Notes in Physics*

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

*Solid State Physics*

*Fokker–Planck Equation*

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 ¯*k**k*^{+}(*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}

*NIST Handbook of Mathematical Functions*

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}

*Nucleation in Condensed Matter: Application in Materials and Biology*