molecules, which seems to be accurate in practice Bornside  _et al._ (1989).
If the initial solute concentration is low and we limit ourselves to the
initial stages of the drying process, this then predicts the volume of our
solution to decrease linearly with time. The resulting monomer density is then
given by

$\rho_{1}(t)/\rho_{1}(0)=V(0)/V(t)\equiv\chi(t)=\frac{1}{1-\alpha t}.$ (76)

Here, $\alpha$ is the reciprocal of the drying time, a measure for the rate of
evaporation. We treat $\alpha$ as a freely adjustable parameter and maintain
our assumption that monomer depletion due to aggregation is negligible. This
simple description of the evaporation kinetics applies, e.g., in the context
of spin-coating of solutions in the initial stages of the drying process van
Franeker  _et al._ (2015a); Bornside  _et al._ (1989), where the drying time
can be expressed in terms of the height $h$ of the film
$\alpha^{-1}=h\left(\mathrm{d}{h}/\mathrm{d}{t}\right)^{-1}$. Experimentally,
drying times for micrometer-thick films can be mere seconds van Franeker  _et
al._ (2015c).

Using these models for monomer attachment and solvent evaporation, we are now
able to investigate for what evaporation rates our analytical result Eq. (60)
remains valid, and under what conditions the constraints Eq. (26) hold. For
this, we first obtain an explicit expression for the scaled time $\tau(t)$ as
a function of the actual, physical time $t$ by invoking Eqs. (22) and (17),
and subsequently calculate the non-stationarity index $\omega$, which we no
longer assume to be constant in time. Using Eq. (75), we find the life time of
a critical cluster Eq. (17) to obey the relation

$\displaystyle\frac{t_{\mathrm{crit}}(t)}{t_{\mathrm{crit}}(0)}$
$\displaystyle=\chi^{-1}(t)\left(1+\Delta\mu_{0}^{-1}\ln\chi(t)\right)^{-2}.$
(77)

scaled to the value at the initial time $t=0$. Here,
$\Delta\mu_{0}=\ln\left[\rho(t=0)/\rho_{\infty}\right]=\ln S(t=0)$ is the
chemical potential difference of a monomer in the dissolved and the condensed
phase at $t=0$.

Next, we insert this expression in the definition of the scaled time Eq. (22),
and find

$\tau(t)=\frac{1}{t_{\mathrm{crit}}(t=0)}\int_{0}^{t}\mathrm{d}{t}^{\prime}\chi(t^{\prime})(1+\Delta\mu_{0}^{-1}\ln\chi(t^{\prime}))^{2},$
(78)

which, combined with Eq. (76), yields

$\tau(t)=\tau_{\mathrm{s}}\left[(1+\Delta\mu^{-1}_{0}\ln\chi(t))^{3}-1\right],$
(79)

expressed in terms of the reciprocal volume $\chi(t)$ defined by Eq. (76) for
notational simplicity. Here,
$\tau_{\mathrm{s}}={\Delta\mu_{0}}/{3\alpha{t_{\mathrm{crit}}}(0)}$ is a time
scale that we interpret below. Unfortunately, inversion of Eq. (79) does not
produce a simple expression of the actual time $t$ in terms of the scaled time
$\tau$. For early times, a Taylor expansion yields
$t/t_{\mathrm{crit}}(0)=\tau_{\mathrm{s}}/2\sqrt{B}\left(\sqrt{\tau/\tau_{\mathrm{s}}+1/4B}-1/2\sqrt{B}\right)$,
with a positively valued $B=(2+\Delta\mu_{0})/6$, showing that the actual time
$t$ increases sub-linearly with the scaled time $\tau$. This is expected
because solvent evaporation speeds up the nucleation process and the time
$\tau$ measures the time that passes from the point of view of the nucleating
solute particles.

We next set out to find the time dependence of the non-stationarity index
$\omega$ defined in Eq. (23), see also the list of symbols List of Symbols.
Here, we focus on $\omega$ instead of $\omega_{*}$, as the former is the
relevant quantity with respect to the constraints in Eq. (26). The two
quantities $\omega$ and $\omega_{*}$ are related via the critical cluster size
as $\omega_{*}=\omega(1-n_{*}^{-1})$, see Eq. (37). Although we have
considered this quantity to be so weakly dependent on time that we can
consider it a constant in our calculations, it can, by its definition Eq.
(23), only be actually independent of time in special cases. Indeed, we find
from the above analysis that

$\omega(t)\equiv
t_{\mathrm{crit}}n_{*}\frac{\mathrm{\partial}{\ln\chi(t)}}{\mathrm{\partial}{t}}=\omega(0)\left(1+\Delta\mu_{0}^{-1}\ln\chi(t)\right)^{-5}=\omega(0)\left(\frac{\tau(t)}{\tau_{\mathrm{s}}}+1\right)^{-5/3},$
(80)

where, $\omega(0)=\alpha{t_{\mathrm{crit}}}(0)n_{*}(0)$. For high nucleation
barriers $\beta\Delta W_{*}$ and sufficiently large critical cluster sizes
$n_{*}$, we can express the time scale $\tau_{\mathrm{s}}$ as the ratio of the
height of the dimensionless nucleation barrier and the non-stationarity index
$\omega$, both evaluated at $\tau=0$, as

$\tau_{\mathrm{s}}={2\beta\Delta
W_{*}(0)}/{3\omega(0)}\approx{2}/{3\varepsilon^{2}\omega(0)}\sim\mathcal{O}\left(\varepsilon^{3}\right),$
(81)

where we neglect a correction term of order $n_{*}^{-2/3}$, use the expression
for the nucleation barrier Eq. (6), and make explicit that $\omega(0)$
implicitly depends on the initial barrier height via $n_{*}$ and
$t_{\mathrm{crit}}$ as $\omega(0)\sim\varepsilon^{-5}$, see Eqs. (7), (17) and
(23).

We are now in a position to interpret the dimensionless timescale
$\tau_{\mathrm{s}}$ as the timescale over which $\omega$ changes with time,
which turns out to scale with the dimensionless nucleation barrier height as
$\tau_{\mathrm{s}}\sim\beta\Delta W_{*}^{-3/2}$. As expressed in Eq. (26), our
analytical results apply if the temporal drift in $\omega$ is small. Combining
the constraint on $\omega$ of Eq. (26) with Eq. (80), we find

$-\frac{1}{\omega}\frac{\mathrm{d}{\omega}}{\mathrm{d}{\tau}}=\frac{5}{3}\left(\frac{1}{\tau_{\mathrm{s}}+\tau}\right)\ll
1,$ (82)

indicating that the non-stationarity index decreases with time, since for
solvent evaporation $\omega>0$. Eq. (82) shows that our predictions should be
valid if $\tau_{\mathrm{s}}\gg 1$. That $\tau_{\mathrm{s}}$ must be large is
actually implicit in Eq. (26), as from Eq. (81) we have $\tau_{\mathrm{s}}\sim
1/\varepsilon^{2}\omega$, which must be (much) larger than unity for our
results to be internally consistent. Note that this also means that the
constraints presented in Eq. (26) can be written as a single constraint:
$\tau_{\mathrm{s}}\gg 1$.

Based on Eq. (82) we might conclude that our results should also apply for
$\tau\gg 1$ even if $\tau_{\mathrm{s}}$ is not much larger than unity. Because
we presume a fixed value for $\omega_{*}\equiv\omega_{*}(0)$ for our
theoretical analysis, this in actual fact is not the case. Indeed, even if
$\omega_{*}$ only drifts very weakly with time, the relative errors between
the actual nucleation flux and our theoretical nucleation flux must accumulate
with time. Eq. (82) actually suggests a natural way to correct for this: we
could a posteriori make the non-stationarity index $\omega_{*}$, the critical
cluster size $n_{*}$ and $\varepsilon$ 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 $\tau-\tau_{\mathrm{i}}(x)$, with
$\tau_{\mathrm{i}}(x)$ the incubation time, as was in fact suggested by
Shneidman Shneidman (2010). As a result of this approach we deduce from Eq.
(82) that our results are expected to become accurate if $\tau\gg 1$ also.

