Article pubs.acs.org/crystal
Statistical Analysis of Series of Detection Time Measurements for the Estimation of Nucleation Rates Giovanni Maria Maggioni,†,‡ Luca Bosetti,†,‡ Elena dos Santos,† and Marco Mazzotti*,† †
Separation Processes Laboratory, Department of Mechanical and Process Engineering, ETH Zurich, Zurich 8092, Switzerland
ABSTRACT: Primary nucleation is a stochastic process; hence, detection times are statistically distributed even though experiments are repeated at the same conditions. This contribution aims at defining and discussing a method to perform an accurate statistical analysis of detection times; it focuses on three main aspects. First, we develop an accurate experimental protocol and set criteria on temperature, T, and supersaturation, S, to accept measurements as part of the experimental series. Applying such protocol to the isonicotinamide−ethanol system at several different supersaturations, we perform multiple series of isothermal nucleation experiments, calculating for each series the associated empirical cumulative distribution function (eCDF) of the detection times. Second, exploiting a set of statistical tools, we investigate whether the repetitions of the experiments at the same supersaturation conditions belong to the same stochastic process, assessing also the possibility to combine them in one single distribution. Third, we estimate for the eCDFs the values of its associated nucleation rate and derive an analytical expression to calculate the propagation of the uncertainty from the eCDF to the nucleation rate estimate.
1. INTRODUCTION Due to the relevance of nucleation from solution for crystallization processes, the experimental and theoretical characterization of its rate is of crucial importance. According to classical nucleation theory,1,2 primary nucleation is an activated process and hence stochastic in nature.2,3 When measuring nucleation times and characterizing nucleation rates, such stochastic nature poses a number of challenges which have been the subject of many papers published in the past4,5 as well as recently.6−10 It is fair to say that characterizing nucleation rates accurately and reliably is still an unsolved scientific problem.11 Issues in measuring nucleation rates can be categorized as follows. First, because the nucleation rate is strongly sensitive to operating conditions,12 namely temperature and solute concentration, experimental operating conditions have to be tightly controlled to guarantee consistency from one experiment to the next. As a consequence, each experiment used to build a representative statistics has to be double-checked before acceptance, as thoroughly discussed recently in the literature.13,14 Second, as each and every nucleation experiment exhibits a different detection time, it is clear that a large number of experiments are needed to build representative statistics. It is therefore necessary to establish quantitative statistical criteria © XXXX American Chemical Society
that allow deciding when different series of experiments carried out at nominally identical conditions can be combined and considered together to estimate nucleation rates. We believe that the best approach involves the use of Kolmogorov−Smirnov statistical tests (KS-test), as proposed by Little, Sear, and Keddie.13 Third, care has to be taken in defining a protocol to extract nucleation rates from distributed nucleation times, and attention has to be given to how to determine the confidence limits of such estimates.4,5,15,16 A recent paper exploited the results of a numerical analysis to quantify the generic inherent uncertainty that has to be associated with experimental estimates of growth rates obtained in small volume experiments.14 Finally, as discussed in some detail in recent experimental papers17,18 and in our previous two contributions,12,19 nucleation is detected after crystals have grown sufficiently. Therefore, nucleation rate and growth rate are necessarily correlated. The estimate of one depends on that of the other, unless one can assume that crystals grow from nucleation to detection in a negligible time as compared to the nucleation time.13,14,17,18 Received: July 21, 2017 Revised: September 13, 2017
A
DOI: 10.1021/acs.cgd.7b01014 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
Figure 1. Experimental setup and protocol. On the left-hand side: schematic representation of one crystallizer (∼1.8 mL) in the C16-TMS with the thermocouple inserted through the cap and the laser sensor. On the right-hand side: typical experimental cycles with INA-EtOH performed at S = 1.72 with dissolution at 60 °C and crystallization at 25 °C. The detection event corresponds to the transmissivity drop and the tiny temperature peak. The saturation temperature corresponds to S = 1.72.
This work aims at addressing the first three issues above (to each one section will be devoted); as to the fourth issue, we assume that the growth time is negligible in all the experimental series presented. The novel aspects with respect to previous literature can be summarized as follows. First, we carried out nucleation experiments with isonicotinamide in ethanol in a Crystal16 setup (Technobis Crystallization Systems), where it was possible to crash-cool the solution in less than 2 min. This device was custom-modified to allow for online temperature monitoring, which provides a vital piece of information for determining whether experiments were conducted exactly at the same conditions. Second, we applied the KS-tests through the whole range of significance levels, and not only at 0.05 as typically done, 13 thus providing a much richer insight on the reproducibility of the experimental series and on the accuracy of the nucleation rate estimate. On the basis of the KS-test, we also developed a statistical method that allows us to decide if one can combine an arbitrary number of repeated series, not just two of them, potentially obtaining much larger pools of data. Third, instead of using a numerical analysis to establish confidence intervals on the estimated nucleation rates,14 we derived and carried out an exact sensitivity analysis for a Poisson model of nucleation, independent of the specific substance used and its nucleation rate parameters.
steadily below such values so as to avoid experimental errors due to signal noise. The modified version of the instrument, named Crystal16 Temperature Measurement System, was equipped with an independent thermocouple in each reactor. The system was designed to avoid interference of the thermocouple with the turbidity measurement and the magnetic bar stirrer (Figure 1). The thermocouples are k-type probes (Iconel 600) 45 mm long and 1.5 mm in diameter, purchased from and certified by Picolog. Because the temperature inside each vial is measured, the supersaturation attained in each experiment is estimated more accurately. It is worth noting that the thermocouples provide evidence that crystallization has occurred in the system, independent of transmissivity, as shown in Figure 1. The drop in transmissivity occurs at the same time when the temperature changes. Due to the isothermal conditions imposed by the thermostat, such changes must be connected to a release of heat due to crystallization, which occurs so suddenly that is not immediately compensated by the thermostat. 2.2. Experimental Protocol. 2.2.1. Experimental Procedure. For each series of experiments, a stock solution (60 g) was prepared by dissolving the desired amount of solute in the solvent at a temperature above saturation (60 °C) in a closed vessel, and used immediately thereafter. For each series of experiments, a new set of 16 glass vials was used as crystallizers and, before each new experimental series, all caps and stirrers were first washed with water and ethanol to dissolve any solids from previous experiments and then dried with compressed air to remove mechanically any possible residual fine particle. Each position in the Crystal16 was assigned to a specific stirrer: given that only the glass vials changed from one series to another, this protocol allowed us to track possible biases due to a specific stirrer, or a specific position, when more repetitions at the same conditions (i.e., same supersaturation) were performed. The caps were not kept in the same position: after accurate tests, their assignment to a specific position was proved neither to influence the results nor improve the procedure. Before filling, the empty vials were weighed, each together with its stirrer and cap. The following precautions were taken during the loading of the vials with clear, undersaturated solution, typically at a temperature above room temperature: (1) the vials were filled as fast as possible to avoid cooling below the saturation temperature, with the stock solution at 60 °C; (2) the syringe used to fill the vials were clean at any moment; in the case crystalline formations are observed on its outer or inner part, and the syringe was substituted with a new one; and (3) after loading, the vials were immediately closed to prevent solvent evaporation. The filling protocol was designed to have conditions as homogeneous as possible inside each crystallizer. The filled and closed vials were then weighed again to determine the actual mass loaded in each position. The coefficient of variation of the loaded masses was less than 2% in all experiments. The vials were placed in their positions in the Crystal16,
2. MATERIALS AND METHODS In this section, we describe first the equipment, the experimental protocol, and the material used in the experiments, and then we summarize the background theory with focus on statistical analysis of data and estimation of nucleation rates. 2.1. Equipment. All nucleation experiments were performed in a customized version of the Crystal16, engineered and built by Technobis Crystallization Systems according to our specifications. The instrument consists of 16 independent reactors, each of them a cylindrical glass vial of about 1.8 mL capacity. The presence of solids in solution was determined by transmissivity measurements. Transmissivity is a measure of the light, which travels unscattered through the solution from a source to a sensor. In the Crystal16, transmissivity is monitored independently in each reactor; the light source is a red laser. The reference state, i.e. 100% transmissivity, is taken as the unsaturated, clear solution at the dissolution temperature. Particles in suspension are detected when transmissivity steadily drops below 90% and remains B
DOI: 10.1021/acs.cgd.7b01014 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
and accurately cleaned thermocouples were immersed into the solutions. A typical cycle consisted of two main steps. First, the solution was kept at the dissolution temperature, Tdis, 60 °C in this work, for a certain amount of time (60 min) to ensure complete crystal dissolution. The dissolution temperature was selected to be sufficiently higher than the saturation temperature but also much lower than the boiling point of the solvent. Then, the solution was cooled (at a rate of 20 °C/min in this work) to the crystallization temperature, Tcry, i.e. 25 °C in all experiments reported in this paper (see Section 2.2.3 for further details about temperature accuracy). The duration of the isothermal step of each cycle varies according to the supersaturation (30 to 240 min from the highest to the lowest value of S in this work). Finally, after the crystallization step, the system was heated at Tdis to dissolve the particles before starting a new cycle. For the entire duration of the experiments, the vials were stirred at 700 rpm with a standard magnetic bar. Once the whole series of cycles was completed, the vials were weighed for a third time to quantify the amount of solvent evaporated during the experiment. Note that, while for each cycle only up to 16 nucleation events and the corresponding detection times can be measured, the number of cycles in a complete series is limited by the extent of evaporation (see below, Section 2.2.3). Therefore, to obtain a large set of data points, typically more than one series at the same nominal conditions were performed. 2.2.2. Filtration of the Solution. The clear, dissolved parent solution was filtered at 60 °C before being loaded in the reactors. Filtration was performed during the loading of the vials: once the syringe was full with parent solution, a 0.22 μm PTFE filter was attached to the tip, and then the vial was filled with the solution passing through the filter. The filter was then removed from the tip of the syringe, which was accurately washed with the solvent and inserted again in the reactor for refilling. All the other steps described for the loading of vials were kept unchanged, but additional care was paid when filters had to be substituted. After loading four vials, the filter was usually replaced with a new one. The change could be necessary before, if crystals were observed in the filter, or if the pressure needed to empty the syringe increased. 2.2.3. Validation of Experimental Measurements. The actual crystallization temperature does not always match exactly the set temperature and that evaporation could occur during the cycles. Temperature changes and solvent evaporation affect the kinetics by modifying the actual supersaturation inside the vials. Thus, we defined criteria to validate each experimental point by verifying that it is indeed measured at the nominal operating conditions or very close to them.13,20 Such validation procedure guarantees that measurements are homogeneous and can be considered as sampled from the same statistical distribution. Note that three numbers characterize each series of experiments: the total number of experiments, Ntot, given simply by the number of cycles times 16; the number of experiments complying with the criteria for reproducibility (as discussed below), N; and the number of experiments where crystal detection was possible before the cooling period is completed and the vials are heated again (see Figure 1), n. Temperature Criterion. Because the thermocouples measure the temperature inside each vial for the entire duration of the experiment, i.e. between establishment of the nominal supersaturation and detection, it is possible to verify that during this whole time interval the temperature remains within a 0.5 °C range around the set value, i.e. it is of 25 ± 0.5 °C. Only experiments fulfilling this criterion yield a valid detection time measurement. Concentration Criterion. In this work, we assume a moderate change of the solute’s activity coefficient in a range between the actual solute concentration, c, and its temperature-dependent solubility, c*(T), i.e. we assume γ(T,c) ≃ γ*(T,c*) . Therefore, the supersaturation is defined as c S= c*(T ) (1)
concentration (hence the supersaturation) in the individual vials. Only the detection times whose corresponding supersaturation is less than 0.02 larger than the set value were considered valid measurements. 2.3. Materials. Isonicotinamide (INA), purity ≥99%, was purchased from Sigma-Aldrich, while ethanol (EtOH) from two suppliers, Fluka (Ethanol Laboratory Reagent, absolute, ≥99.5%) and EMD Millipore (Ethanol Absolute For Analysis Emsure Acs, Iso, Reag. Ph. Eur.). For the experiments where the solution was filtered, 13 mm syringe filters (PTFE hydrophobic, nonsterile, pore size 0.22 μm) were used. It is known that INA has different polymorphs.21,22 The supersaturation in this work is calculated for the form stable at 25 °C, named Form 1 by Hansen and co-workers.21 We measured the solubility of INA in ethanol at 25 °C gravimetrically and found it to be 0.093 g/g (gram of substance per gram of solvent); this is the average value of 6 independent measurements. 2.4. Empirical Cumulative Distribution Functions. Each series of Ntot experiments yields N acceptable experiments, fulfilling the criteria discussed above, and n values of the detection time, τ, where n is the number of experiments in which crystals were indeed detected (hence, N − n is the number of acceptable experiments in which no crystal was detected before the end of the experiment). Because nucleation is an activated phenomenon,1 the nucleation process is stochastic; hence, the measured detection time τ (which is the sum of the nucleation time, tN, and the growth time, tG, the time between nucleation and detection) is itself a stochastic variable distributed according to an intrinsic underlying cumulative probability function, F(τ). The n measured detection times form the ordered vector τ: τ = [τ1, τ2 , ..., τn]T : 0 ≤ τ1 ≤ τ2 ≤ ... ≤ τn < +∞
(2)
through which the following empirical cumulative distribution function (eCDF) F*N (τ) is defined as FN*(τ ) =
1 N
n
∑ 1[τj ,+∞)(τ)
(3)
j=1
where 1[a , +∞) is the indicator function on the interval [a, +∞), i.e. equal to 1 if a ≤ τ < ∞, and 0 otherwise. The vector τ constitutes a discrete, finite sample of the stochastic variable τ, whereas FN*(τ) is an estimate of F(τ). To illustrate the quantities and concepts above, we assume that nucleation is a Poisson process and that the nucleation time, i.e. tN = τ − tG, is distributed according to the following Poisson cumulative probability function of τ:
F(τ ) = 1 − exp(− JV (τ − tG))
(4)
We plotted in Figure 2 the function F(τ) (blue line) and the two eCDFs **(τ), computed according to eq 3 from two randomly FN*(τ) and FM sampled sets of detection times, with N = 20 (red) and M = 80 (yellow), respectively (the two sets were generated according to the inverse transform sampling algorithm using the uniform random number generator implemented in Matlab). It is apparent that, according to both intuition and theory,23 the quality of the estimation of F(τ) improves with increasing number of experimental measurements. 2.5. Estimation of Nucleation Rates. The ultimate goal of detection time experiments is to measure the nucleation rate, J. In this study, J was estimated for each supersaturation value S according to the method presented by Jiang et al.20 and Xiao et al.,14 which is based on two hypotheses. First, nucleation is a Poisson process, as to the cumulative distribution function defined by eq 4. Second, the growth time tG depends on the experimental conditions (supersaturation and temperature) and can be estimated for each series of experiments as the shortest detection time measured at such conditions. The value of J at each supersaturation for a given system volume V can be estimated by using the linearized form of eq 4 and by minimizing the following sum of the squares of the residues between experiments and model:
Note that because of the temperature criterion above, supersaturation has always been computed using the value of solubility at 25 °C, i.e. 0.094 g/g. Solvent evaporation was observed during the experiments and is thus to be considered. It was assumed that the same amount of ethanol evaporates during each cycle, progressively increasing the
n
Ψ=
∑ (log(1 − FN*(τj)) + JV (τj − tG))2 j=1
C
(5) DOI: 10.1021/acs.cgd.7b01014 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
the form of the corresponding eCDF (blue, red, yellow, and violet for series 1, 2, 3 and 4; the number of cycles and of retained measurements for each series is reported in the inset of the same figure). When confronted with this set of data or similar ones, the question is whether to consider them self-consistent or reproducible, i.e., whether to consider them as estimates of the same underlying cumulative probability function. If this were the case, then one could combine the four sets of data and obtain a single set of data with a much larger overall number of data points. The latter could be the basis for a more accurate estimate of the nucleation rate corresponding to the selected supersaturation level. In this section, we explore when and how multiple series of nucleation experiments can be combined by first providing the statistical basis for such investigation and then applying it to the four series shown in Figure 3. 3.2. Statistical Basis. The proper measure of the distance in the space of distribution functions between either an eCDF FN*(τ) and a cumulative distribution function F(τ), or a pair of independent eCDFs **(τ), is the infinite or Chebyshev norm defined as the FN*(τ) and FM maximum of the absolute values of the difference between the two distributions. Such a norm is in turn a random variable, whose distribution was studied first by Kolmogorov24 and Smirnov25 and then by many others. More specifically, we define the following two norms:
Figure 2. Plots of the cumulative probability function, F(τ), of eq 4 (with J = 10 000 #/(m3 s), tG = 0 s, and V = 1 mL), blue, and of two ** with M = 80, eCDFs sampled from it, i.e. FN* with N = 20, red, and FM yellow. Note that in this case n* = N and n** = M. In the case of the examples in Figure 2, while F(τ) was obtained with V = 1 mL and J = 10 000 [#/(m3 s)], the values of J estimated with 20 or 80 measurements are J = 9 000 [#/(m3 s)] and J = 12 400 [#/(m3 s)], respectively.
DN = FN*(τ ) − F(τ )
* *(τ ) − FN*(τ ) DM , N = FM
3. COMBINATION OF EXPERIMENTAL SERIES
(6)
∞
(7)
∞
We call both of them DQ, where Q = N in the former case, and Q = MN/ (M + N) in the latter case, as shown by Kolmogorov.24 We also define the null hypothesis that FN*(τ) is sampled from the continuous cumulative distribution function F(τ) in the first case, and that F*N (τ) **(τ) are sampled both independently from the same continuous and FM cumulative distribution function in the second case. Under the null hypothesis, DQ is a random variable, and the quantity √QDQ (called for brevity K-number in the first case and KS-number in the second) exhibits a limiting cumulative distribution function Φ(λ) given by
3.1. Repeated Experiments. In a deterministic process, reproducibility is the possibility to replicate the experimental results for the same set of conditions and operating parameters within a certain accuracy, which is controlled by the background noise affecting the measured signals. Thus, in an ideal, zero-noise deterministic system, the values of the variables measured in different repetitions of the same experiment do not change. On the contrary, the outcome of a stochastic process such as nucleation is by definition not reproducible in this sense. With reference to the experiments reported here, in principle, each series samples the same underlying distribution, but its measurements are expected to be different from the measurements obtained in a different series performed under the same operating conditions. This point is illustrated in Figure 2 using simulated data and, more importantly, in Figure 3, where 4 series of experiments carried out at the same temperature (25 °C) and supersaturation (S = 1.41) are shown in
∞
Φ(λ) = 1 − 2 ∑ (− 1)i − 1exp(− 2i 2λ 2) i=1
(8)
This should be interpreted as follows. Under the null hypothesis, the probabilities that the product √QDQ is smaller than λ (this is called confidence level and indicated as P) or larger than λ (the significance level, α) are Q →+∞
P = Pr { Q DQ ≤ λ} ⎯⎯⎯⎯⎯⎯⎯⎯→ Φ(λ)
(9)
Q →+∞
α = Pr { Q DQ > λ} ⎯⎯⎯⎯⎯⎯⎯⎯→ 1 − Φ(λ)
(10)
The probabilities on the l.h.s. of these equations are in general functions of Q and λ, but for Q ≥ 35 (a condition that is always fulfilled in the experiments reported here), they are indistinguishable from the asymptotic expression on the r.h.s.26 The asymptotic relationships P = Φ(λ) = 1 − α are illustrated in Figure 4. Note that the definitions of confidence and significance levels provided by eqs 9 and 10 are mathematically rigorous and should be referred to whenever interpreting data using P and α. For data yielding values of the K- or KS-number corresponding to point A in the figure, the null hypothesis is fulfilled with a confidence level as low as 0.2 and with a significance as high as 0.8. For data corresponding to point B in the figure, the null hypothesis is fulfilled at a higher confidence level but at a lower significance level. Point R in the figure corresponds to the value which is used as a standard in statistics, namely a confidence level of 0.95 at a significance level of 0.05, i.e. corresponding to λ = 1.358. It is also worth noting that the following inequality on α (DKWM inequality, from Dvoretzky, Kiefer, Wolfowitz,27 and Massart28) holds true for every Q and λ:
Figure 3. Nucleation experiments of INA in EtOH. Cumulative distribution functions of 4 experimental series at the same nominal conditions, namely S = 1.41. For series 3, the corresponding uncertainty area was drawn: its boundaries were computed with P̃ = 0.95 according to the procedure detailed in Section 3.3. Note that all series are mostly within the uncertainty area. D
DOI: 10.1021/acs.cgd.7b01014 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
corresponding KS-number λ = √QDN,M, and (iv) determine the corresponding P and α values through eqs 9 and 10 (or by reading the values in Figure 4). The KS-numbers for all 6 possible combinations of 2 series from the 4 in Figure 3 span the range 0.69 (pair 2−3) to 1.19 (pair 1−4) and are used to plot the crosses in Figure 5a on the background of the function
Figure 4. Kolmogorov distribution Φ(λ) (eq 8). Both the confidence level P (lower horizontal-axis) and the significance level α (upper horizontal-axis) are indicated. A and B represent two arbitrary points with P = 0.20 and 0.70, respectively, while R represents the standard value P̃ = 0.95, corresponding to α̃ = 0.05 and λ = 1.358. α = Pr { Q DQ > λ} ≤ ΦDKWM(λ) = 2exp(− 2λ 2)
(11)
Using ΦDKWM(λ) instead of Φ(λ) yields an error of less than 0.1% on α for α for λ ≥ 1 and α ≤ 0.271. The function ΦDKWM(λ) can be easily inverted, thus yielding:
λ=
1 ⎛⎜ 2 ⎞⎟ log 2 ⎝α⎠
(12)
which is a fully acceptable approximation when α ≤ 0.271. 3.3. Statistical Tests. We consider again the 4 experimental series illustrated in Figure 3, having a number of data points between 84 and 220, with the objective of merging them into a single series with a larger number of data points by combining as many of them as possible, i.e. by combining those that can be considered well reproduced under the chosen experimental conditions. Such discussion shall be based on two statistical tests, which are in turn based on the two norms defined in the previous section. The first test considers a single eCDF FN*(τ), the yellow series in Figure 3 having N = 220, and evaluates the probabilistic distance of the eCDF from the continuous cumulative distribution function F(τ) given by eq 6, under the null hypothesis that the former is indeed sampled from the latter.26 To this aim, one must (i) select a significance level α, e.g. α = 0.05, (ii) calculate the corresponding λ value either from eq 12 if α ≤ 0.271 or by inverting eq 10, i.e. as λ = Φ−1(1 − α); otherwise, e.g. λ = 1.358 in the case α = 0.05, (iii) determine the distance between F*N (τ) and F(τ) as DN = λ/√N (obtained by substituting the inequality √QDQ > λ in eq 10 with the equality, and using Q = N), e.g. DN = 0.0916 in the case α = 0.05, and (iv) use it to draw a band around the eCDF, which is defined as F*N (τ) ± DN and is shown in Figure 3 in light yellow color. The interpretation of this first test is that under the null hypothesis the continuous cumulative distribution function F(τ) will lie within the light yellow band of Figure 3 with a probability P = 1 − α. It is worth making a couple of additional remarks. First, the amplitude of such uncertainty band depends only on the number of data points, N, and on the significance level, α. Second, in the example of Figure 3, all eCDFs appear to lie within the yellow band, which is defined for α = 0.05 in this case. Could this feature be considered proof of good reproducibility? In principle, yes, but with the caveat that a different (larger) value of α would yield a different (narrower) band and hence a different assessment of reproducibility. To deal with this conundrum, the second statistical test is useful. 3.4. Combining Series. The second test, a 2-KS test, analyzes a pair of independent eCDFs F*N (τ) and F*M*(τ) and checks the null hypothesis that they are sampled independently from the same continuous cumulative distribution function.13 To this aim, one must (i) define Q = NM/(N + M), (ii) calculate DN,M through eq 7, (iii) compute the
Figure 5. Statistical analysis and combination of experimental series (at S = 1.41) with different numbers of points N (see Section 2.5). (a) Results of the KS-tests illustrated with respect to the Kolmogorov distribution. Crosses, triangles, and boxes have vertical coordinates corresponding to the KS-number for pairs, triples, and quadruples, respectively, and horizontal coordinate as discussed in the main text. The boundaries of the three shaded areas are defined by the KS-number, λ (vertical coordinate), and the confidence level, P (horizontal coordinate), associated with the three combined series: yellow for the pair, blue for the triple, and red for the quadruple. The staircase line gives for the series 3, pair 2−3, triple 1−2−3, and quadruple 1−2−3−4 the corresponding number of data points. (b) The four curves plot the function Ej(P) defined in the interval [Pj,1), as defined in the main text; colors are black, yellow, blue, and red, for the series 3, 2−3, 1−2−3, and 1−2−3−4, respectively. Φ(λ) shown in Figure 4. Their coordinates are the KS-number itself (λ coordinate) and for the sake of comparison and clarity the value of P associated with the smallest KS-number among them (yellow cross), which is P = 0.271 for the pair 2−3. Because uncertainty, e.g., the uncertainty band in Figure 3, scales with 1/√N, the reciprocal of the square root of the number of data points, combining experimental series reduces uncertainty. Combining the data points of series 2 and 3 into one, thus forming what we call the pair 2−3 with a total number of N23 = 338 data points, constitutes the safest, most conservative choice. This is also the one that imposes the strictest constraint on reproducibility between series, E
DOI: 10.1021/acs.cgd.7b01014 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
(see the beginning of Section 3.3) is obtained by dividing λ by the square root of a different number, namely Mj; we call such uncertainty Ej. We now call Pj the minimum confidence level at which the jth eCDF is feasible, and λj the corresponding K-number from eq 9 (upper border of the shaded regions in Figure 5a); note that P1 = λ1 = 0. The uncertainty associated with the jth combined eCDF is a function of the confidence level as calculated above or equivalently of λ; it defines a curve in the open interval [Pj, 1), which is plotted in Figure 5b and is called Ej(P). It is worth noting that the smallest uncertainty levels for the three combined series, i.e. Ej(Pj), are rather similar, although they correspond to very different confidence levels, whereas at a given confidence level, e.g. P̃ = 0.95, where all combinations are feasible, the uncertainty Ej(P̃ ) decreases for increasing number of data points, as expected. This effect highlights once again the importance of considering the whole spectrum of KS-number and confidence levels when considering how to combine single eCDFs into a larger set of data points as a basis for the estimation of physical properties such as the nucleation rate.
namely a confidence level of 0.271 and a significance level of 0.739, which are rather far from the typical values of 0.95 and 0.05, respectively, used in statistical tests. The yellow-shaded region in Figure 5a highlights the KS-number and the corresponding confidence level (0.271) associated with the pair 2−3. The next step in increasing the number of data points and decreasing the uncertainty is that of forming a triple, e.g. in this case the triple 1−2− 3. To this aim, one considers the KS-numbers of all pairs involved, i.e. 1− 2, 2−3 and 1−3, and those of the 2KS tests between these three pairs and the third single series, i.e. (1−2)-3, (2−3)-1, and (1−3)-2. All such ternary 2KS tests and the corresponding KS-numbers were used to plot the blue triangle symbols in Figure 5a, while the corresponding 2KS tests on pairs above are plotted as the two blue crosses and as the yellow cross in the same figure (the black triangles and the black crosses correspond to the other possible combinations involving also series 4). Let us consider the four possible triples that can be obtained from the four series in Figure 3, i.e. 1−2−3, 1−2−4, 1−3−4, and 2−3−4, as well as the corresponding six KS-numbers, which are associated with each triple as described above for the triple 1−2−3; among these, let us use the largest KS-number for each triple. In the case of the data in Figure 3, the triple 1−2−3 exhibits the smallest largest KS-number (0.91) among the four possible triples. As a consequence, the triple 1−2−3, having a total of N123 = 443 data points, is the safest, most conservative choice for a triple, and it was thus chosen as new combination. The blue-shaded region in Figure 5a highlights its KS-number and the corresponding confidence level (0.626). The triple 1−2−3 would fulfill the null hypothesis for all possible 2KS tests made with its three component series with a confidence level of 0.626 or larger (hence with a significance level of 0.374 or smaller). Finally, all four series in Figure 3 can be combined after evaluating the KS-numbers of all pairs, of all merged pairs with singles series, of merged pairs with other merged pairs (tests involving four original series overall), and of all triples with single series (involving four series as well). All KS-numbers obtained by including four series corresponds to the box symbols in Figure 5a. In this particular case, it turns out that the largest KS number of all is that of the 1−4 pair (1.19, corresponding to the red cross in Figure 5a). Such value determines the red-shaded region in Figure 5a corresponding to the quadruple 1−2−3−4 (defined by a total of N1234 = 527 data points), bounded by λ = 1.19 and P = 0.881. With reference to the quadruple 1−2−3−4 (red region), the triple 1− 2−3 (blue region), and the pair 2−3 (yellow region), Figure 5a illustrates effectively where each of these combinations of single experimental series is feasible, i.e. for decreasing values of KS-number, λ, and of confidence level, P, in going from the quadruple to the pair. It is obvious that in choosing how many individual series to combine, one has to trade off inclusiveness (in going from left to right in the diagram, more series are included while the confidence level increases and the significance level decreases) and strictness (in going from right to left, the λ value decreases, and the requirements for experimental series to be considered as sampling the same continuous cumulative distribution function become more and more stringent). Contrary to what typically is done in the literature, where fixed values of P and α are considered (typically 0.95 and 0.05, respectively) and based on these values, the null hypothesis is tested, we argue that Figure 5a, where the whole range of values of P and α are explored, provides a more complete picture of the statistical features of the different experimental series and their interrelation. Figure 5a shows also the number of data points on which the three combinations above, as well as the single series with the largest number of data points, are based. Using these numbers, the uncertainty associated with each of these four series (the single series 3, the pair 2−3, the triple 1−2−3, and the quadruple 1−2−3−4) can be calculated for any confidence level larger than the smallest confidence level, at which that series is feasible (left border of the shaded regions in Figure 5a). We labeled each of these four eCDFs with the index j, chosen as the number of original series included: the series 3 has j = 1, the pair 2−3 has j = 2, the triple 1−2−3 has j = 3, and the quadruple has j = 4. Accordingly, their number of data points is called Mj, so that M1 = N3, M2 = N23, M3 = N123, and M4 = N1234. Because they have different numbers of data points, the uncertainty with which the underlying continuous cumulative distribution function is estimated
4. ESTIMATION OF NUCLEATION RATE 4.1. Nucleation Rates from Combined eCDFs. The ultimate goal of this work was to determine an estimate of the nucleation rate from the eCDF at a given supersaturation through the procedure described in Section 2.5 and assess its accuracy using the statistical tools presented in the previous sections. The most important result of the previous section is that such estimation can be based on individual experimental series or on different combinations thereof, each combination being feasible beyond a minimum confidence level and corresponding to a different uncertainty (depending on the number of data points, based on which it is constructed). The goal of this section is to link the statistical features of the eCDF to the uncertainty to be attributed to the corresponding estimate of the nucleation rate. Let us consider again the four eCDFs measured at S = 1.41 and plotted in Figure 3 as well as all the eCDFs that can be formed by combining them in all possible ways. In Figure 6, the values of the
Figure 6. Nucleation rates estimated from the four experimental series at S = 1.41. Symbols refer to all individual series (•), and all possible combinations of two (×), three (▲), or four (■) series; their horizontal coordinate is the number of data points of the individual series of or the combines ones. Red error bars represent the asymmetric confidence intervals (eq 16, exact form) estimated at P̃ = 0.95, whereas the blue error bars represent those at P = Pj, computed using eqs 17 and 18. nucleation rate estimated by applying the method described in Section 2.5 to the four individual eCDFs (circles), the six possible pairs (crosses), the four triple (triangles), and the quadruple (box symbol) are plotted vs the number of data points forming that eCDF. It is immediately apparent how different the estimates obtained from the individual series are, namely from about 140 #/(s m3) to about 240 #/(s m3), despite the relative similarity of the four eCDFs in Figure 3. As series are combined, the variation between estimates from pairs and F
DOI: 10.1021/acs.cgd.7b01014 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
Figure 7. Propagation of uncertainty and confidence intervals for nucleation rate estimates. (a) Plot of the function ΔF*(η) defined by eq 16. For ΔF* ± DN, the curve gives upper and lower relative bounds of the confidence interval on the nucleation rate estimates, η+ and η− (see main text). (b) Plot of the (relative) confidence interval, η+ and η−, vs number of data points, N, at two confidence levels P̃ = 0.95 (see eq 17, solid line) and a generic value P = 0.75 (dashed line). Note that the confidence interval is always asymmetric (|η+| > |η−|) and that they agree qualitatively and semiquantitatively with those calculated by Xiao et al.14 and J. The relative change in the value of the nucleation rate, η, is defined as
from triples becomes smaller and smaller, whereas of course there is only one quadruple and hence only one nucleation rate estimate that is based on all the data points of all the four original eCDFs. It is thus obvious that if we had to select one of the single series, even that with the largest number of points, we would risk a large error. But this is also true if we decided to use a pair (or a triple) of experimental series to retain only the points belonging to the two series whose KSnumber is the smallest, i.e. the pair 2−3 associated with the cross at the rightmost vertex of the yellow-shaded triangle in Figure 6 (in fact the difference between the nucleation rate estimated from the pair 2−3 and that estimated from the pair 1−4 is about 50% of the former value). On the other hand, the variability exhibited by the estimates obtained from the single series might reflect an experimental variability that is intrinsic and that for the sake of best experimental practices should be accounted for. Based on this argument, one should combine all data points and estimate the nucleation rate from the quadruple. By doing so, however, one risks including in the analysis one or more real outliers caused by an undetected experimental error, e.g. the single series yielding the largest nucleation rate value, thus spoiling the final estimate of the nucleation rate. The discussion above demonstrates the existence of a real issue to be addressed by the researcher. It appears natural to exploit the statistical analysis outlined in the previous section to provide a rationale for the choice of the experimental series that have to be included in the combined eCDF, on which the final nucleation rate estimate is based. 4.2. Propagation of Uncertainty. This section aims at establishing a quantitative relationship between the intrinsic uncertainty of the eCDF, measured by the distance DN (see eq 6) and related to its number of data points N, and the uncertainty that characterizes the nucleation rate estimated from the eCDF via the regression presented in Section 2.5. Recently, such relationship has reportedly been determined through a large number of numerical calculations, yielding a 0.95 confidence interval of J as a function of N (see Figure 2 of Xiao et al.14). Assuming a general validity of the result thus obtained, the authors directly applied such numerical confidence interval to their experimental data. Here, we applied a different approach based on the following observations. On the one hand, any uncertainty on the eCDF propagates to a related uncertainty on the value of J extracted from that eCDF. On the other hand, any change in the value of J causes a related change in the cumulative distribution function given for instance by eq 4. The latter effect is easy to compute, as shown below, and can be used to determine the former. Let us consider two cumulative Poisson distributions, F(τ; J′) and F(τ; J), calculated for the same system (same system volume, V, and same growth time tG) but using two different nucleation rates, namely J′
η=
J′ − J J
(13)
Using eq 4, the difference between the two distributions (defined for τ ≥ tG) is ΔF(τ ) = exp(− JV (τ − tG)) − exp(− J ′V (τ − tG))
(14)
This function is 0 when τ = tG and τ →∞, and it is positive or negative if J′ > J or J′ < J, respectively. The stationary point of ΔF(τ) (maximum if positive, minimum if negative) is where d ΔF(τ)/dτ = 0, namely where JV (τ − tG) = log J ′ − log J
(15)
Substituting this into eq 14 yields the following equation for the stationary value of ΔF(τ), denoted ΔF*: η η ΔF * = ≃ 1 + 1/ η e (1 + η) (16) where e is the Euler number. The functional relationship between ΔF* and η is strictly monotonically increasing. The approximated linear relationship holds for small values of η (i.e., ΔF*/η approaches 1/e when η →0), causing an error on the value of ΔF* of up to 10% for |η| ≤ 0.2. The Kolmogorov statistics establishes the intrinsic uncertainty, DN, with which the eCDF F*N (τ) based on N data points estimates the underlying continuous cumulative distribution function for any given confidence level P (see eq 9). More precisely, according to the first test presented in Section 3.3, such uncertainty is calculated from the Knumber associated with P, i.e. λ obtained from eq 9, as DN = λ/√N. Such uncertainty on the underlying continuous cumulative distribution function corresponds to an uncertainty on the associated nucleation rate, and such correspondence must exactly be given by eq 16, where ΔF* = ±DN. Kolmogorov statistics applies in fact to a distance calculated using the infinite norm of eq 6. Note that eq 16 indicates that the relative uncertainty on the nucleation rate is about 2.7 times larger than that on the eCDF (in terms of absolute values). The exact formula yields asymmetric confidence intervals for the nucleation rate estimate with the distance of the upper bound from the estimated value larger than that of the lower bound. Hence, for the same absolute value of the left-hand side of eq 16 but different sign, one obtains two values of η, namely η+ and η−, with η+ > 0 > η− and |η+| > |η−|. The asymmetric behavior of η can be seen in Figure 7a, where η is plotted as a function of ΔF*, and in Figure 7b, where η+ and η− are reported as a function of N for two different confidence levels P̃ = 0.95 (solid line) and P = 0.75 (dashed line). Note that the relative uncertainty G
DOI: 10.1021/acs.cgd.7b01014 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
Figure 8. Statistical analysis and combination of series at different supersaturation levels. These four plots are obtained and have the same meaning as the plot of Figure 5b (which is reproduced here as Figure 8a for the sake of comparison with the others). (1 + ηj−̃ )Jj ≤ Jj ≤ (1 + ηj+̃ )Jj (j = 1, 2, 3, 4)
on the nucleation rate is still rather large even for a large number of points (|η| ≈ 0.2 for both confidence levels when N = 500, as can be seen in Figure 7b); the shape of the curves in Figure 7b suggest that further improving the accuracy of the nucleation rate estimate requires a number of experiments which might be experimentally challenging to obtain. These results are consistent both qualitatively and semiquantitatively with those reported by Xiao et al.14 (see Figure 2 of their work). Thus, P confidence intervals on the value of J estimated from a specific eCDF as in Section 2.5 are given by the following inequalities (1 + η−) J ≤ J ≤ (1 + η+) J; if and only if the linearized formula is used and then an approximate symmetric confidence interval is obtained. Finally, although eq 16 is valid for a Poisson process only, the logic followed here could be applied to any stochastic process, whatever its underlying statistics. 4.3. Estimating Nucleation Rates and Their Confidence Intervals. In this section, we consider the different possible estimates of the nucleation rate reported in Figure 6 and attribute to each of them proper confidence intervals. We focus on the four combined eCDFs highlighted in Figure 5b and the corresponding nucleation rate estimates plotted in Figure 6 (the rightmost circle, the rightmost cross, the triangle symbol at the right vertex of the blue shaded region, and the square symbol; see also Table 2 for the quantitative values). Consistent with the notation introduced in Section 3.3, we call such estimates Jj with j = 1, ..., 4. We consider first P̃ = 0.95 because of its broad use in the context of parameter estimation; the associated K-number is λ̃ = 1.358. For the jth combined eCDF, the left-hand side of eq 16 is calculated as ΔF* = ± Ej(P̃ ) ; by solving this equation, one obtains the two values η̃+j and η̃−j and hence the following confidence interval:
(17)
Such confidence intervals are indicated as red error bars in Figure 6, while the linear approximation is reported in Table 2, last column, as relative error. It is obvious that the amplitude of the confidence interval at a fixed value of the confidence level is controlled by the number of data points on which the nucleation rate estimate is based. Therefore, such amplitude increases when going from the quadruple 1−2−3−4 to the single series 3, namely from 0.34 to 0.56 in relative terms. In the case of the single series, such amplitude is quantitatively similar to that calculated as the difference between the maximum and minimum nucleation rates estimated from the 4 single series divided by the average value (about 0.53). While evaluating the 0.95 confidence intervals is the standard approach, in this work, we consider also the confidence intervals calculated at the minimum confidence level, at which each combined eCDF is feasible (such approach is meaningless in the case of a single series, where such minimum confidence level is zero). In other words, for the jth eCDF (j > 1), the left-hand side of eq 16 is calculated as ΔF* = ±Ej(Pj) (j > 1), where Pj is the mentioned minimum confidence level, which is highlighted in Figure 5b and reported in the fifth column of Table 2. By solving eq 16, one obtains the two values η+j and η−j , using which the following confidence interval is obtained:
(1 + ηj−)Jj ≤ Jj ≤ (1 + ηj+)Jj (j = 2, 3, 4)
(18)
These confidence intervals are indicated as blue error bars in Figure 6, whereas their linear approximations are reported in Table 2, second to last column, as relative errors. It is apparent that the amplitude of these confidence intervals (blue error bars) is smaller than that of the previous ones (red error bars) because the former are calculated at a lower H
DOI: 10.1021/acs.cgd.7b01014 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
confidence level than the latter. It is also evident that, differently from what was observed in the previous approach, the relative amplitude increases when going from the single series 3 to the quadruple 1−2−3− 4, namely from 0.22 to 0.30, and such an increase is less pronounced. It is worth noting that in the case of the combined eCDF consisting of two series, for instance the combination 2−3, the blue error bar is half of the red one. Before discussing these results any further, it is worth underlining that we analyzed specifically the 4 experimental series at supersaturation 1.41 but without any loss of generality. We in fact observed exactly the same trends with the experimental series measured at all the other supersaturation levels (see Section 5 below). 4.4. Choice of a Nucleation Rate Estimate and of Its Confidence Interval. Let us first underline that the statistical meaning of the confidence interval on an estimated quantity at a given confidence level P is the following: if the same statistical test was repeated many times (in our case, if the series 3 was repeated) at the same experimental conditions performing the same number of experiments and the P confidence intervals were calculated at every repetition, such intervals would bracket the true value of the estimated quantity in a fraction P of the repetitions.29 In looser terms, this definition can be rephrased as follows: the confidence interval of the nucleation rate estimate obtained in a single repetition brackets its true value with probability P. For the sake of brevity, we use the latter formulation in the following. Because the amplitude of the confidence interval increases with increasing confidence level, then reducing both reduces also the range of values attained by the physical properties that depend on the estimated nucleation rate (e.g., the amount of fines in a crystallization process) when the nucleation rate is taken within its confidence intervals. Such a reduction is obtained at the price of a reduced probability that the true value of the nucleation rate belongs to the considered confidence interval and that the corresponding physical property belongs to the associated range of values obtained. In more quantitative terms, because the amplitude of the confidence interval for a given eCDF scales with λ and with the reciprocal of the square root of its number of data points, to halve the 0.95 confidence interval (λ = 1.358), the confidence level has to be reduced to 0.20 (λ = 0.645), whatever the underlying number of data points. On the basis of the discussion in the previous section, it is clear that when considering a single series (series 3 in our case, for which j = 1) confidence intervals on the estimated nucleation rate for any confidence level P between 0 and 1 can be calculated. The choice of the relevant confidence level P is determined exclusively by the trade off just described between inclusiveness (high values of P) and strictness in constraining experimental variability (low values of P). The lowest value of P is 0; hence, the smallest amplitude is also 0, and therefore, only the confidence interval at P̃ = 0.95 is plotted in Figure 6 (red error bar). It is obvious that larger values of P could be considered but also that letting P approach 1 would be meaningless considering that the confidence interval would become infinitely large. In the cases of combined eCDFs instead (j = 2,3,4), there is a minimum confidence level Pj, below which the jth combination is not feasible. This yields a minimum amplitude of the confidence interval (blue error bar in Figure 6); between the blue and the red error bars in Figure 6 lie all possible confidence intervals for the nucleation rate Jj calculated at confidence levels between Pj and P̃ = 0.95. Again, the choice of the relevant confidence level P is determined by the trade off just described, but in this case there is a lower bound that cannot be overcome to avoid loss of feasibility of the considered combined eCDF.
As one can see in Figure 8, for each supersaturation, our method identifies the minimum confidence level at which two, three, four, or five series can be combined; such levels are different for different combinations. Furthermore, such minimum levels vary widely for different supersaturations: they are evenly spaced for S = 1.41 and 1.71 but differ strongly from one another for S = 1.58 and 1.91. As discussed in Section 3.4, the spacing of Pj provides an immediate, qualitative insight on how close and similar the experimental series are with each other. Another important observation concerns the feasibility of combining all series at the standard confidence level P̃ = 0.95. Indeed, for the results reported in this work, only S = 1.41 could form a combination with all series because for its quadruple P4 < P̃; for the remaining three supersaturations, the combination of all distributions can be obtained only for a confidence level larger than the reference value. We recall that the choice P̃ = 0.95, albeit commonly used and accepted, bears in itself no specific, privileged statistical meaning. This was originally pointed out already by Fisher,30 who introduced such concept and recommended to choose the confidence level on a case by case basis based on the knowledge one has of the specific experimental system. Our approach complies with two different requirements: it allows not only a sound statistical analysis of the data but also clear and thorough detail of all the steps leading to the choice of whether to combine different experimental series, thus enhancing the transparency of the results and their portability and verifiability. Finally, we highlight a practical point: collecting data at high supersaturation is cumbersome because often experiments are rejected by the criteria on temperature and concentration, as discussed in Section 2.2 (see Table 1). In fact, at high Table 1. Effect of Supersaturation Level on the Relationship between the Number of Experiments and the Number of Acceptable Experiments as Defined in Section 2.2.3a S
series
Ntot
N
N/Ntot
1.41 1.58 1.71 1.91
1−2−3−4 2−3−4 1−2−3−5 2−3−4−5
704 896 800 1,072
527 728 481 507
0.74 0.81 0.60 0.47
All detection times were measured at 25 ± 0.5 °C (see Section 2.2). At each supersaturation, the largest combined series feasible at P̃ = 0.95 was considered. a
supersaturation, crystals can form more likely during the cooling phase before the system attains the desired crystallization temperature. Investigating such conditions would thus require several independent repetitions, further stressing the importance of a clear procedure to combine different series of experiments together. 5.2. Nucleation Rates. After performing the statistical analysis of the experimental series and combining them accordingly, we estimated for each supersaturation the values of the nucleation rate associated with each series (individual and combined), Jj; these are plotted in Figure 9 vs the confidence level Pj, corresponding to the individual or combined eCDF. Note that the choice of the linear scale for J penalizes the visibility of the data at S = 1.41, which are clearly illustrated in Figure 6. In light of the results in Section 5.1, we also computed the inherent uncertainty of each nucleation rate, reported in Figure 9 as red and blue error bars estimated according to eqs 17 and 18, respectively (see also Table 2 for the exact values).
5. RESULTS 5.1. Supersaturation. The statistical analysis discussed in Sections 3 and 4 was applied to all values of supersaturation that we investigated, namely S = 1.58, 1.71, and 1.91, beside S = 1.41; either four or five experimental series were carried out for each supersaturation level, yielding results qualitatively comparable to those shown in Figure 3. Here, we illustrate and discuss the outcome of the statistical analysis by analyzing the experiments in terms of diagrams created, as in Figure 5b. I
DOI: 10.1021/acs.cgd.7b01014 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
S = 1.91 because these combinations could not be formed at the reference confidence level P̃ = 0.95, as highlighted in Section 5.1. Note that for all values of supersaturations investigated, the amplitude of the confidence interval of the nucleation rate decreases by combining more series together because the number of points Mj increases (see Table 2). However, by comparing the uncertainties at different supersaturations, the data appear to be heteroscedastic (i.e., their variance is not constant but depends on the supersaturation) because such amplitude increases with S; this behavior was reported also by Xiao et al.14 The values of J obtained so far are point estimates, i.e. each value of S yields one value of J. However, we would rather be able to fit a functional form dependent on supersaturation and temperature, which could also be used for modeling and predicting the behavior of the system at conditions not yet explored experimentally. Several alternative models have been discussed and reported in the literature to interpret detection time experiments,14,20,31 and we investigated this issue in two previous papers.12,19 A proper fitting can be based on classical least-squares methods or on criteria derived from the Kolmogorov−Smirnov statistics,12 but both approaches require nucleation and growth rate models, physical conditions for detection, and a proper understanding of secondary nucleation, its role, and mechanism. Finally, we are well aware that the system investigated here, likewise others presented in the literature,21,22 exhibit more polymorphic forms. While the nucleation rate estimate in terms of nuclei formed per unit time and volume has a clear meaning also in the case of polymorphic systems, its interpretation requires determination of which polymorph forms under which conditions and how stochasticity affects the formation of the different polymorphs. All this goes beyond the scope of this work and will be the subject of a follow-up work.
Figure 9. Nucleation rate at different values of S with their associated asymmetric confidence intervals computed according to eqs 17 and 18, represented as blue and red error bars, respectively. The bullets (•), crosses (×), triangles (▲), square (■), and stars (★) correspond to values of J calculated from the single distributions, pairs, triples, quadruple, and quintuples, respectively. For each value of S, the dashed line corresponds to the value of J of the distribution combined at the highest value of P. Nucleation rates are plotted vs the confidence level Pj of the corresponding combination of series. Because the combinations of all series at S ≥ 1.58 are feasible at a confidence level larger than P̃ = 0.95, the red error bars are missing in these cases.
Table 2. Summary of Results Obtained by Applying the Statistical Analysis and Procedure to Estimate Nucleation Rates and Confidence Intervals to the Series Measured at Different Supersaturationsa S (−)
series
Mj (#)
Jj [#/(s m3)]
Pj (−)
ηj (−)
η̃j (−)
1.41
3 2−3 1−2−3 1−2−3−4 2 3−2 2−3−4 1−2−3−4 4 1−2 1−2−3 1−2−3−5 1−2−3−4−5 2 3−5 2−3−5 2−3−4−5 1−2−3−4−5
220 338 422 527 316 412 709 1,025 154 261 373 481 635 260 213 389 649 767
167 157 167 170 628 711 667 709 920 953 973 1023 994 2424 1694 1516 1744 1689
0 0.27 0.63 0.88 0 0.16 0.75 0.98 0 0.18 0.57 0.76 0.99 0 0.01 0.75 0.93 0.99
0 0.11 0.13 0.15 0 0.09 0.11 0.14 0 0.13 0.14 0.14 0.20 0 0.08 0.13 0.15 0.18
0.28 0.22 0.20 0.17 0.23 0.20 0.15
1.58
1.71
1.91
6. CONCLUSIONS Estimating the rate of nucleation and the rate of a stochastic process from measurements presents challenges foreign to deterministic processes. Even though it is obvious that more experiments should result in a better estimate, provided such experiments are consistent with each other and carried out under exactly the same experimental conditions, it is also clear that any experimental technique will be intrinsically limited in the number of data points that can be collected in a single series of measurements. In this paper, we addressed the issue of when and how multiple experimental series, each yielding a so-called eCDF, can be merged to obtain a combined series with a number of data points given by the sum of those of each individual series. We investigated the problem starting from a well established statistical framework, namely the classical studies by Kolmogorov and Smirnov. We then presented a few novel contributions to the lively debate and growing literature on primary nucleation. First, we applied a modified version of a commercial tool to improve its experimental accuracy and enhance its productivity because rigorous control of experimental conditions is crucial. Then, we utilized Kolmogorov and Kolmogorov−Smirnov tests to assess under which conditions two individual series can be combined and, building on this analysis, discussed how an arbitrary number of experimental series can be analyzed and possibly combined into one. While avoiding conducting the whole discussion in terms of standard values of the confidence level and related
0.35 0.26 0.21 0.18 0.26 0.29 0.21 0.16
a
The same individual and combined series selected in Figures 8 and 9 are considered. The corresponding number of points, Mj, and estimated nucleation rate, Jj, are given. Note that the confidence levels are calculated for the sake of simplicity using the linear approximation given by eq 16 (symmetric intervals); ηj is computed at the minimum feasible confidence level Pj (column 5) and η̃j at the standard confidence level P̃ = 0.95.
Even though one can arbitrarily choose a different confidence value P ≥ Pj, the uncertainties on the nucleation rates will not be smaller than those computed at Pj. No uncertainty is reported for the quadruple at S = 1.58 and for the quintuples at S = 1.71 and at J
DOI: 10.1021/acs.cgd.7b01014 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
(7) Yang, H.; Svärd, M.; Zeglinski, J.; Rasmuson, Å. C. Cryst. Growth Des. 2014, 14, 3890−3902. (8) Liu, J.; Svärd, M.; Rasmuson, A. C. Cryst. Growth Des. 2014, 14, 5521−5531. (9) Brandel, C.; ter Horst, J. H. Faraday Discuss. 2015, 179, 199−214. (10) Pons Siepermann, C. A.; Huang, S.; Myerson, A. S. Cryst. Growth Des. 2017, 17, 2646−2653. (11) Sear, R. P. CrystEngComm 2014, 16, 6506−6522. (12) Maggioni, G. M.; Mazzotti, M. Cryst. Growth Des. 2017, 17, 3625− 3635. (13) Little, L. J.; Sear, R. P.; Keddie, J. L. Cryst. Growth Des. 2015, 15, 5345−5354. (14) Xiao, Y.; Tang, S. K.; Hao, H.; Davey, R. J.; Vetter, T. Cryst. Growth Des. 2017, 17, 2852−2863. (15) Forsyth, C.; Mulheran, P. A.; Forsyth, C.; Haw, M. D.; Burns, I. S.; Sefcik, J. Cryst. Growth Des. 2015, 15, 94−102. (16) Forsyth, C.; Burns, I. S.; Mulheran, P. A.; Sefcik, J. Cryst. Growth Des. 2016, 16, 136−144. (17) Kadam, S. S.; Kulkarni, S. A.; Ribera, R. C.; Stankiewicz, A. I.; ter Horst, J. H.; Kramer, H. J. Chem. Eng. Sci. 2012, 72, 10−19. (18) Kulkarni, S. A.; Kadam, S. S.; Meekes, H.; Stankiewicz, A. I.; Ter Horst, J. H. Cryst. Growth Des. 2013, 13, 2435−2440. (19) Maggioni, G. M.; Mazzotti, M. Faraday Discuss. 2015, 179, 359− 382. (20) Jiang, S.; ter Horst, J. H. Cryst. Growth Des. 2011, 11, 256−261. (21) Hansen, T. B.; Taris, A.; Rong, B. G.; Grosso, M.; Qu, H. J. Cryst. Growth 2016, 450, 81−90. (22) Kulkarni, S. A.; McGarrity, E. S.; Meekes, H.; ter Horst, J. H. Chem. Commun. 2012, 48, 4983. (23) DeGroot, M. H.; Schervish, M. J.; Sheet, C. Probability and Statistics 2011, 911. (24) Kolmogoroff, A. Ann. Math. Stat. 1941, 12, 461−463. (25) Smirnov, N. Ann. Math. Stat. 1948, 19, 279−281. (26) Massey, F. J. J. Am. Stat. Assoc. 1951, 46, 68−78. (27) Dvoretzky, A.; Kiefer, J.; Wolfowitz, J. Ann. Math. Stat. 1956, 27, 642−669. (28) Massart, P. Ann. Probab. 1990, 18, 1269−1283. (29) Ross, S. M. Introduction to Probability and Statistics for Engineers and Scientists; Academic Press: Cambridge, 2014. (30) Fisher, R. A. Biol. Monogr. Manuals, 5th ed.; Oliver & Boyd: Edinburgh, 1925. (31) Sullivan, R.; Davey, R.; Sadiq, G.; Dent, G.; Back, K.; Ter Horst, J.; Toroz, D.; Hammond, R. Cryst. Growth Des. 2014, 14, 2689−2696.
concepts (e.g., 95% confidence level), we explored the whole available range of confidence levels to enhance our understanding. By doing so, we introduced a spectrum of the characteristics of the series that we would like to merge and showed the informative abilities of the tool. While we did not give any one-fits-all recommendation for the choice of a specific estimate of the nucleation rate from the several that were presented and discussed, we believe that the tools and insight provided here will help researchers make informed choices. In fact, because combining all nucleation time measurement might be misleading, the procedure used for combining individual series must be clearly documented: it is each researcher’s choice whether to combine the data together and at which confidence level, but use of the methods illustrated here will enhance transparency and clarity. We also discussed how to estimate nucleation rates from the combined series as well as the trade-off between inclusiveness in using as many original series as possible and strictness in limiting the acceptable experimental variability from one series to the other. Such analysis was heavily based on an exact relation for the propagation of the uncertainty from the eCDFs to the nucleation rates, which was derived here. Different functional forms would be obtained by applying statistics different from the Poisson distribution, but our methodology would still apply, and we believe our main conclusions would not be altered. We did not use the nucleation rates obtained at different supersaturation levels to estimate the parameters of a specific nucleation rate expression. We made this choice not because we consider this unimportant; in fact, we investigated this topic in previous works.12,19 However, we believe that the statistical analysis described and applied here is a prerequisite to estimate nucleation rates at specific operating conditions independent of the functional form of the nucleation rate on physical properties and operating conditions. Moreover, estimating the nucleation rate parameters requires making assumptions on or investigating more deeply crystal growth and detection methods and, for the system considered here, its polymorphism. With respect to all of these issues, the analysis developed in this contribution is indeed a prerequisite.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]; Phone: +41 44 632 24 56; Fax: +41 44 632 11 41. ORCID
Marco Mazzotti: 0000-0002-4948-6705 Author Contributions ‡
G.M.M. and L.B. contributed equally to the manuscript.
Notes
The authors declare no competing financial interest.
■
REFERENCES
(1) Becker, R.; Döring, W. Ann. Phys. 1935, 416, 719−752. (2) Kashchiev, D. Nucleation: Basic Theory with Application; Butterworth-Heinemann: Oxford, 2000. (3) Toschev, S.; Milchev, A.; Stoyanov, S. J. Cryst. Growth 1972, 13− 14, 123−127. (4) Pino-García, O.; Rasmuson, Å. C. Ind. Eng. Chem. Res. 2003, 42, 4899−4909. (5) Goh, L.; Chen, K.; Bhamidi, V.; He, G.; Kee, N. C. S.; Kenis, P. J. A.; Zukoski, C. F.; Braatz, R. Cryst. Growth Des. 2010, 10, 2515−2521. (6) Yang, H.; Rasmuson, Å. C. Cryst. Growth Des. 2013, 13, 4226− 4238. K
DOI: 10.1021/acs.cgd.7b01014 Cryst. Growth Des. XXXX, XXX, XXX−XXX