The question now arises how strongly solvent evaporation actually affects
nucleation under experimentally relevant conditions. As discussed in Sec. II,
this depends on the conditions present at the moment that the first
(supercritical) nuclei form and is hence related to our final assumption
(iii): how the choice of the initial supersaturation affects our results.
Considering that the induction time appears to vanish very close to the
binodal, which is a consequence of the non-stationarity index $\omega_{*}$
becoming large very close to the binodal as discussed in the previous section,
we only need to consider the effect of the initial conditions on the tQSS
flux. To study how the choice for the initial conditions affect our results
requires a reasonable estimate for the initial nucleation barrier height and a
reasonable good understanding of the nucleation kinetics, so that we can
estimate the attachment frequencies. These attachment frequencies can be
determined by combining Eq. (75) with parameter estimates from either
experimental or atomistic numerical studies Kelton (1991). 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$, and, consequently, neither the magnitude of the evaporation-related
effects on nucleation under experimentally relevant conditions.

In order to nevertheless make headway, we estimate a reasonable value for the
smallest initial supersaturation required and subsequently discuss how
changing this initial supersaturation changes our results qualitatively.
Nucleation barriers for homogeneous nucleation observable on experimental time
scales are generally estimated between 10 - 100 $k_{\mathrm{B}}T$ Kashchiev
(2000). For higher nucleation barriers, nucleation almost exclusively
commences via heterogeneous nucleation Kashchiev (2000). Hence, let us take
100 $k_{\mathrm{B}}T$ as a reasonable estimate for the maximum barrier that
results in nucleation on experimental timescales and we use a value for the
dimensionless surface tension, see Eq. (6), $\sigma\approx 8.0\pm 0.4$,
corresponding with the nucleation of fullerene clusters in carbon disulfide
Aksenov  _et al._ (2005). This results in a maximum value for the initial
supersaturation of $S(0)\approx 2.18-2.47$ with a maximum critical cluster
size estimate of $n_{*}\approx 239-275$.

As long as the constraint $\tau_{\mathrm{s}}\gg 1$ presented in Eq. (81) is
valid for the conditions corresponding to the initial supersaturation, the
expressions for the nucleation rate Eqs. (60) and (61) and the result for the
induction time in Eq. (72) are expected to be accurate. Upon changing the
initial supersaturation for fixed values of the other kinetic and
thermodynamic parameters, we find that 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 scales linearly with the barrier height,
$t_{\mathrm{crit}}\sim\beta\Delta W_{*}$, indicating that if we were to double
the initial barrier height, the initial nucleation time scale would also
double. Since all results are presented in terms of the scaled time $\tau$
instead of the actual time $t$, this change is actually implicit in the
results presented. Second, the non-stationarity index $\omega_{*}$ depends on
the initial nucleation barrier as $\omega_{*}\sim\beta\Delta W_{*}^{5/2}$, so
for fixed volumetric evaporation rates $\alpha$ at different initial barrier
heights, the non-stationarity index can change considerably. For example, if
the initial barrier height changes from 100 $k_{\mathrm{B}}T$ to 10
$k_{\mathrm{B}}T$, $\omega_{*}$ decreases by a factor of about $3\times
10^{2}$. As a result, the impact of solvent evaporation on nucleation
decreases rapidly (superlinearly) with increasing initial degree of
supersaturation.

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 pre-empted by heterogeneous nucleation. This might occur
if solvent evaporation is very rapid, resulting in accumulation of solute near
the solution-air interface Schaefer  _et al._ (2017). 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 wet film thickness. As we discuss in Appendix B, the extension of our
work to heterogeneous nucleation is straightforward to carry out and all
results in this manuscript apply with only minor modifications. We refer to
Appendix B for the details.

To further investigate the limitations of our theory, we return to numerics,
reminding the reader that we compared our theory with numerical calculations
in Fig. 4, without explaining how we set up the latter. This we do now, and
subsequently re-examine the results presented in Fig. 4, returning in detail
to the origins of the decreasing accuracy for the highest values of
$\omega_{*}$ in Fig. 4. Finally, as we found in the previous section that the
induction time vanishes at very small supersaturation, we expect that only the
tQSS nucleation flux of Eq. (61) is relevant for experimental situations.
Hence, we next focus attention on the nucleation flux at times much later than
the induction time and discuss the differences between the tQSS and the QSS
fluxes.

For the numerical calculations we solve the discrete version of our nucleation
equations, which are numerically more simple to solve. The discrete version of
the nucleation Master equation Eq. (1) reads

$\frac{\mathrm{\partial}{\rho_{n}(t)}}{\mathrm{\partial}{t}}=k^{+}_{n-1,n}(t)\rho_{n-1}(t)-\left(k^{-}_{n,n-1}(t)+k^{+}_{n,n+1}(t)\right)\rho_{n}(t)+k^{-}_{n+1,n}(t)\rho_{n+1}(t)+K_{n}(t).$
(83)

Here, the subscripts indicate the cluster size, with $k^{+}_{n,n+1}(t)$ the
attachment frequency for a monomer to an $n$-mer, $k^{-}_{n,n-1}(t)$ the
detachment frequency of a monomer from an $n$-mer and $K_{n}(t)$ the
evaporation source term given by Eq. (9). We associate the detachment
frequencies to the attachment frequencies using the same detailed balance
argument we used for the continuum nucleation equation. See the discussion of
Eq. (2) in Section II. To apply these within the discrete description, we
write

$k^{+}_{n,n+1}(t)F_{\mathrm{eq}}(n,t)=k^{-}_{n+1,n}(t)F_{\mathrm{eq}}(n+1,t),$
(84)

with $F_{\mathrm{eq}}$ again the “equilibrium” cluster densities (3). We solve
the system of equations using a standard backward Euler method with time step
size $\Delta\tau=0.05$, assuming $n_{\mathrm{max}}=10000\gg n_{*}$ as the
maximum cluster size. This value for $n_{\mathrm{max}}$ is sufficiently large
not to affect our results for the nucleation flux.

In order to solve Eq. (83), we need to set appropriate initial conditions.
Here, we again use the arbitrarily chosen nucleation barrier of $\beta\Delta
W_{*}(0)=70$ and critical cluster size of $n_{*}(0)=152$, and compare for
different values for the non-stationarity index $\omega_{*}$. These values for
$\beta\Delta 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 in principle be accurate. As long as the constraint
presented in Eq. (82) 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_{*}\sim\beta\Delta W_{*}^{5/2}$ at fixed
evaporation rate, as discussed above.

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 the discrete nucleation equations Shneidman and Weinberg
(1992); Wu (1992b). Specifically, they may result in slightly different
induction times. The magnitude of these differences turn 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 Shneidman and Weinberg (1992); Shneidman (2010).

We can now re-examine the theoretical and numerical results shown in Fig. 4.
As discussed in Section. VI, the theoretical nucleation fluxes shown in the
left panel of Fig. 4 are in very good agreement with the numerical results for
the values for $\omega_{*}$ shown, although the agreement decreases somewhat
for the larger values of $\omega_{*}$. We can now explain this by considering
the time scale $\tau_{\mathrm{s}}$. For the smallest values of $\omega_{*}$
shown ($\omega_{*}=5\times 10^{-3}$), the time scale
$\tau_{\mathrm{s}}=1.1\times 10^{5}$ is very large, meaning that the
assumptions underlying our analytical results are warranted. For the largest
value of $\omega_{*}$ shown, the time scale $\tau_{\mathrm{s}}=1.1\times
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 $\tau>12$.

Again, 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 question might arise how much more accurate
the tQSS flux is compared to the QSS flux. This 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 $\beta\Delta W(0)=70$ and
compare them with the QSS flux Eq. (57) (dotted, color-matched), as well as
the tQSS flux of Eq. (61) (dashed, color matched). As discussed above, we let
$\omega_{*}$, $n_{*}$ and $\varepsilon$ depend on the shifted time
$\tau-\tau_{\mathrm{i}}(x)$ to account for the weak drift in these quantities.
In Fig. 7 we do not include the initial time for $\tau<8$, but focus on the
theoretical expressions that are accurate at later times. For these later
times, we find that for low values of $\omega_{*}$ (red and green curves), the
QSS and tQSS fluxes are, for all intends and purposes, identical. This we also
expect based on the small $\omega_{*}$ limit of Eq. (61), see Eq. (62). At
$\tau=12$ we are (far) beyond the induction time, and find that at this
instant in time for higher values of $\omega_{*}$, the QSS flux overpredicts
the nucleation flux by about 70 % for $\omega_{*}=0.5$ (orange), and by about
90 % for $\omega_{*}=1.0$ (blue). The tQSS flux also overpredicts the
nucleation flux, but does so to a significantly smaller extent, about 10 % for
$\omega_{*}=0.5$, and about 20 % for $\omega_{*}=1.0$.

Figure 7: The nucleation flux $J(x,\tau)$ evaluated at the scaled cluster size
$x=1.5$, normalized by the $J_{\mathrm{QSS}}(\tau=0)$ (see Eq. (57)) as a
function of the scaled time $\tau$. Crosses are the numerical solutions with
time interval between the data points of $\Delta\tau=1/4$, dotted lines are
the QSS flux Eq. (57), dashed lines the tQSS flux given by Eq. (61). We set
the critical cluster size $n_{*}=152$ and initial nucleation barrier
$\beta\Delta W_{*}(\tau=0)=70$. We show results for $\omega_{*}(0)=5\times
10^{-3}$, $\tau_{\mathrm{s}}=1.1\times 10^{5}$ (red), $\omega_{*}(0)=5\times
10^{-2}$, $\tau_{\mathrm{s}}=1.1\times 10^{4}$ (green), $\omega_{*}(0)=5\times
10^{-1}$, $\tau_{\mathrm{s}}=1.1\times 10^{3}$ (orange) and
$\omega_{*}(0)=1.12$, $\tau_{\mathrm{s}}=5.6\times 10^{1}$ (blue).

Clearly, the tQSS nucleation flux (taking into account the time-drifting
$\omega_{*}$, $n_{*}$ and $\varepsilon$) 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 $\tau_{\mathrm{s}}$, see Eq. (82). The temporal
drift in the non-stationarity index $\omega_{*}$, the critical cluster size
$n_{*}$ and in $\varepsilon$ can in this case not be accounted for a
posteriori. For $\tau_{\mathrm{s}}$ small, the temporal drift should actually
be taken into account in the calculations presented in Sec. III, V and VI.
This significantly complicates finding analytical solutions to the nucleation
Master equation, and we are at present not aware of any method that allows one
to achieve this.

##  IX 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 for an instantaneous quench in the metastable
region upon which nucleation commences under the conditions of on-going
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 Kashchiev (2000). 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 and 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 Kashchiev (2000). 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 Kashchiev (2000). 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 factor 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 Shneidman
(2010) and should therefore have consequences in other fields where nucleation
is non-stationary.

If the non-stationarity index is large, then solvent-evaporation strongly
affects nucleation. The nucleation cascade cannot 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 increasing evaporation rate. This agrees with our
expectation that an increasing supersaturation due to solvent evaporation must
result in continuously 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 the 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, 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. The reason is that the former is inversely
proportional to the non-stationarity index.

Experimentally, producing an instantaneous, non-zero supersaturation upon
which the solvent evaporation commences, seems difficult to achieve. Indeed,
the inspiration of 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 Mei  _et al._ (2013), where the solution
quickly but gradually changes from undersaturated to supersaturated Mei  _et
al._ (2013). For early times just after crossing the binodal, nuclei arguably
do not emerge as nucleation is exceedingly slow 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 co-
ordinate (the cluster size) scales that in general both also depend on time,
and do so differently Trinkaus and Yoo (1987). This causes the kinetic
nucleation equations to be exceedingly hard to solve analytically Kashchiev
(2000); Trinkaus and Yoo (1987); Shneidman (1987, 1988). This is also true for
the relevant time and reaction co-ordinate scales within our model. We solved
this problem by making use of a co-ordinate system that is co-moving with
respect to both the reaction co-ordinate and time, following the earlier work
of Shneidman Shneidman (2010). To obtain analytical predictions within the
prescription of a co-moving co-ordinate 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 done to obtain the quasi-
steady-state nucleation flux Kashchiev (2000); Shneidman (2010). 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 and the
non-stationarity index. If this timescale is sufficiently large (on the scale
of the nucleation time) 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 Hu  _et al._ (2020); Zeng  _et al._
(2020). 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 Habraken  _et al._ (2013). For example,
for perovskite thin film solar cells, a higher evaporation rate results in
more and smaller crystalline domains Hu  _et al._ (2020); Zeng  _et al._
(2020), 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 van Franeker  _et al._ (2015c). 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 Michels  _et al._ (2021) or incorporate it in a nucleation-and-
growth model such as the Johnson-Mehl-Avrami-Kolmogorov theory Cahn (1995),
which is often applied in combination with in-situ diagnostics of the
formation of functional polycrystalline thin films Levitsky  _et al._ (2021).
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
polycrystalline final 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 van Franeker  _et al._ (2015c, b). Such additional
components are known to affect the nucleation-and-growth process van Franeker
_et al._ (2015c), 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.

## 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) for helpful discussions.

## Data availability

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

## List of Symbols

Table 1: List of symbols Symbol | Explanation | Equation | Units
---|---|---|---
$\rho$ | dimensionless cluster densities |  | [ - ]
$F_{\mathrm{eq}}$ | dimensionless “equilibrium” cluster distribution |  | [ - ]
$\nu=\rho/F_{\mathrm{eq}}$ | reduced cluster densities | (12) | [ - ]
$v_{0}$ | molecular volume of solute molecules |  | [ m3 ]
$J$ | nucleation flux | (65) | [ s-1 ]
$K$ | evaporation source nucleation master equation |  | [ s-1 ]
$n$ | cluster size reaction coordinate (number of molecules in cluster) | (9) | [ - ]
$n_{*}$ | critical cluster size | (7) | [ - ]
$x$ | reduced cluster size reaction coordinate |  | [ - ]
$\tau$ | reduced time co-ordinate | (22) | [ - ]
$k^{+}$ | attachment frequencies |  | [s-1]
$\bar{k}^{+}$ | dimensionless attachment frequencies |  | [ - ]
$\beta$ | reciprocal thermal energy |  | $[J^{-1}]$
$\Delta W$ | nucleation free energy | (6) | $[J]$
$\Delta\mu$ | chemical potential difference between phases | (5) | $[J]$
$\gamma_{\infty}$ | surface tension |  | $[J]$
$A$ | area |  | [ m2 ]
$\rho_{\infty}$ | solubility limit |  | [ - ]
$Z$ | Zeldovich Factor | (8) | [ - ]
$\sigma$ | dimensionless surface tension |  | [ - ]
$S$ | supersaturation | (5) | [ - ]
$V$ | volume in which the solution is contained |  | [ m3 ]
$N$ | number of solute molecules in the solution |  | [ - ]
$\chi$ | dimensionless reciprocal volume | (9) | [ - ]
$\zeta$ | non-stationarity function | (23) | [ - ]
$\omega$ | unmodified non-stationarity index | (23), (25) | [ - ]
$\omega_{*}$ | non-stationarity index, defined as $\omega_{*}=\omega(1-n_{*}^{-1})$ |  | [ - ]
$g$ | (deterministic) cluster growth rate | (16), (20) | [ - ]
$t_{\mathrm{crit}}$ | lifetime of a cluster | (17) | [ s ]
$X$ | Stretched reduced cluster size |  | [ - ]
$\tau_{\mathrm{d}}$ | subcritical cluster growth time | (33) | [ - ]
$\tau_{\omega}$ | timescale of the non-stationarity effects on the subcritical growth | (34) | [ - ]
$\psi$ | cluster size dependence of $\zeta$ | (23) | [ - ]
$J_{\mathrm{tQSS}}$ | renormalized, time-shifted quasi-steady-state nucleation flux | (61) | [ - ]
$J_{\mathrm{QSS}}$ | quasi-steady-state nucleation flux | (57) | [ - ]
$J_{\mathrm{SS}}$ | steady-state nucleation flux | (57) | [ - ]
$J_{\mathrm{stat}}$ | nucleation flux under stationary conditions | (58) | [ - ]
$\tau_{\mathrm{i}}$ | incubation time | (59) | [ - ]
$\tau_{\mathrm{stat}}$ | part of the incubation time | (42) | [ - ]
$\tau_{\mathrm{non-stat}}$ | non-stationary contribution to the lag time | (45) | [ - ]
$\tau_{\mathrm{cor}}$ | time shift due to no-stationary effects |  | [ - ]
$\tau_{\mathrm{s}}$ | time scale over which $\omega$ drifts with time | (81) | [ - ]

##  Appendix A Matched Asymptotic Expansions

The dimensionless nucleation master equation Eq. (27) we solve subject to the
boundary and initial conditions in Eq. (29). We introduce the Laplace
transform $\hat{\nu}(x,s)$ of $\nu(x,\tau)$ with respect to the time variable,
presuming that $\varepsilon$, $\zeta(x)$ and $\bar{k}^{+}(x)$ do not depend on
time. This yields Eq. (30). Let us treat $\varepsilon\ll 1$ as an
asymptotically small parameter, implying a high nucleation barrier. The method
we use is known as matched asymptotic expansions Nayfeh (2011); Shneidman
(1987, 1988); Shi  _et al._ (1990). As discussed in the main text, we can
distinguish three regions: Two outer regions for $x<1$ and $x>1$, and an inner
region valid near $|x-1|<\mathcal{O}(\varepsilon)$. We seek three piece-wise
solutions to Eq. (30) that are valid in the three regions separately.

We start with the two outer solutions, where we neglect the diffusive
contribution in Eq. (30) that is proportional to $\varepsilon^{2}$. We
introduce the order-by-order expansion in $\varepsilon$ as

$\hat{\nu}(x,s)=\hat{\nu}_{0}(x,s)+\varepsilon\hat{\nu}_{1}(x,s)+\varepsilon^{2}\hat{\nu}_{2}(x,s)+\cdots.$
(85)

Inserting this into Eq. (30) and only retaining the lowest order terms yields

$\left(s+\zeta(x)\right)\hat{\nu}_{0}(x,s)=g(x)\frac{\mathrm{\partial}{\hat{\nu}_{0}(x,s)}}{\mathrm{\partial}{x}}$
(86)

for both outer solutions. The solution of this differential equation is an
exponential, reading

$\hat{\nu}_{0}(x_{1},s)=\hat{\nu}_{0}(x_{0},s)\exp\left[s\int_{x_{0}}^{x_{1}}\frac{\mathrm{d}{x}^{\prime}(s+\zeta(x^{\prime}))}{g(x^{\prime})}\right].$
(87)

Using the boundary conditions given in Eq. (29), we find the leading order
solutions for the two different outer regions to satisfy Eqs. (31) and (32).

Next we focus on the inner solution $\hat{\nu}_{\mathrm{inner}}(x,s)$, valid
in the region where $|x-1|\ll 1$. In this region we cannot neglect the terms
involving $\varepsilon$ as in this region, as we discuss in detail in the main
text in Sec. V. Standard technique is now to introduce a so-called stretched
variable, to magnify the changes in the solution in this inner region of width
of order $\varepsilon$. In our case, the appropriate stretched variable is
$X=(x-1)/\varepsilon$ Shi  _et al._ (1990); Shneidman (1987, 1988). In terms
of the stretched variables, Eq. (30) reads

$\displaystyle(s+\zeta(X))\hat{\nu}_{\mathrm{inner}}(X,s)$ $\displaystyle=$
(88)
$\displaystyle\quad\left[\frac{\mathrm{\partial}{}}{\mathrm{\partial}{X}}\bar{k}^{+}(X)+\bar{k}^{+}(X)\frac{3}{\varepsilon}\left(1-(1+\varepsilon
X)^{-1/3}\right)\right]\frac{\mathrm{\partial}{\hat{\nu}_{\mathrm{inner}}(X,s)}}{\mathrm{\partial}{X}}+$
$\displaystyle\quad\bar{k}^{+}(X)\frac{\mathrm{\partial}^{2}{\hat{\nu}_{\mathrm{inner}}(X,s)}}{\mathrm{\partial}{X}^{2}},$

where we have inserted explicitly the growth rate
$g(X)=\bar{k}^{+}(X)\frac{3}{\varepsilon}\left(1-(1+\varepsilon
X)^{-1/3}\right)$ expressed in stretched variables. Next, we introduce the
order-by-order expansion for $\hat{\nu}_{\mathrm{inner}}(X,s)$, mirroring the
expansion in $\varepsilon$ for the outer solutions. We apply the same
expansion to the reduced attachment frequencies $\bar{k}^{+}(X)$ and the
measure for the non-stationarity $\zeta(X)$. For both, only the zeroth-order
terms survives, if we make use of the reasonable assumption that the
attachment frequencies vary smoothly with the cluster size $n$. The resulting
expression for the inner solution reads

$\displaystyle\frac{\mathrm{\partial}^{2}{\hat{\nu}_{\mathrm{inner,(0)}}(X,s)}}{\mathrm{\partial}{X}^{2}}+2X\frac{\mathrm{\partial}{\hat{\nu}_{\mathrm{inner,(0)}}(X,s)}}{\mathrm{\partial}{X}}-2\left[s+\zeta(X=0)\right]\hat{\nu}_{\mathrm{inner,(0)}}(X,s)=0,$
(89)

where we use

$\left(1-(1+\varepsilon X)^{-1/3}\right)\approx\frac{1}{3}\varepsilon X.$ (90)

The ordinary differential equation (89) is a special form of the differential
equation for hypergeometric functions Olver  _et al._ (2010), which for some
values for $m=s+\zeta(X=0)$ may yield analytical solutions. In our case, the
only relevant solution reads

$\hat{\nu}_{\mathrm{inner,(0)}}(X,s)=A(m)\mathrm{i}^{m}\mathrm{erfc}(X)+B(m)\mathrm{i}^{m}\mathrm{erfc}(-X),$
(91)

for $m=s+\zeta(X=0)$ is a natural number, and the repeated integral of the
complementary error function is defined in Eq. (39). This restriction on $m$
to natural numbers is justified, as only the solutions for integer $m$
contribute when we return from the Laplace domain to the time domain.

In order to obtain expressions for the functions $A(m)$ and $B(m)$, we apply a
matching procedure, i.e., demand that at some suitably chosen intermediate
point the inner and outer solution must be equal to each other. This we do for
both outer solutions. The matching conditions are given in the main text in
Eq. (40) and Eq. (41). The matching condition Eq. (41) straightforwardly
yields $B(m)=0$. We can obtain $A(m)$ by using the asymptotic limit for the
repeated integral of the complementary error function

$\lim\limits_{X\to-\infty}\mathrm{i}^{m}\mathrm{erfc}(X)\sim\frac{2(-X)^{m}}{\Gamma(m+1)},$
(92)

and the asymptotic limit of the outer solution for $x<1$ as

$\lim\limits_{x\to
1^{-}}\hat{\nu}_{0}(x<1,s)\sim(1-x)^{m}\exp\\{s\tau_{\mathrm{stat}}+\omega_{*}\tau_{\mathrm{non-
stat}}\\},$ (93)

with $\tau_{\mathrm{stat}}$ given by Eq. (43) and $\tau_{\mathrm{non-stat}}$
by Eq. (44). Inserting these expressions in Eq. (41) yields

$A(m)=s^{-1}\varepsilon^{m}\Gamma(m+1)\exp\\{s\tau_{\mathrm{stat}}+\omega_{*}\tau_{\mathrm{non-
stat}}\\},$ (94)

resulting in the inner solution

$\hat{\nu}_{\mathrm{inner}}(x,s)=\frac{1}{2s}\varepsilon^{m}\Gamma(m+1)\exp\\{s\tau_{\mathrm{stat}}+\omega_{*}\tau_{\mathrm{non-
stat}}\\}\mathrm{i}^{m}\mathrm{erfc}\left(\frac{x-1}{\varepsilon}\right),$
(95)

which is only valid near $x=1$.

Next, we focus on the inverse Laplace transform of the nucleation flux derived
in the main text in Eq. (55). For ease of notation, we write it as

$\widehat{J}(X,s)=G(X,s)\widehat{J}_{\mathrm{stat}}(X,s).$ (96)

with

$G(X,s)=\left[\frac{s}{(s-\omega_{*})}\exp\\{-\omega_{*}\left(\tau_{\mathrm{stat}}-\tau_{\mathrm{non-
stat}}\right)\\}\right].$ (97)

To execute the inverse Laplace transform, we make use of the convolution
theorem, which allows us two write the Laplace inverted Eq. (96) as

$J(X,\tau)=\int_{0}^{\tau}g(X,\tau^{\prime})J_{\mathrm{stat}}(X,\tau-\tau^{\prime})\mathrm{d}{\tau}^{\prime},$
(98)

where $g(X,t)$ reads

$g(X,\tau)=\left(\delta(\tau)+\omega_{*}\mathrm{e}^{\omega_{*}\tau}\right)\exp\\{-\omega_{*}\left(\tau_{\mathrm{stat}}-\tau_{\mathrm{non-
stat}}\right)\\},$ (99)

with $\delta(\tau)$ the delta function, and
$J_{\mathrm{stat}}(X,\tau-\tau^{\prime})$ is the stationary nucleation flux
defined in Eq. (58).

We first rewrite Eq. (96) as

$J(X,\tau)\exp\\{\omega_{*}\left(\tau_{\mathrm{stat}}-\tau_{\mathrm{non-
stat}}\right)\\}=J_{\mathrm{stat}}(X,\tau)+\int_{0}^{\tau}\omega_{*}\mathrm{e}^{\omega_{*}\tau}J_{\mathrm{stat}}(X,\tau^{\prime})\mathrm{d}{\tau}^{\prime}.$
(100)

Next, we use integration by parts to obtain

$\displaystyle
J(X,\tau)\exp\\{\omega_{*}\left(\tau_{\mathrm{stat}}-\tau_{\mathrm{non-
stat}}\right)\\}=$ (101) $\displaystyle\quad
J_{\mathrm{stat}}(X,\tau)+\left[\mathrm{e}^{\omega_{*}\tau^{\prime}}J_{\mathrm{stat}}(X,\tau-\tau^{\prime})\right]_{\tau^{\prime}=0}^{\tau^{\prime}=\tau}-\int_{0}^{\tau}\mathrm{e}^{\omega_{*}\tau^{\prime}}\frac{\mathrm{\partial}{}}{\mathrm{\partial}{\tau^{\prime}}}J_{\mathrm{stat}}(X,\tau-\tau^{\prime})\mathrm{d}{\tau}^{\prime}.$

This reduces to

$\displaystyle
J(X,\tau)\exp\\{\omega_{*}\left(\tau_{\mathrm{stat}}-\tau_{\mathrm{non-
stat}}\right)\\}=\left[\mathrm{e}^{\omega_{*}\tau}J_{\mathrm{stat}}(X,0)\right]-\int_{0}^{\tau}\mathrm{e}^{\omega_{*}\tau^{\prime}}\frac{\mathrm{\partial}{}}{\mathrm{\partial}{\tau^{\prime}}}J_{\mathrm{stat}}(X,\tau-\tau^{\prime})\mathrm{d}{\tau}^{\prime}.$
(102)

We set the first term equal to zero, as the initial stationary flux just after
an instantaneous quench into a metastable state should be zero. In reality, we
find that $J_{\mathrm{stat}}(X,0)$ contains a (minor) inconsistency as it is
not actually equal to zero. It is sufficiently small that we are nevertheless
justified in setting it to zero. See the discussion near Eq. (58). The
remaining integral reads

$J(X,\tau)=\exp\\{-\omega_{*}\left(\tau_{\mathrm{stat}}-\tau_{\mathrm{non-
stat}}\right)\\}\int_{0}^{\tau}\mathrm{e}^{\omega_{*}\tau^{\prime}}\frac{\mathrm{\partial}{}}{\mathrm{\partial}{\tau^{\prime}}}J_{\mathrm{stat}}(X,\tau-\tau^{\prime})\mathrm{d}{\tau}^{\prime}.$
(103)

We now use that $J_{\mathrm{stat}}(X,\tau)=J_{\mathrm{SS}}f(X,\tau)$, with
$J_{\mathrm{SS}}$ the steady-state nucleation flux and $f(X,\tau)$ a function
that accounts for the initial induction period, see Eq. (58). Next, we define
$J_{\mathrm{QSS}}(0)=J_{\mathrm{SS}}$ to correspond the the (quasi-)steady-
state flux using the initial conditions. This allows us to write

$J(X,\tau)=\int_{0}^{\tau}{J_{\mathrm{QSS}}(0)}\mathrm{e}^{\omega_{*}\tau^{\prime}}\exp\\{-\omega_{*}\left(\tau_{\mathrm{stat}}-\tau_{\mathrm{non-
stat}}\right)\\}\frac{\mathrm{\partial}{}}{\mathrm{\partial}{\tau^{\prime}}}\frac{J_{\mathrm{stat}}(X,\tau-\tau^{\prime})}{J_{\mathrm{QSS}}(0)}\mathrm{d}{\tau}^{\prime}.$
(104)

Using Eq. (57) and presuming that the time drift of the pre-exponential in the
quasi-steady-state nucleation flux is small relative to the contributions due
to the exponential, this can be rewritten as

$J(X,\tau)=\int_{0}^{\tau}J_{\mathrm{QSS}}\left(\tau^{\prime}-\tau_{\mathrm{stat.}}+\tau_{\mathrm{non-
stat}}\right)\frac{J_{\mathrm{stat}}(X,\tau-\tau^{\prime})}{J_{\mathrm{QSS}}(0)}\mathrm{d}{\tau}^{\prime}.$
(105)

After re-introducing the non-stretched variable $x$ and using the definition
$\tau_{\mathrm{cor}}=\tau_{\mathrm{non-stat}}-\tau_{\mathrm{stat.}}$, we find
Eq. (56).

##  Appendix B 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 Turnbull and Vonnegut (1952); Kelton and Greer (2010);
Djikaev  _et al._ (2003). Hence, we use it for illustrative purposes.

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

Focusing first on the thermodynamics of heterogeneous nucleation, where we
neglect any contribution from a line tension, the nucleation free energy is
given by Kashchiev (2000)

$\beta\Delta W_{\mathrm{Het}}=-n\beta\Delta\mu+f(\phi)\sigma n^{2/3},$ (106)

where $\sigma$ is the dimensionless surface tension between the solution and
the nucleus (see also Eq. (4) in the main text), and the scaling function
$f(\phi)$ 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 type of interfaces have different surface tensions
Kashchiev (2000); Kelton and Greer (2010). In general, additional physics,
such as that arising from the presence of a line tension, can be incorporated
in the scaling function $f(\phi)$, in which case the scaling function might
also depend on the cluster size $n$. Eq. (106) shows that the difference
between homogeneous and heterogeneous nucleation can, in general, be accounted
for by a renormalized surface tension. This can lower Turnbull and Vonnegut
(1952); Kelton and Greer (2010) or increase Djikaev  _et al._ (2003) the
nucleation barrier with respect to homogeneous nucleation. We refer to
Refs.Turnbull and Vonnegut (1952); Kashchiev (2000); Kelton and Greer (2010)
for details on the scaling function $f(\phi)$.

Note that we do not need to subtract the self-consistent correction from the
free energy to set the free energy to form a monomer to zero as we did in the
case of homogeneous nucleation, see Section. II in the main text. For
heterogeneous nucleation, we must distinguish between a monomer in free
solution and a monomeric cluster, that is, a monomer that is adsorbed onto a
surface, adding an adsorption free energy to the free energy barrier. The
change in chemical potential $\Delta\mu$ we presume similar to that in
homogeneous nucleation, that is, we presume the chemical potential to be given
by Eq. (5), where the density $\rho(t)$ acts as the local solute density at
substrate or solution-gas interface and $\rho_{\infty}$ 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(\phi)$ and if we
exclude the self-consistent monomer correction. These quantities are given by

$\beta\Delta
W^{*}_{\mathrm{Het}}=\frac{1}{3}n_{*}^{2/3}f(\phi)\sigma=3\varepsilon^{-2},$
(107) $n^{*}_{\mathrm{Het}}=\left(\frac{2f(\phi)\sigma}{3\ln
S(t)}\right)^{3},$ (108)

and

$Z_{\mathrm{Het}}=\sqrt{-\frac{1}{2\pi}\frac{\mathrm{\partial}^{2}{\Delta
W}}{\mathrm{\partial}{n}^{2}}}\Bigg{|}_{n=n_{*}}=\sqrt{\frac{f(\phi)\sigma}{\pi}}\frac{1}{3}{n_{*}}^{-2/3}.$
(109)

The expression for the equilibrium distribution does change somewhat, and is
now given by

$F_{\mathrm{eq}}(n,t)=\rho_{0}e^{-\beta\Delta W_{\mathrm{Het}}},$ (110)

where $\rho_{0}$ is the surface density of nucleation sites, which for a flat
interface is constant in time as evaporation does not affect this quantity.

The kinetics of heterogeneous nucleation are slightly different, since the
continuity equation does not need to be supplemented by a source term. The
reason is of course that the number of clusters per surface area does not
change in time due to solvent evaporation. Hence, we find

$\frac{\mathrm{\partial}{\rho(n,t)}}{\mathrm{\partial}{t}}=-\frac{\mathrm{\partial}{J(n,t)}}{\mathrm{\partial}{n}},$
(111)

where $J(n,t)$ is still given by the Zeldovich form Eq. (2). Although solvent
evaporation is not explicitly present in the continuity equation, it remains
present as a supersaturation in the nucleation barrier and hence a model for
solvent evaporation is still required.

Combining the thermodynamic and kinetic descriptions of heterogeneous
nucleation in the same way as we did in the main body of the article, we find

$\displaystyle\frac{\mathrm{\partial}{\nu(x,\tau)}}{\mathrm{\partial}{\tau}}$
$\displaystyle+\nu(x,\tau)\zeta(x,\tau)=$ (112)
$\displaystyle\quad\left[\frac{\varepsilon^{2}}{2}\frac{\mathrm{\partial}{\bar{k}^{+}(x)}}{\mathrm{\partial}{x}}-\varepsilon^{2}\bar{k}^{+}(x)\frac{\mathrm{\partial}{\beta\Delta
W_{\mathrm{Het}}(x,\tau)}}{\mathrm{\partial}{x}}+x\frac{\mathrm{\partial}{\ln
n_{*}(\tau)}}{\mathrm{\partial}{\tau}}\right]\frac{\mathrm{\partial}{\nu(x,\tau)}}{\mathrm{\partial}{x}}+$
$\displaystyle\quad\frac{\varepsilon^{2}}{2}\bar{k}^{+}(x)\frac{\mathrm{\partial}^{2}{\nu(x,\tau)}}{\mathrm{\partial}{x}^{2}},$

expressed in the co-moving $\\{x,\tau\\}$ co-ordinate system, see Sec. III,
and must be supplemented by the same boundary conditions as in the homogeneous
nucleation case, Eq. (29). Here, we define the function $\zeta(x,\tau)$ that
describes the non-stationarity

$\zeta(x,\tau)=-\frac{\mathrm{\partial}{}}{\mathrm{\partial}{\tau}}\beta\Delta
W_{\mathrm{Het}}(x,\tau)=\omega\psi(x),$ (113)

see Eq. (23) in the main text. For the non-stationarity coefficient $\omega$
we must account not only for the changing supersaturation, but also consider
that the function $f(\phi)$ 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. (106) 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 can 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

$\omega=n_{*}\frac{\mathrm{\partial}{\ln\chi(\tau)}}{\mathrm{\partial}{\tau}}-\sigma
n_{*}^{2/3}\frac{\mathrm{\partial}{f}}{\mathrm{\partial}{\phi}}\frac{\mathrm{\partial}{\phi}}{\mathrm{\partial}{\tau}},$
(114)

Here, the first term is due to the non-stationarity originating from the time-
dependent supersaturation and is also present in homogeneous nucleation, see
(23). The second term accounts for the change in the contact angles. Similar
contributions can be included if the scaling function $f(\phi)$ turns out to
depend on time.

Presuming that $\omega$ remains approximately constant within the co-moving
co-ordinate system described in the main text, the solution to Eq. (112) for
the nucleation flux is equivalent to the expression found for homogeneous
nucleation Eqs. (60) and Eq. (61). We must, however, replace the critical
nucleus size, the critical nucleation free energy, the Zeldovich factor and
the non-stationarity factor $\omega$ to by their respective expressions
presented in this Appendix. Moreover, the nucleation attachment mode might be
different, so another description for the attachment frequencies might be
required. The same holds for the induction time Eq. (72).

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

## References

## References

  * Richter  _et al._ (2017) L. J. Richter, D. M. DeLongchamp, and A. Amassian, Chemical Reviews 117, 6332 (2017).
  * Mei  _et al._ (2013) J. Mei, Y. Diao, A. L. Appleton, L. Fang, and Z. Bao, Journal of the American Chemical Society 135, 6724 (2013).
  * Ling  _et al._ (2018) H. Ling, S. Liu, Z. Zheng, and F. Yan, Small Methods 2, 1800070 (2018).
  * Di Carlo Rasi and Janssen (2019) D. Di Carlo Rasi and R. A. J. Janssen, Advanced Materials 31, 1806499 (2019).
  * Gaspar  _et al._ (2018) H. Gaspar, F. Figueira, L. Pereira, A. Mendes, J. Viana, and G. Bernardo, Materials 11, 2560 (2018).
  * McDowell  _et al._ (2018) C. McDowell, M. Abdelsamie, M. F. Toney, and G. C. Bazan, Advanced Materials 30, 1707114 (2018).
  * Zhang  _et al._ (2020) K. Zhang, Z. Wang, T. Marszalek, M. Borkowski, G. Fytas, P. W. M. Blom, and W. Pisula, Materials Horizons 7, 1631 (2020).
  * Michels  _et al._ (2021) J. J. Michels, K. Zhang, P. Wucher, P. M. Beaujuge, W. Pisula, and T. Marszalek, Nature Materials 20, 68 (2021).
  * Negi  _et al._ (2016) V. Negi, A. Lyulin, and P. Bobbert, Macromolecular Theory and Simulations 25, 550 (2016).
  * van Franeker  _et al._ (2015a) J. J. van Franeker, D. Westhoff, M. Turbiez, M. M. Wienk, V. Schmidt, and R. A. J. Janssen, Advanced Functional Materials 25, 855 (2015a).
  * van Franeker  _et al._ (2015b) J. J. van Franeker, M. Turbiez, W. Li, M. M. Wienk, and R. A. Janssen, Nature Communications 6, 6229 (2015b).
  * Lee  _et al._ (2020) T. Lee, A. V. Sanzogni, P. L. Burn, and A. E. Mark, ACS Applied Materials & Interfaces 12, 40548 (2020).
  * Schaefer  _et al._ (2016) C. Schaefer, J. J. Michels, and P. van der Schoot, Macromolecules 49, 6858 (2016).
  * Schaefer  _et al._ (2019) C. Schaefer, S. Paquay, and T. C. B. McLeish, Soft Matter 15, 8450 (2019).
  * Negi  _et al._ (2018) V. Negi, O. Wodo, J. J. van Franeker, R. A. J. Janssen, and P. A. Bobbert, ACS Applied Energy Materials 1, 725 (2018).
  * van Franeker  _et al._ (2015c) J. J. van Franeker, G. H. L. Heintges, C. Schaefer, G. Portale, W. Li, M. M. Wienk, P. van der Schoot, and R. A. J. Janssen, Journal of the American Chemical Society 137, 11783 (2015c).
  * Zhu  _et al._ (2020) L. Zhu, M. Zhang, G. Zhou, T. Hao, J. Xu, J. Wang, C. Qiu, N. Prine, J. Ali, W. Feng, X. Gu, Z. Ma, Z. Tang, H. Zhu, L. Ying, Y. Zhang, and F. Liu, Advanced Energy Materials 10, 1904234 (2020).
  * Virkar  _et al._ (2010) A. A. Virkar, S. Mannsfeld, Z. Bao, and N. Stingelin, Advanced Materials 22, 3857 (2010).
  * Wu  _et al._ (2020) Y. Wu, S. Schneider, C. Walter, A. H. Chowdhury, B. Bahrami, H.-C. Wu, Q. Qiao, M. F. Toney, and Z. Bao, Journal of the American Chemical Society 142, 392 (2020).
  * Tennyson  _et al._ (2019) E. M. Tennyson, T. A. S. Doherty, and S. D. Stranks, Nature Reviews Materials 4, 573 (2019).
  * Kashchiev (2000) D. Kashchiev,  _Nucleation_, 1st ed. (Elsevier, Oxford, 2000).
  * Auer and Frenkel (2001) S. Auer and D. Frenkel, Nature 409, 1020 (2001).
  * Filion  _et al._ (2010) L. Filion, M. Hermes, R. Ni, and M. Dijkstra, The Journal of Chemical Physics 133, 244115 (2010).
  * Filion  _et al._ (2011) L. Filion, R. Ni, D. Frenkel, and M. Dijkstra, The Journal of Chemical Physics 134, 134901 (2011).
  * Reiss  _et al._ (1997) H. Reiss, W. K. Kegel, and J. L. Katz, Physical Review Letters 78, 4506 (1997).
  * Reiss  _et al._ (1998) H. Reiss, W. K. Kegel, and J. L. Katz, The Journal of Physical Chemistry A 102, 8548 (1998).
  * Drossinos and Kevrekidis (2003) Y. Drossinos and P. G. Kevrekidis, Physical Review E 67, 026127 (2003).
  * Aasen  _et al._ (2020) A. Aasen, D. Reguera, and O. Wilhelmsen, Physical Review Letters 124, 045701 (2020).
  * Aasen  _et al._ (2023) A. Aasen, O. Wilhelmsen, M. Hammer, and D. Reguera, The Journal of Chemical Physics 158, 114108 (2023).
  * Pan and Chandler (2004) A. C. Pan and D. Chandler, The Journal of Physical Chemistry B 108, 19681 (2004).
  * Peters and Trout (2006) B. Peters and B. L. Trout, The Journal of Chemical Physics 125, 054108 (2006).
  * Lechner  _et al._ (2011) W. Lechner, C. Dellago, and P. G. Bolhuis, The Journal of Chemical Physics 135, 154110 (2011).
  * Oxtoby (1998) D. W. Oxtoby, Accounts of Chemical Research 31, 91 (1998).
  * Oxtoby and Evans (1988) D. W. Oxtoby and R. Evans, The Journal of Chemical Physics 89, 7521 (1988).
  * Reguera  _et al._ (2005) D. Reguera, J. M. Rubí, and J. M. G. Vilar, The Journal of Physical Chemistry B 109, 21502 (2005).
  * Durán-Olivencia  _et al._ (2018) M. A. Durán-Olivencia, P. Yatsyshin, S. Kalliadasis, and J. F. Lutsko, New Journal of Physics 20, 083019 (2018).
  * Lutsko and Durán-Olivencia (2013) J. F. Lutsko and M. A. Durán-Olivencia, The Journal of Chemical Physics 138, 244908 (2013).
  * Lutsko (2018) J. F. Lutsko, New Journal of Physics 20, 103015 (2018).
  * Kuhnhold  _et al._ (2019) A. Kuhnhold, H. Meyer, G. Amati, P. Pelagejcev, and T. Schilling, Physical Review E 100, 052140 (2019).
  * Vehkamaki (2006) H. Vehkamaki,  _Classical Nucleation Theory in Multicomponent Systems_ (Springer-Verlag, Berlin/Heidelberg, 2006).
  * Wöhler and Schilling (2022) W. Wöhler and T. Schilling, Physical Review Letters 128, 238001 (2022).
  * Shneidman and Hänggi (1994) V. A. Shneidman and P. Hänggi, Physical Review E 49, 641 (1994).
  * Shneidman (2007) V. A. Shneidman, The Journal of Chemical Physics 127, 041102 (2007).
  * Trinkaus and Yoo (1987) H. Trinkaus and M. H. Yoo, Philosophical Magazine A 55, 269 (1987).
  * Shneidman (2011) V. A. Shneidman, Physical Review E 84, 031602 (2011).
  * Shneidman (2010) V. A. Shneidman, Physical Review E 82, 031603 (2010).
  * Bornside  _et al._ (1989) D. E. Bornside, C. W. Macosko, and L. E. Scriven, Journal of Applied Physics 66, 5185 (1989).
  * Schaefer  _et al._ (2017) C. Schaefer, J. J. Michels, and P. van der Schoot, Macromolecules 50, 5914 (2017).
  * Reguera and Rubi (2003a) D. Reguera and J. M. Rubi, The Journal of Chemical Physics 119, 9877 (2003a).
  * Reguera and Rubi (2003b) D. Reguera and J. M. Rubi, The Journal of Chemical Physics 119, 9888 (2003b).
  * Kalikmanov (2013) V. Kalikmanov,  _Nucleation Theory_, Lecture Notes in Physics, Vol. 860 (Springer Netherlands, Dordrecht, 2013).
  * Ruckenstein and Nowakowski (1990) E. Ruckenstein and B. Nowakowski, Journal of Colloid and Interface Science 137, 583 (1990).
  * Wu (1996) D. T. Wu, in  _Solid State Physics_, Vol. 50, edited by H. Ehrenreich and F. Spaepen (Academic Press, San Diego, 1996) pp. 37–197.
  * Slezov and Schmelzer (1998) V. V. Slezov and J. W. Schmelzer, Journal of Physics and Chemistry of Solids 59, 1507 (1998).
  * Schenter  _et al._ (1999) G. K. Schenter, S. M. Kathmann, and B. C. Garrett, Physical Review Letters 82, 3484 (1999).
  * Note (1) The validity and accuracy of this equation, specifically related to treating $n$ as a continuous variable, are discussed by Wu (1992b).
  * Zeldovich (1943) Y. B. Zeldovich, Acta Physicochimica U.R.S.S. 18, 1 (1943).
  * Danglad-Flores  _et al._ (2018) J. Danglad-Flores, S. Eickelmann, and H. Riegler, Chemical Engineering Science 179, 257 (2018).
  * Gumpert  _et al._ (2023) F. Gumpert, A. Janßen, C. J. Brabec, H.-J. Egelhaaf, J. Lohbreier, and A. Distler, Engineering Applications of Computational Fluid Mechanics 17, 2242455 (2023).
  * Kožíšek  _et al._ (2011) Z. Kožíšek, K. Sato, S. Ueno, and P. Demo, The Journal of Chemical Physics 134, 094508 (2011).
  * Tiwari and van der Schoot (2016) N. S. Tiwari and P. van der Schoot, The Journal of Chemical Physics 144, 235101 (2016).
  * Dhont (2003) J. K. Dhont,  _An Introduction to Dynamics of Colloids_ (Elsevier, Amsterdam, 2003).
  * Crank (1975) J. Crank,  _The Mathematics of Diffusion_ (Clarendon Press, Oxford, 1975).
  * Risken (1996) H. Risken,  _Fokker-Planck Equation_ (Springer, Berlin, 1996) pp. 63–95.
  * Note (2) Due to the second-order reaction kinetics, we expect the attachment frequencies to be given by $k^{+}(x,\tau)=\rho(1,\tau)\kappa(x,\tau)$ with $\kappa(x,\tau)$ a reaction coefficient associated with monomer attachment to a cluster of size $x$. As a result, the reduced attachment frequencies are $\bar{k}^{+}(x,\tau)=\kappa(x,\tau)/\kappa(1,\tau)$. If $n_{*}$ and $\varepsilon$ can be treated as essentially independent of time, then, for all relevant models for the aggregation kinetics known to us, this holds for $\kappa(x,\tau)/\kappa(1,\tau)$ as well Turnbull and Fisher (1949); Kashchiev (2000); Kalikmanov (2013).
  * Cahn (1995) J. W. Cahn, MRS Proceedings 398, 425 (1995).
  * Nayfeh (2011) A. H. Nayfeh,  _Introduction to Perturbation Techniques_ (John Wiley & Sons, 2011).
  * Shi  _et al._ (1990) G. Shi, J. H. Seinfeld, and K. Okuyama, Physical Review A 41, 2101 (1990).
  * Hoyt and Sundar (1993) J. Hoyt and G. Sundar, Scripta Metallurgica et Materialia 29, 1535 (1993).
  * Turnbull and Fisher (1949) D. Turnbull and J. C. Fisher, The Journal of Chemical Physics 17, 71 (1949).
  * Shneidman (1987) V. Shneidman, Soviet Physics. Technical Physics 37, 76 (1987).
  * Shneidman (1988) V. Shneidman, Soviet Physics. Technical Physics 33, 1338 (1988).
  * Olver  _et al._ (2010) F. W. Olver, D. W. Lozier, R. F. Boisverrt, and C. W. Clark, eds.,  _NIST Handbook of Mathematical Functions_ (Cambridge University Press, New York, 2010).
  * Wu (1992a) D. T. Wu, The Journal of Chemical Physics 97, 2644 (1992a).
  * Shneidman and Weinberg (1992) V. A. Shneidman and M. C. Weinberg, The Journal of Chemical Physics 97, 3621 (1992).
  * Aksenov  _et al._ (2005) V. Aksenov, T. Tropin, M. V. Avdeev, V. Priezzhev, and J. Schmelzer, Physics of Particles and Nuclei 36, S52 (2005).
  * Wedekind  _et al._ (2007) J. Wedekind, R. Strey, and D. Reguera, The Journal of Chemical Physics 126, 134103 (2007).
  * Spinella  _et al._ (1998) C. Spinella, S. Lombardo, and F. Priolo, Journal of Applied Physics 84, 5383 (1998).
  * Rasmussen (1982) D. H. Rasmussen, Journal of Crystal Growth 56, 45 (1982).
  * Wedekind  _et al._ (2009) J. Wedekind, G. Chkonia, J. Wölk, R. Strey, and D. Reguera, The Journal of Chemical Physics 131, 114506 (2009).
  * Binder (1984) K. Binder, Physical Review A 29, 341 (1984).
  * Bartell and Wu (2007) L. S. Bartell and D. T. Wu, The Journal of Chemical Physics 127, 174507 (2007).
  * Sosso  _et al._ (2016) G. C. Sosso, J. Chen, S. J. Cox, M. Fitzner, P. Pedevilla, A. Zen, and A. Michaelides, Chemical Reviews 116, 7078 (2016).
  * de Hoog and Lekkerkerker (1999) E. H. A. de Hoog and H. N. W. Lekkerkerker, The Journal of Physical Chemistry B 103, 5274 (1999).
  * van der Schoot (1999) P. van der Schoot, The Journal of Physical Chemistry B 103, 8804 (1999).
  * Asherie  _et al._ (1996) N. Asherie, A. Lomakin, and G. B. Benedek, Physical Review Letters 77, 4832 (1996).
  * Hagen  _et al._ (1993) M. H. J. Hagen, E. J. Meijer, G. C. A. M. Mooij, D. Frenkel, and H. N. W. Lekkerkerker, Nature 365, 425 (1993).
  * Tumbek and Winkler (2012) L. Tumbek and A. Winkler, Surface Science 606, L55 (2012).
  * Yu  _et al._ (2019) L. Yu, D. Qian, S. Marina, F. A. A. Nugroho, A. Sharma, S. Hultmark, A. I. Hofmann, R. Kroon, J. Benduhn, D.-M. Smilgies, K. Vandewal, M. R. Andersson, C. Langhammer, J. Martín, F. Gao, and C. Müller, ACS Applied Materials & Interfaces 11, 21766 (2019).
  * Bangsund  _et al._ (2019) J. S. Bangsund, T. R. Fielitz, T. J. Steiner, K. Shi, J. R. Van Sambeek, C. P. Clark, and R. J. Holmes, Nature Materials 18, 725 (2019).
  * Kelton (1991) K. Kelton, Solid State Physics 45, 75 (1991).
  * Wu (1992b) D. T. Wu, The Journal of Chemical Physics 97, 1922 (1992b).
  * Hu  _et al._ (2020) H. Hu, M. Singh, X. Wan, J. Tang, C.-W. Chu, and G. Li, Journal of Materials Chemistry A 8, 1578 (2020).
  * Zeng  _et al._ (2020) L. Zeng, S. Chen, K. Forberich, C. J. Brabec, Y. Mai, and F. Guo, Energy & Environmental Science 13, 4666 (2020).
  * Habraken  _et al._ (2013) W. J. E. M. Habraken, J. Tao, L. J. Brylka, H. Friedrich, L. Bertinetti, A. S. Schenk, A. Verch, V. Dmitrovic, P. H. H. Bomans, P. M. Frederik, J. Laven, P. van der Schoot, B. Aichmayer, G. de With, J. J. DeYoreo, and N. A. J. M. Sommerdijk, Nature Communications 4, 1507 (2013).
  * Levitsky  _et al._ (2021) A. Levitsky, S. A. Schneider, E. Rabkin, M. F. Toney, and G. L. Frey, Materials Horizons 8, 1272 (2021).
  * Turnbull and Vonnegut (1952) D. Turnbull and B. Vonnegut, Industrial & Engineering Chemistry 44, 1292 (1952).
  * Kelton and Greer (2010) K. Kelton and A. Greer,  _Nucleation in condensed matter: Application in Materials and Biology_ , 1st ed. (Pergamon, Amsterdam, 2010).
  * Djikaev  _et al._ (2003) Y. S. Djikaev, A. Tabazadeh, and H. Reiss, The Journal of Chemical Physics 118, 6572 (2003).