A Population Balance Model for Chiral Resolution via Viedma

ARTICLE pubs.acs.org/crystal

A Population Balance Model for Chiral Resolution via Viedma Ripening Martin Iggland and Marco Mazzotti* Institute of Process Engineering, ETH Zurich, Sonneggstrasse 3, 8092 Zurich, Switzerland ABSTRACT: It has recently been observed that a mixture of two enantiomeric crystals subjected to attrition in a solution containing a racemizing agent undergoes symmetry breaking, with only one of the enantiomers eventually remaining in solid form. This process is believed to occur due to the interplay between racemization in solution, attrition, agglomeration and crystal growth and dissolution caused by the crystal size dependence of solubility. In this work, we present a mathematical model of this process based on population balances, incorporating all these elementary mechanisms. This model is able to explain the experimentally observed behavior. A sensitivity analysis on several parameters is carried out, leading to the conclusion that racemization, growth and dissolution, and agglomeration are necessary for deracemization to occur. While deracemization occurs also in the absence of attrition phenomena, these enhance it by reducing the time needed to achieve complete chiral purity.

’ INTRODUCTION Chirality plays a paramount role in nature, besides being an important issue in the pharmaceutical industry, since two enantiomers of the same chiral substance can have vastly different effects, for instance when interacting with biological receptors in the human body. Enantiomers can be obtained in pure form by either chiral synthesis or enantiomeric resolution, or by a combination of both. Chiral resolution can be achieved by chromatographic methods or by preferential crystallization, where one enantiomer is allowed to crystallize while the other should remain in solution. Sodium chlorate (NaClO3) is an achiral molecule which can crystallize in two enantiomeric forms. Recently, Viedma1 observed experimentally the symmetry breaking of a mixture of both enantiomeric forms of crystalline sodium chlorate to only one enantiomer, achieved simply by stirring their aqueous suspension in the presence of glass beads. Noorduin et al.2,3 have observed the same behavior in a system containing a chiral molecule (N-(2-methylbentylidene)-phenylglycine amide) and a racemizing agent (1,8-diazabicyclo[5.4.0]undec-7-ene (DBU)). Deracemization of the solid phase has also been observed in a system where the racemization in solution takes place in two steps.4,5 It has been proposed that this mechanism of deracemization be given the name Viedma ripening.6 The parameters affecting the rate of change of the enantiomeric excess in the solid phase have been investigated experimentally. It has been found that a small enantiomeric excess, or a size imbalance, is sufficient to direct the system toward the chiral form initially in excess. Increased attrition by addition of glass beads or an increased stirring rate leads to a decrease in the time required to achieve complete chiral conversion,1,3,7
9 as does an increased racemization rate.3 It has been reported that the evolution can be directed toward a specific enantiomer using chiral additives.2,10 r 2011 American Chemical Society

All reported examples are for systems forming conglomerates, that is, where no solid racemic compound forms. The process has also been demonstrated to work on a system exhibiting epitaxial growth.11 This fact, along with the requirement of solution-phase racemization, limits the number of systems to which Viedma ripening could be applied industrially. Only a fraction of all pharmaceutically relevant compounds form conglomerates. A variation of the process considered here is to start with a clear solution and induce nucleation by, for example, cooling. Compared to the corresponding process without attritionenhancement and racemizing agents (i.e., standard preferential crystallization techniques), this will lead to a pure solid product more often, with reduced sensitivity to experimental parameters.12
14 Since this involves nucleation, the mechanisms may be quite different from those considered in this work. The mechanism of this symmetry breaking process has been subject of intense debate recently, and there is still uncertainty about the details, as recently discussed in the literature.15 The mechanism we assume in this work is based on the propositions of several authors1,6,16,17 and includes racemization in solution, a size-dependent solubility, agglomeration and attrition. The influence of these mechanisms on particle size and solution concentration is shown schematically in Figure 1a. Racemization in solution is due to the action of a catalyst and is responsible for the conversion of solute molecules of one enantiomer into the other. The case of achiral molecules yielding enantiomeric crystals is equivalent to that of a chiral compound undergoing an infinitely fast racemization in solution. The size-dependence of solubility is a thermodynamic effect, which is also responsible for the process known as Ostwald Received: July 8, 2011 Revised: August 12, 2011 Published: September 19, 2011 4611

dx.doi.org/10.1021/cg2008599 | Cryst. Growth Des. 2011, 11, 4611–4622

Crystal Growth & Design

ARTICLE

Figure 1. Process scheme showing actions of various mechanisms. (a) Agglomeration between crystals increases their size and reduces the number of crystals, while attrition leads to the creation of new, small crystals. Growth and dissolution affect the size of crystals and are responsible for transfer of solute material between phases. (b) The black curve shows the solubility as a function of size, as given by eq 2. The gray lines show the effects of the various mechanisms on one particle: agglomeration and attrition affect only the size of crystals, while growth and dissolution affect both the size of crystals and the concentration in solution.

ripening. Crystal formation from a supersaturated solution is associated with a Gibbs free energy gain proportional to the volume of the newly formed crystal and to a Gibbs free energy loss associated with the creation of its surface. Since the surface dominates at small sizes, whereas the volume dominates at large sizes, the energy gain balances the loss at a critical size xC, beyond which new nuclei are stable and can grow.18 The critical size depends on the supersaturation in solution S = c/c∞, according to19 xC ðSÞ ¼

2ka γVm α ¼ ln S 3kv RT ln S

ð1Þ

where c is the solute concentration in concentration and c∞ is the bulk solubility, that is, the solubility of an infinitely large particle. The capillary length α combines all the physical parameters of the crystals (surface tension γ, molar volume Vm) and their surface and volume shape factors ka and kv; R and T are the universal gas constant and the absolute temperature. The fact that the stability of a crystal depends on its size implies that the solubility of a crystalline material is dependent on the crystal size, as can be seen by using eq 1 to obtain the solubility c*(x) of a particle of size x, that is, the Gibbs
Kelvin or Gibbs
Thomson equation:
 cðxÞ α ¼ exp ð2Þ SðxÞ ¼ c∞ x Equations 1 and 2 define the curve shown in Figure 1b, which gives the critical size as a function of concentration (eq 1) or the solubility of a particle as a function of its size (eq 2). Every point in the plane represents a crystal of size x suspended in a solution having concentration c = Sc∞. Points underneath the curve represent crystals that dissolve either because they are small or because the supersaturation is low. Points above the curve correspond to crystals that grow because they are large or because the solution is highly supersaturated. Points on the line correspond to crystals having exactly the critical size at that solution concentration and therefore being stable. Crystals with sizes below the critical size will dissolve, thus increasing the solute concentration in solution. In turn, crystals larger than the critical size will grow, thus consuming the supersaturation.

Attrition leads to a reduction in crystal size and the creation of small crystals having the same chirality as their parent crystal. Agglomeration of small crystals with larger crystals of the same chirality leads to an increase in size and a decrease in the total number of crystals. Neither attrition nor agglomeration affect the concentration in solution directly but can cause a crystal to jump between the dissolution and the growth regime. Since agglomeration involves two particles, it has a higher order dependence on particle concentration than growth or attrition. Agglomeration is much faster for the population whose crystals are present in a larger number.16 The relatively simple experimental protocol of this process gives it potential to become a valuable alternative method for chiral resolutions. In order to make control, optimization, and scale-up possible, a deeper understanding of the mechanisms and of the effects of process parameters is required. More importantly, a detailed mathematical model which can successfully reproduce experimental observations is a way to demonstrate that the mechanism is plausible. Several models have been presented earlier in the literature. Uwaha20,21 proposed a reaction-type model to explain the evolution toward a homochiral state, considering only the total concentration of crystals of each configuration, monomers in solution, and chiral clusters. This model does reproduce some experimental observations, but, due to its assumptions, its ability to describe size-dependent effects is very limited. Noorduin et al.6 have adapted this model in an attempt to explain experiments carried out without the racemization reaction. Saito and Hyuga22 have also presented a reaction-type model, including crystals of various sizes. However, the largest cluster is made up of only six monomer units. Wattis23 used these models, making various simplifications to be able to analytically investigate the model behavior. Noorduin et al.17 modeled the chiral conversion process using a Monte Carlo type model including attrition effects and an exchange of monomers between two crystals, representing Ostwald ripening. However, the critical size as a function of concentration is not considered. All monomer exchange takes place directly and not via the solution, even when there is a racemization step in between. Agglomeration is not included in the model. Another Monte Carlo type model was presented by Cartwright et al.24 In this model, crystals are not present initially but are created by primary nucleation. Attrition reduces the size of clusters, by creating new crystals of the same chirality as the parent crystal. Ostwald ripening is assumed to cause the dissolution only of crystals in the smallest size range. Uwaha and Katsuno25 use a Monte Carlo model in order to investigate whether Ostwald ripening is essential for attrition-enhanced chiral resolution. In their model, attrition is assumed to limit the size of crystals. Having reached this maximum size, the crystals can break, thus creating fragments of various sizes. Other mechanisms modeled include the reversible attachment of monomer and dimer units to clusters. This model has then been extended to cover an open system, in which some of the largest chiral clusters are taken out of the system and replaced by achiral molecules or by a racemic mixture of clusters.26 In another paper,27 the same authors consider a different description of Ostwald ripening, incorporating the Gibbs
Thomson effect for the size-dependence of solubility. Saito and Hyuga have used a series of Monte Carlo type models to simulate the interaction of monomers and clusters on surfaces,28
31 by also considering the action of grinding and the influence of the number of kink sites 4612

dx.doi.org/10.1021/cg2008599 |Cryst. Growth Des. 2011, 11, 4611–4622

Crystal Growth & Design

ARTICLE

on the attachment and detachment rates. Hatch et al.32 have also used a Monte Carlo simulation to model chiral symmetry breaking, focusing on diffusion on a surface and on the kinetics of the subsequent formation of chiral species. The assumption that growth and dissolution of particles due to size-dependent solubility leads to the direct exchange of substance molecules between crystals, without passing through the solution, is equivalent to assuming that growth and dissolution do not lead to any change in mass of the crystal population. These assumptions are not always justified when considering Ostwald ripening. Furthermore, several of the models mentioned neglect the concentration dependence of the critical size, which is the basis for size-dependent solubility in the first place. Attrition and agglomeration phenomena are also described in a very simplified manner in the models mentioned above. Attrition is assumed to limit the size of crystals to an arbitrarily chosen maximum size. Crystals that reach this size are broken into parts. Agglomeration is described by the incorporation of growth clusters, a cluster being a collection of a few, usually two, monomers. The models referenced above are able to handle systems having either chiral or achiral molecules in solution and are able to qualitatively reproduce the experimental observations. However, they do not accurately reflect the mechanisms which are proposed to be the cause of the observed behavior. In this paper, Viedma ripening is modeled using population balance equations (PBEs) for the first time. This is possible and necessary because all mechanisms involved in Viedma ripening are well-known and can be modeled rigorously in a population balance framework. We believe that this detailed model overcomes the assumptions and limitations of previous models. The objective is to demonstrate that the proposed mechanism is correct. In the first part of this paper, we discuss the individual mechanisms and present the mathematical model. In the second part, simulation results are presented and a sensitivity analysis on various kinetic parameters is carried out. Some conclusions are drawn in the final section.

’ MATHEMATICAL MODEL Mechanisms. The mechanisms accounted for by the mathematical model are racemization in solution, crystal attrition and agglomeration, as well as growth and dissolution driven by sizedependent solubility. We consider a perfectly symmetric system in terms of mechanisms and their rates, where the only asymmetry between the two enantiomers stems from the initial conditions. Accordingly, we have not considered the presence of chiral additives. Note that D and L are the two enantiomers of the chiral compound under consideration, having a concentration cD and cL in solution and that two populations of enantiopure crystals are present, characterized by particle size distributions nD and nL, where ni(x) dx is the concentration of particles of the population i (i = D, L) of characteristic size between x and x + dx. Racemization. Molecules in solution can change their stereochemical configuration due to the racemization reaction. In practical applications, racemizing agents such as DBU2,3,9 or sodium methoxide33 are used in order to promote this reaction. The reaction is first-order in both directions and proceeds according to kD

D h L, kL

kD ¼ kL ¼ kr

Figure 2. Agglomeration kernels A1 and A2 for typical kinetic parameters.

where kD and kL are kinetic rate constants. They must be equal for symmetry reasons, since the reaction is at equilibrium for a racemic mixture. If there is only one type of solute molecule, as in the case of NaClO31 or ethylediammonium sulfate,7 the racemization reaction can be seen as being infinitely fast. Growth and Dissolution of Particles Subject to Size-Dependent Solubility. Because of the competition between energy gained by creation of volume and energy lost due to creation of surface, the solubility of a particle depends on its size.18 The function for the size-dependent solubility of crystals is given by eq 2 but can be simplified as34,35
 α SðxÞ ¼ 1 þ ð4Þ x Because of the racemization reaction, the enantiomeric excess in solution will always be low. Therefore, we can assume that the bulk solubility c∞ does not depend on the concentrations of the two enantiomers, and it is equal to the solubility of the enantiopure crystals in a racemic solution. Of course, it is the same for both enantiomers. The effect of the size-dependent solubility is incorporated into the growth rate expression, giving an expression such as the following for surface-integration limited growth:
 ci α
1
ð5Þ Gi ðxÞ ¼ kg x c∞ This growth-rate expression can be used in a population balance model to describe Ostwald ripening.19,34
36 Note that kg and α are the same for both enantiomers, again for symmetry reasons. Agglomeration. When modeling agglomeration, only collisions between two crystals of the same chirality are considered; that is, we assume that crystals of different enantiomers cannot agglomerate. The rate of agglomeration depends on the collision frequency and on the sticking probabibility, and can be sizedependent.37
39 In this work, all the terms which do not depend on the solids concentration are summarized in the agglomeration kernel A, which is a function of the sizes of both agglomerating particles. Two slightly different forms of the agglomeration kernel are used:39
 x1 þ x2 3 a1, 2 ð6Þ A1 ðx1 , x2 Þ ¼ 2 1 þ a1, 1 Xeq ðx1 , x2 Þ
 x1 þ x2 3 A2 ðx1 , x2 Þ ¼ a2, 2 expð
a2, 1 Xeq ðx1 , x2 ÞÞ 2

ð3Þ 4613

ð7Þ

dx.doi.org/10.1021/cg2008599 |Cryst. Growth Des. 2011, 11, 4611–4622

Crystal Growth & Design

ARTICLE

following two PBEs: ∂ni ∂ðGi ðxÞni ÞÞ þ ∂x ∂t Z ∞ ¼ ni ðη, tÞbðηÞgðx, ηÞ dη
ni ðx, tÞbðxÞ x pffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z qffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2 x Aðξ, 3 x3
ξ3 Þ 3 þ n ðξ, tÞn ð x3
ξ3 , tÞ dξ i i 2 0 ðx3
ξ3 Þ2=3 Z ∞
ni ðx, tÞ Aðξ, xÞni ðξ, tÞ dξ, ði ¼ L, DÞ ð11Þ 0

Concentration changes in the solution are due to growth and dissolution, as well as to racemization. The mass balance for the solutes L and D are written as

Figure 3. Daughter distribution for q = 10.

ðx1 x2 Þ2 Xeq ðx1 , x2 Þ ¼ 2 x1 þ x22
x1 x2

ð8Þ

For both kernels, agglomeration is fastest for a combination of large and small particles. In eq 6, agglomeration is slowest for a combination of two small particles, while in eq 7 agglomeration is slowest for a combination of two large particles. In Figure 2, the agglomeration kernels are visualized for typical values of the kinetic parameters, chosen so that similar values of the kernels are obtained in the small size range. The parameters used in eqs 6 and 7 are the same for both crystal populations for symmetry reasons. Attrition. Attrition occurs due to collisions of crystals with other crystals, with the walls of the reaction vessel, with the stirrer, or with glass beads, if present. These collisions can result in different numbers of fragments having different sizes. The distribution of these fragments is modeled by defining a daughter distribution.38,40
42 The daughter distribution used in this work is U-shaped and symmetric in volume-space, thus yielding only two fragments and is characterized by a single adjustable, integer parameter q, which affects the curvature of the distribution: the larger the value of q, the larger the difference in size between the two resulting fragments. Figure 3 shows the daughter distribution, whose equation reads as !2q
2q þ 1 2 η3 2 3 x
ð9Þ gðx, ηÞ ¼ 3x ð2q þ 1Þ 3 η 2 where η and x are the sizes of the original particle and of the formed fragment, respectively. The rate of breakage by attrition of particles of size x is given by the following rate function:
β x ð10Þ bðxÞ ¼ kb x0 where kb and β are adjustable parameters and x0 is a scaling factor (see eq 21 below). All the parameters introduced in this section, that is, q, kb, β, and x0, are the same for both crystal populations.

ð12Þ

dμL, 3 dcL ¼
kv F þ kr ðcD
cL Þ dt dt

ð13Þ

where the time derivative of the third moment of the distribution ni, μi,3, is given by the following equation. Z ∞ dμi, 3 ¼3 x2 Gi ðxÞni ðx, tÞ dx, ði ¼ L, DÞ ð14Þ dt 0 The initial and boundary conditions are as follows (i = L,D): ni ðx, t ¼ 0Þ ¼ n0i ðxÞ

ð15Þ

ci ðt ¼ 0Þ ¼ c0i

ð16Þ Z

μi, 3 ðt ¼ 0Þ ¼ μ0i, 3 ¼

∞ 0

ni ð0, tÞ ¼ 0

x3 n0i ðxÞ dx

ð17Þ ð18Þ

where c0i and n0i are the initial concentration of solute i and initial distribution of its crystals, respectively. The initial concentrations c0L and c0D are assumed to be equal because of the racemization reaction. A derivation of the boundary condition (eq 18) and of (eq 14) is given in the appendix. The progress of the resolution is measured using the solid phase enantiomeric excess of D, which is defined as ee ¼

μ3, D
μ3, L μ3, D þ μ3, L

ð19Þ

Formulation of the Model Using Dimensionless Variables. In order to reduce the number of parameters in the model, the following dimensionless variables are introduced:

Si ¼

’ MODEL Model Equations. The evolution of the particle size distribution of each population of enantiopure crystals under the effects of ripening, attrition, and agglomeration is described by the

dμD, 3 dcD ¼
kv F þ kr ðcL
cD Þ dt dt

ci x t ,y ¼ ,τ ¼ x0 t0 c∞

New distributions fi are introduced, in terms of the dimensionless size y; since fi dy = ni dx, fi = x0ni. The parameters t0 and x0 are chosen in such a way that the growth rate is simplified and one of 4614

dx.doi.org/10.1021/cg2008599 |Cryst. Growth Des. 2011, 11, 4611–4622

Crystal Growth & Design

ARTICLE

the distributions (we have chosen D, but this choice is arbitrary) has its initial mean size at y = 1 x0 ð20Þ t0 ¼ kg μ0D, 1 x0 ¼ 0 μD, 0

rj, 2

kr x0 kb x0 , Kbg ¼ , rj, 1 ¼ aj, 1 x20 , kg kg

aj, 2 x40 kv F ¼ , ν¼ c∞ kg

ð22Þ

where j = 1,2 depending on the agglomeration kernel used. With these definitions, eq 11 can be cast as 

∂fi ∂ðGi ðyÞfi Þ þ ∂y ∂τ ! Z ∞ β β ¼ Kbg fi ðε, τÞε gðy, εÞ dε
y fi ðy, τÞ y

pffiffiffiffiffiffiffiffiffiffiffiffiffi ffiffiffiffiffiffiffiffiffiffiffiffiffi  Ad ð 3 y3
y3a , ya Þ p 3 f y3
y3a , τ fi ðya , τÞ dya i ðy3
y3a Þ2=3 0 Z ∞
fi ðy, τÞ Ad ðy, ya Þfi ðya , τÞ dya , ði ¼ D, LÞ ð23Þ 2Z y

parameter

ð21Þ

The other distribution, for crystals of form L, has an adjustable initial mean y0L. For simplicity, we introduce the following parameters: Krg ¼

Table 1. Values of the Parameters Which Were Kept Constant for All Simulations

y þ 2

0

where G*i(y) = Gi(yx0)/kg and the agglomeration kernel Ad is given by one of the following two equations:
 y1 þ y2 3 r1, 2 ð24Þ Ad, 1 ðy1 , y2 Þ ¼ 2 1 þ r1, 1 Yeq ðy1 , y2 Þ Ad, 2 ðy1 , y2 Þ ¼ r2, 2


 y1 þ y2 3 expð
r2, 1 Yeq ðy1 , y2 ÞÞ 2

ð25Þ

where Yeq ðy1 , y2 Þ ¼

ðy1 y2 Þ2 y21 þ y22
y1 y2

ð26Þ

The derivative of the third moment of population i is given by Z ∞ dμi, 3 dji, 3  ¼ 3x30 ð27Þ y2 Gi ðyÞfi dy ¼ x30 dτ dτ 0 Hence the mass balances of eqs 12 and 13 can be recast as djD, 3 dSD ¼
νx30 þ Krg ðSL
SD Þ dτ dτ

ð28Þ

djL, 3 dSL ¼
νx30 þ Krg ðSD
SL Þ dτ dτ

ð29Þ

The parameter t0 can be interpreted as the characteristic time of growth of crystals of size x0. The parameter ν represents the ratio between the density of the substance in the solid phase and its density in the liquid phase.

value

x0

1  10
5 m

q β

10 1

S0D

1

S0L

1

Details about the numerical method used to solve the model equations are provided in Appendix B.

’ RESULTS The model presented above has several parameters that can be varied. We have chosen to keep some constant, such as the parameter q in the daughter distribution function, the initial mean size x0, and the initial supersaturation for both enantiomers, SD and SL; moreover, the initial crystal size distributions are all normally distributed. The constant parameters and their values are summarized in Table 1. The values of all other parameters for all simulations are given in Table 2. The results are presented in three parts. In the next section, we present a base case simulation with a set of parameters leading to complete chiral resolution, which demonstrates the capabilities of the model to describe deracemization via Viedma ripening. Then we present simulations using the same parameters as the base case but with varying initial conditions; these simulations are compared qualitatively to experimental observations, in order to validate the model. Finally, a sensitivity analysis to establish the role of different kinetic parameters is presented. Base Case. We have identified a set of parameters which leads to complete chiral resolution, that is, complete deracemization. In this base case, the two initial distributions have the same width and the same mean size. There are initially more crystals of form D than of form L, with an initial solid phase enantiomeric excess ee0 = 2%. See Table 2, simulation #1 for the values of all other parameters. The evolution of the solid phase enantiomeric excess is shown in Figure 4a. It follows an exponential increase, as observed in several experiments1,7,9,11,43 and reaches 99% enantiomeric excess at a dimensionless time τ = 460. The evolution of the mass can be followed through the third moments, plotted in Figure 4b. Initially, all crystals of both types start to dissolve, thus increasing supersaturation above the initial value S0i = 1. Dissolution continues until the critical size has decreased to a value low enough to allow net growth of D crystals, while L crystals continue to dissolve. Eventually, all L crystals are dissolved and L is converted to D through racemization. The supersaturation profiles, shown in Figure 4c, confirm this analysis. The strong initial increase is caused by the dissolution of crystals; the concentrations of both enantiomers are kept nearly equal by the racemization reaction. In the next phase, the concentration increase due to net dissolution of L crystals and the concentration decrease due to net growth of D crystals are of equal magnitude, leading to a nearly constant supersaturation. When the deracemization nears completion, the dissolution can no longer replenish the liquid phase, and the supersaturation keeps decreasing slowly due to ripening of D crystals. From the distributions shown in Figure 5, it can be seen that there is a rapid initial change for both distributions. Attrition leads to the appearance of a second 4615

dx.doi.org/10.1021/cg2008599 |Cryst. Growth Des. 2011, 11, 4611–4622

Crystal Growth & Design

ARTICLE

Table 2. Initial Parameters for All Simulationsa ID 1

Krg

α/x0

Aj

rj,1

rj,2/[m3]

Kbg

y0L

10

10
3

A1

1010

10
5

10

10

10
5

10


3

μ03,D

μ03,L

σ0D

1

0.05

0.048

0.1

1

0.048

0.05

0.1

σ0L

ν

outcome

0.1

1

D

0.1

1

L

2

10

10

A1

10

3a,b

10

10
3

A1

1010

10
5

10

1

0.05

0.03,0.045

0.1

0.1

1

D

4

10

10
3

A1

1010

10
5

10

0.5

0.05

0.05

0.1

0.1

1

D

10

10
5

10

1.5

0.05

0.03

0.1

0.1

1

D


3

5

10

10

A1

10

6

10

10
3

A1

1010

10
5

10

1

0.05

0.05

0.2

0.1

1

D

7

10

10
3

A1

1010

10
5

10

1

0.05

0.05

0.1

0.1

1




8a
g 9a
d

0 10

10
3 10
3

A1 A1

1010 1010

10
5 10
5

10 10

0.5 1

0.04
0.07 0.05

0.04
0.01 0.03

0.1 0.1

0.1 0.1

1 1
104




10a
i

10
4
104

10
3

A1

1010

10
5

10

1

0.05

0.03

0.1

0.1

1

D

11a
k

10

10
3

A1

1010

10
5

10
5
105

1

0.05

0.03

0.1

0.1

1

D/


12

10

A1

1010

10
5

0

1

0.025

0.015

0.1

0.1

1

D

A1

10

10

10
5

10

1

0.05

0.03

0.1

0.1

1

D/


13a
i

10
3
8


10


3

D

10

10

13j
n

10

10
2
103

A1

1010

10
5

10

1

50

30

0.1

0.1

1

D/


14a
l

10

10
3

A1

1010

10
10
10
1

10

1

0.05

0.03

0.1

0.1

1

D/


15a
k

10

10
3

A2

110

5  10
15
10
4

10

1

0.05

0.03

0.1

0.1

1

D/


a

The outcome denotes the type of crystals still present at the end of the process; a dash denotes that complete deracemization does not take place. When two outcomes are given (e.g. D/
), see the text and figures for details. Aj refers to the agglomeration kernel used, as given in eqs 24 and 25. All parameters, except rj,2, are dimensionless.

Figure 4. Simulation #1, base case. Evolution in time of (a) solid phase enantiomeric excess, (b) third moments of both crystal populations, (c) solution concentration. In (c), the dashed line represents the concentration of the D enantiomer, which is slightly lower than that of the L enantiomer.

Figure 5. Simulation #1, base case. Distribution at different times for population (a) D and (b) L. (c) Mean and critical sizes as a function of time. The dashed line represents the mean particle size of L crystals.

peak at smaller sizes. For the initially larger population of crystals (D), the combined effects of dissolution, growth, agglomeration,

and continued breakage cause the two peaks to grow into each other. The smaller population (L) continues to dissolve. 4616

dx.doi.org/10.1021/cg2008599 |Cryst. Growth Des. 2011, 11, 4611–4622

Crystal Growth & Design If the masses of the two initial populations are reversed (simulation #2), the evolution of all properties is symmetrical to that shown here, and in the final state only crystals of the L enantiomer are present. If the value of the parameter q in the daughter distribution is changed, the results obtained are similar and are therefore neither shown nor discussed any further. Contrary to the simulation illustrated in Figures 4 and 5 where racemization in solution is active from the start, in many experiments reported in the literature, the suspension of crystals is subjected to grinding for a length of time before the racemizing agent is added. In order to rule out that this difference has any substantial effects, we have repeated the base case simulation, disabling the racemization reaction up to a time τ = 100. We have observed that the only difference is that the initial stage is longer; hence, the whole deracemization process takes longer. Model Validation. The model presented here is able to reproduce the experimental behavior reported in the literature. However, since kinetic data are not available for all mechanisms included in the model, all results in this paper are qualitative. In order to validate the model, simulations have been run using parameters which reflect different experimental conditions reported in the literature. The first few simulations presented here concern the direction of evolution of the system, that is, toward more D or more L crystals, depending on the initial conditions. The two initial populations of crystals can differ in mass, mean size, or width (standard deviation), or in a combination of these factors. Initial Enantiomeric Excess. Experiments have shown that a slight initial enantiomeric excess in the solid phase can direct the evolution of the system toward the major enantiomer.1,2,7,9 Our simulations confirm that chiral purity is obtained for initial solid phase enantiomeric excesses of 2% (base case, above) as well as 5.2% and 25% (see simulations #3a,b in Table 2), with the only difference being that the larger the initial enantiomeric excess, the sooner complete chiral purity is reached. Initial Size Imbalance. Additionally to the enantiomeric excess, the size of the crystals in the two initial populations plays an important role. If the initial masses of the two populations are equal, the system will evolve toward the population of crystals which has the larger initial mean size (simulation #4). Experiments have also been reported in which there is a difference in both the initial size and the initial mass.44 In particular, a large mass of small crystals of one enantiomer were subjected to Viedma ripening together with a smaller mass of larger crystals, and at the end of the process, only crystals of the enantiomer which had the larger initial mass remained. Of course, this result can be reproduced by our simulations, as shown in simulation #5, where D crystals initially in excess, but with a smaller mean size than the L crystals, prevail at the end of the process. Since there are two opposing effects involved in this case, a generalization of this result is not trivial. Unless one population is dominant with respect to both mass and crystal size, the direction of evolution is not clear a priori. In the simulation cited above, the population with the larger initial mass prevailed. However, there should be a critical point at which the larger initial mean size of the other crystal population dominates and causes the direction of evolution to change. Initial Width of the Particle Size Distribution. The direction of the evolution of chirality during Viedma ripening can also be determined by a larger variation in the particle size distribution of

ARTICLE

Figure 6. Liquid phase enantiomeric excess of L as a function of solid phase enantiomeric excess of D, from simulations without racemization (#8a
g). The final time of all simulations is equal. Lines are to guide the eye.

one of the initial populations (i.e., higher standard deviation of the initial particle size distribution). As shown by simulation #6, the system evolves toward a state of complete chiral purity containing only crystals of the population which is initially broader. Equal Initial Conditions. Viedma1 and Cheung et al.7 have observed that systems with equal initial masses of the two enantiopure crystals evolve toward chiral purity, with the end result seemingly being random. In a simulation starting with the exact same mass and mean size for both populations (simulation #7), they evolve in exactly the same way, and thus there is no change in the enantiomeric excess. This shows that a system exhibiting perfect initial symmetry will conserve such symmetry. We believe that the experimentally observed deracemization is a consequence of the fact that the system was asymmetric from the start, although to such a small extent that could not be measured nor controlled.45 Other sources of asymmetry during ripening could be the presence of chiral impurities at trace concentrations.46 No Racemization. Noorduin et al.6 have performed experiments without racemizing agents and measured the solutionphase enantiomeric excess as a function of the solid-phase enantiomeric excess. They found that, without racemization, the solution is enriched in the enantiomer which is less abundant in the solid phase. Our model is able to reproduce these results, as shown in Figure 6. These simulations, in which the racemization reaction was suppressed, were carried out with different values for the initial solid-phase enantiomeric excess of D, with parameters as listed in Table 2 for simulations #8a
g. In all these simulations, the total initial mass was 0.08, and the final time was τ = 100, when the system has reached a steady state where the mass of both particle populations remains constant. The final liquid-phase enantiomeric excess of L is plotted against the final solid-phase enantiomeric excess of D. Since there is no racemization, the two populations evolve independently of each other. This means that each population evolves due to the effect of agglomeration, attrition and growth and dissolution. A steady state is reached when the size decrease due to attrition or dissolution balances the size increase due to agglomeration and growth, which depends on the particle 4617

dx.doi.org/10.1021/cg2008599 |Cryst. Growth Des. 2011, 11, 4611–4622

Crystal Growth & Design

Figure 7. Time needed to reach 90% solid phase enantiomeric excess as a function of the parameter ν, which is in turn a function of solubility. Lines are to guide the eye.

concentration (i.e., particle mass) and leads to a constant solute concentration in solution. An increase in particle mass leads to an increased agglomeration rate, and thus the steady state size is larger. This means that the critical size must also be larger, in order for dissolution to balance the size increase caused by agglomeration, thus leading to a lower concentration in solution. When considering two populations, a large difference in mass between the two populations means a large solid-phase enantiomeric excess. Thus, a high solid-phase enantiomeric excess of one enantiomer will, in the absence of racemization, lead to a liquidphase enantiomeric excess of the opposite enantiomer. Solubility Effect. It has been observed that an increase in solubility leads to a faster deracemization process.3 Experimentally, it is not possible to decouple the change in solubility (obtained by using different solvents) from the change in racemization rate also observed in different solvents; of course this is not a problem in the simulations where model parameters can be varied independently. Figure 7 shows results of simulations carried out with different values of the parameter ν = kvF/c∞ (simulations #9a
d). At lower values of ν (ν e 10
1), that is, high solubility, both populations dissolve. As observed in the figure, the simulations predict that the deracemization process becomes slower for increasing values of ν (decreasing solubilities), in accordance with the experimental observations. Sensitivity Analysis. For the sensitivity analysis, the kinetic parameters for racemization, agglomeration, and attrition were varied in a large range, in order to investigate the region in which chiral resolution is possible, and the effect of the different mechanisms on the process. The relative importance of growth and dissolution was varied by adjusting the capillary length α. The maximum solid phase enantiomeric excess reached is used as a measure of the extent of deracemization. In all simulations where complete chiral purity is obtained, the final enantiomeric excess exceeds 90% (100% is not always reached due to numerical difficulties). In other cases, the solid phase enantiomeric excess reaches a maximum and then starts decreasing again. Racemization Rate. The only parameter affecting the racemization rate is the kinetic rate constant, Krg. In simulations #10a
i, the value of Krg was varied over a wide range of values, spanning from 10
4 to 104. As shown in Figure 8a, even slow racemization leads to complete deracemization. Figure 8b shows that the time needed

ARTICLE

Figure 8. (a) Maximum solid phase enantiomeric excess reached and (b) time needed to reach a solid phase enantiomeric excess of 90% as a function of the racemization rate parameter Krg. Lines are to guide the eye.

Figure 9. (a) Maximum solid phase enantiomeric excess reached and (b) time needed to reach a solid phase enantiomeric excess of 90% as a function of the dimensionless capillary length α/x0. Lines are to guide the eye. Stars: lower initial mass; circles: higher initial mass.

to obtain chiral purity decreases strongly with an increasing racemization rate, as already observed in experiments.3,4 With an increase in the racemization rate, the flux of material from one enantiomeric form to the other becomes faster, which leads to faster symmetry breaking. This observation is valid as long as the racemization is rate-limiting. Above a threshold, which is located at Krg ≈ 1 in Figure 8b, an increase in the racemization rate does not reduce the deracemization time any further, as other mechanisms become controlling. The position of this threshold and the corresponding characteristic time may well depend on the other parameters, such as the agglomeration or breakage rates. As shown above, deracemization is not possible when there is no racemization in solution. Capillary Length. Growth and dissolution leads to very complex dynamics, which are affected by all the other mechanisms in the process: agglomeration and attrition have effects on the particle size distribution, while racemization changes the supersaturation in solution. Furthermore, the degree of growth and dissolution is affected by the capillary length, α. The smaller the capillary length, the smaller the critical size is at a given supersaturation. In the extreme case of α approaching zero, the critical size is virtually zero as well and there is no size-dependence of solubility. In the other extreme of very large values of α, the critical size is large even at high supersaturations. It is intuitive that for low values of α, dissolution does not occur to a 4618

dx.doi.org/10.1021/cg2008599 |Cryst. Growth Des. 2011, 11, 4611–4622

Crystal Growth & Design

Figure 10. (a) Maximum solid phase enantiomeric excess reached and (b) time needed to reach a solid phase enantiomeric excess of 90% as a function of the second parameter in the agglomeration rate expression. Lines are to guide the eye. Stars: eq 24; circles: eq 25

measurable extent, while an increase in α leads to an enhancement of growth and dissolution. In this section, we vary the value of α in order to investigate its influence on the deracemization rate. Figure 9a shows the maximum enantiomeric excess reached with different values of the capillary length α, plotted against the dimensionless capillary length α/x0. In Figure 9 the results of simulations #13a
i, with μ03,D = 0.05 and μ03,L = 0.03, are plotted as stars, and simulations #13j
n, with μ03,D = 50 and μ03,L = 30, are plotted as circles. In the simulations with the smaller initial mass, all the solid phase dissolves when α/x0 g 10
2. For higher initial mass, complete dissolution is observed when α/x0 g 102. The solid phase enantiomeric excess in these simulations reaches a maximum shortly before all the solid dissolves. The times required to reach 90% enantiomeric excess in the solid phase are shown in Figure 9b. As expected, the deracemization process is faster when the capillary length is larger. For very small values of α/x0, that is, very small capillary lengths, dissolution hardly takes place, and thus the mass of both populations remains constant. When the initial particle concentration is increased, agglomeration becomes more important due to the second-order dependence of the agglomeration rate on the particle concentration. As it will be shown later, faster agglomeration leads to a slower deracemization process; thus, a higher initial particle concentration will lead to a longer time until complete deracemization is reached, as shown in Figure 9b, where the circles represent simulations #13j
n with a larger initial mass than simulations #13a
i, shown as stars in the figure. Agglomeration. The two agglomeration kernels used in this work, eqs 24 and 25, have two parameters each. In both cases, the first parameter, rj,1, influences the shape of the agglomeration kernel, whereas the second changes the absolute value of the agglomeration rate. We have chosen to vary the second parameter, rj,2, in order to investigate the influence of the agglomeration rate on the deracemization rate. The parameter r1,2 in eq 24 was varied between 10
10 and
1 10 (simulations #14a
l), whereas the parameter r2,2 in eq 25 was varied between 5  10
15 and 10
4 (simulations #15a
k). As it can be seen in Figure 10b, the agglomeration rate has an influence on the time required to reach a solid phase enantiomeric excess of 90%. The trend is, however, not monotonic. An agglomeration rate which is either too low or too high leads to an

ARTICLE

Figure 11. (a) Maximum solid phase enantiomeric excess reached and (b) time needed to reach a solid phase enantiomeric excess of 90% as a function of the attrition rate parameter Kbg. Lines are to guide the eye.

Figure 12. Third moment of distributions and solid phase enantiomeric excess of D as a function of time for simulation #11j, Kbg = 104.

increased deracemization time, and, as shown in Figure 10a, a low value can also hinder the deracemization process completely. This behavior can be explained by considering the effect of agglomeration on the distribution of particle sizes. The deracemization relies on exploiting differences in the two particle populations. Mechanisms which are first-order in the particle concentration are affected by differences in the particle sizes only, whereas second-order processes are affected by differences in particle sizes and in particle concentration. Agglomeration is the only mechanism considered here which is second-order; thus, it is the only way that a difference in mass between the two populations can be exploited to reach chiral resolution. If agglomeration is slow, it takes longer before a mass difference between the two populations leads to complete deracemization. The flux of material from one chiral form to the other relies on dissolution and on the subsequent racemization. If agglomeration is too fast, the small crystals will rather agglomerate than dissolve, thus impeding the conversion between enantiomers and thus the whole deracemization process. Attrition. The attrition rate parameter Kbg was varied from 10
5 to 105 (simulations #11a
k), that is, covering a very broad range of attrition rates. These rates are not necessarily achievable experimentally, and very high attrition rates would probably lead to other problems not considered in this model, such as phase transitions.9 As readily observed in Figure 11a, chiral purity is obtained for all but very high attrition rates. The attrition rate has a strong effect on the time taken to reach chiral purity, which decreases with an increasing attrition rate parameter, as shown in Figure 11b. When the attrition rate is too high, very many small 4619

dx.doi.org/10.1021/cg2008599 |Cryst. Growth Des. 2011, 11, 4611–4622

Crystal Growth & Design

ARTICLE

words, there is no loss of generality by limiting the analysis to the constitutive equations presented in this work.

Figure 13. Time evolution of the solid phase enantiomeric excess of D for three simulations with different values of Kbg.

particles appear. This leads to complete dissolution of all particles, as seen in Figure 12a, which shows the evolution of the third moment of both populations in simulation #11j. The enantiomeric excess increases initially but reaches a maximum and starts decreasing again toward the end of the process, when nearly all L crystals have dissolved and the mass of D crystals is still decreasing, as shown in Figure 12b. The behavior is similar when Kbg = 105 (simulation #11k). For low values of Kbg (simulations #11a
d, 10
5 e Kbg e 10
2), the solid phase enantiomeric excess no longer follows an exponential increase. Rather, the initial increase is fast, but then it slows down. When there is no attrition (Kbg = 0, simulation #12), the curve is again different, with the initial fast increase being even more pronounced (and initially even faster than for Kbg = 10
4). The rate of change of the enantiomeric excess decreases steadily up to ee = 100%. A comparison of three cases is shown in Figure 13. Experiments in a system without glass beads showed enantiomeric excess profiles similar to those obtained here for low values of Kbg. 4,8,43 It is worth noting that complete deracemization is possible even when there is no attrition.

’ DISCUSSION OF ASSUMPTIONS The constitutive equations eqs 5
10 used to describe the kinetics of the elementary mechanisms involved in Viedma ripening represent a choice among a relatively large number of possible alternatives. Other growth rate expressions than eq 5 could be considered, for example, diffusion-limited growth, or birth and spread growth. We have used two different agglomeration kernels (eqs 6 and 7), but more could, in principle, be considered. An unlimited number of daughter distributions for breakage instead of eq 9 would be possible. In some cases we have explored these alternatives and discussed the results above. In other cases, a discussion is not included for the sake of brevity. All of our calculations, as well as our physical intuition, show that the phenomena observed depend on a few key features of the physical mechanisms and of the related models: (1) presence of a racemization reaction; (2) size-dependence of solubility; (3) the occurrence of agglomeration of crystals of the same chirality at a rate which is secondorder in the particle concentration; (4) the breakage of crystals at a rate which is first-order in the particle concentration. In other

’ CONCLUSIONS Viedma ripening is an important, new process which is interesting in its own right. In addition, it has been considered for potential practical applications3,9,16 and for its possible role in the evolution of homochirality in biologically relevant molecules.1,2,7,16,46,47 This work presents a comprehensive model based on population balances which is able to describe the behavior observed experimentally during deracemization by accounting for all the relevant mechanisms, namely, racemization, agglomeration, and attrition as well as growth and dissolution. The model combines PBEs which have been used successfully for each mechanism separately in the past and avoids many simplifications and assumptions made in previous models of Viedma ripening. One of the advantages of a model such as the one presented here is that it allows the decoupling of the different mechanisms or the selective suppression of one or more of them, something that is not possible experimentally but is very useful in order to understand the role played by each individual mechanism separately in such a complex process. The analysis carried out here shows that all the mechanisms considered have an effect and that their interplay is complex. Furthermore, it shows that racemization, growth, and dissolution due to a size-dependent solubility and agglomeration are necessary mechanisms; that is, without any of them Viedma ripening does not occur. Each of these mechanisms has an effect on the deracemization rate. On the contrary, attrition is not necessary but does indeed increase the deracemization rate; in other words, Viedma ripening occurs also when breakage of crystals is suppressed but in that case it is very slow. The model presented in this work fills crucial gaps in understanding Viedma ripening and contributes to making its exploitation possible. ’ APPENDIX A: DERIVATION OF BOUNDARY CONDITIONS AND MOMENTS EXPRESSIONS Since the PBE, eq 11, includes size-dependent growth rates, the standard boundary conditions applied to PBEs are not valid in this case. Furthermore, the correct expressions have to be used for the time-derivatives of the distribution moments. To find the boundary conditions and moment derivatives for the system presented here, eq 11 is multiplied by xj and integrated over x from 0 to ∞: Z ∞ Z ∞ ∂ni ∂ðGi ðxÞni Þ dx xj dx ¼
xj ∂x ∂t 0 0 Z ∞ Z ∞ þ xj ni ðη, tÞbðηÞgðx, ηÞ dη dx 0 Z ∞ x
xj ni ðx, tÞbðxÞ dx 0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z ∞ 2 þ jZ x x Aðξ, 3 x3
ξ3 Þ þ ni ðξ, tÞ 2 0 ðx3
ξ3 Þ2=3 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3  ni ð x3
ξ3 , tÞ dξ dx Z ∞ Z ∞
xj ni ðx, tÞ Aðξ, xÞni ðξ, tÞ dξ dx ð30Þ 0

4620

0 dx.doi.org/10.1021/cg2008599 |Cryst. Growth Des. 2011, 11, 4611–4622

Crystal Growth & Design

ARTICLE

In the term on the left-hand side, the order of the integration and the differentiation can be switched, giving the time derivative of the jth moment. The growth term, that is, the first term on the right-hand side, can be integrated by parts (j ¼ 6 0): Z ∞ ∂ðGi ðxÞni Þ dx ¼ ½xj Gi ðxÞni ðx, tÞ∞ xj 0 ∂x 0 Z ∞
j xj
1 Gi ðxÞni ðx, tÞ dx ð31Þ 0

The attrition terms, that is, the second and third terms on the righthand side of eq 30, can be simplified to Z ∞ Z η Z ∞ ni ðη, tÞbðηÞ xj gðx, ηÞ dx dη
xj bðxÞni ðx, tÞ dx 0

0

becomes

Z ∞ dμi, 3 ¼ ½x3 Gi ðxÞni ðx, tÞx ¼ 0 þ 3 x2 Gi ðxÞni ðx, tÞ dx dt 0 ð38Þ

The first term on the right-hand side of eq 38 is also zero, since ) x3 f 0 asx f 0 ð39Þ 0 > Gi ð0Þni ð0, tÞ >
∞ leading to Z ∞ dμi, 3 ¼3 x2 Gi ðxÞni ðx, tÞ dx dt 0

0

ð32Þ Now, we consider the zeroth moment, j = 0. The agglomeration terms, the last two terms in equation eq 30, can be written as Z Z 1 ∞ ∞
Aðx, ξÞni ðx, tÞni ðξ, tÞ dx dξ ð33Þ 2 0 0 Using eqs 31, 32, and 33, eq 30 becomes Z ∞ dμi, 0 ¼ ½Gi ðxÞni ðx, tÞx ¼ 0 þ ðNF
1Þ bðxÞni ðx, tÞ dx dt 0 Z Z 1 ∞ ∞
Aðx, ξÞni ðx, tÞni ðξ, tÞ dx dξ ð34Þ 2 0 0 since there are no infinitely large particles and the integral over the daughter distribution must give the number of fragments:38,42,48 ni ðx f ∞, tÞ ¼ 0 Z

which can be calculated when Gi(x) and ni(x,t) are known.

’ APPENDIX B: NUMERICAL SOLUTION Equation 23 is solved using the method of pivots.38 The growth term is discretized using a grid with M points. The distributions are discretized assuming all particles are concentrated at the gridpoints. M

∑ Fi, m δðy
ym Þ m¼1

fi ðy, tÞ ¼

ð40Þ

The zeroth and third moments of the distributions are conserved in the attrition and agglomeration terms.38 The discretized PBE becomes !   dFi, m Gi ðyb ÞFi, b
Gi ðya ÞFi, a ¼
dτ yb
ya ! M

Fi, k Qm, k yβk
Fi, m yβm ∑ k¼m

þ Kbg

ð35Þ


 1 þ 1
δjk 2 ym
1 e ðyj þ yk Þ < ym
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3
1Þ y3j þ y3k Ad ðyj , yk ÞFi, j Fi, k  γðm m jgk

∞

gðx, ηÞ dx ¼ NF

∑

ð36Þ

0

If we assume no growth, the only change in number of particles is due to attrition and agglomeration. The second term in eq 34 must be positive and finite, since attrition can only lead to an increase in number of particles but cannot lead to an infinite increase. The term given in eq 33 must be negative and finite. If we assume no attrition (b(x) = 0) and no agglomeration (A(x,ξ) = 0), the change in total number of particles is given by the growth term Gi(x)ni(x,t). Since, for this problem, Gi(x) must be such that small particles dissolve, this term must be negative. Also, it must be finite; otherwise, all particles would disappear instantaneously. Thus, 0 > Gi(0)ni(0,t) >
∞. Since Gi(0) f
∞, then ni(0,t) must be zero. Now we consider the third moment, j = 3. Since conservation of mass on attrition requires that38,42,48 Z ∞ x3 gðx, ηÞ dx ¼ η3 ð37Þ


 1 1
δjk 2 ym e ðyj þ yk Þ < ymþ1
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 y3j þ y3k Ad ðyj , yk ÞFi, j Fi, k  γðmÞ m

∑


Fi, m Z Qm, k ¼

ym

ym
1

Z

þ

0

the two attrition terms in eq 32 cancel. One can also show that the agglomeration terms in eq 33 cancel, as expected since agglomeration cannot lead to a change in mass. Thus, eq 30

jgk

þ

M

∑ Ad ðym, yj ÞFi, j j¼1

ð41Þ


1Þ γðm ðuÞgðu, yk Þ du m ymþ1

ym

γðmÞ m ðuÞgðu, yk Þ du

ð42Þ

where a and b are set depending on which scheme is used at the current point and δjk is the Kronecker delta. A central difference scheme is used at all points, except at the edges of the grid. At the 4621

dx.doi.org/10.1021/cg2008599 |Cryst. Growth Des. 2011, 11, 4611–4622

Crystal Growth & Design

ARTICLE

lower edge, a forward difference is used, and at the last point, a second-order backward difference is used, that is: 1 < m < M : a ¼ m
1, m ¼ 1: a ¼ 1, m¼M: a ¼ m
2,

b ¼ m þ 1, b ¼ 2, b ¼ M:

The quantities Qm,k which define a matrix Q are time-independent and can be calculated once, then stored and reused in later iterations. A derivation of the expression for Q and formulas for γ are given by Ramkrishna.38 For the sake of numerical efficiency, wherever possible the terms in the sums of eq 41 are calculated once only and stored in a look-up table. The time derivative of the third moment (eq 27) in discretized form is given by M dji, 3  ¼3 y2m Gi ðym ÞFi, m dτ m¼1

∑

ð43Þ

Equation 41 is solved for both populations and all pivots. Equation 43 is solved for both populations along with eqs 28 and 29; hence 2M + 4 ODEs have to be solved. The scheme is solved using the ODE15s solver which is built-in in Matlab.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT This work is a contribution to the research project INTENANT and is therefore partially funded by the European Commission within the 7th Framework Programme (Project No. 214 129). ’ REFERENCES (1) Viedma, C. Phys. Rev. Lett. 2005, 94, 065504. (2) Noorduin, W. L.; Izumi, T.; Millemaggi, A.; Leeman, M.; Meekes, H.; Van Enckevort, W. J. P.; Kellogg, R. M.; Kaptein, B.; Vlieg, E.; Blackmond, D. G. J. Am. Chem. Soc. 2008, 130, 1158–1159. (3) Noorduin, W. L.; Meekes, H.; van Enckevort, W. J. P.; Millemaggi, A.; Leeman, M.; Kaptein, B.; Kellogg, R. M.; Vlieg, E. Angew. Chem., Int. Ed. 2008, 47, 6445–6447. (4) Tsogoeva, S. B.; Wei, S.; Freund, M.; Mauksch, M. Angew. Chem. 2009, 121, 598–602. (5) Wei, S.; Mauksch, M.; Tsogoeva, S. B. Chem.—Eur. J. 2009, 15, 10255–10262. (6) Noorduin, W. L.; van Enckevort, W. J. P.; Meekes, H.; Kaptein, B.; Kellogg, R. M.; Tully, J. C.; McBride, J. M.; Vlieg, E. Angew. Chem., Int. Ed. 2010, 49, 8435–8438. (7) Cheung, P. S. M.; Gagnon, J.; Surprenant, J.; Tao, Y.; Xu, H. W.; Cuccia, L. A. Chem. Commun. 2008, 987–989. (8) Viedma, C.; Ortiz, J. E.; de Torres, T.; Izumi, T.; Blackmond, D. G. J. Am. Chem. Soc. 2008, 130, 15274
15275. (9) Noorduin, W. L.; van der Asdonk, P.; Bode, A. A. C.; Meekes, H.; van Enckevort, W. J. P.; Vlieg, E.; Kaptein, B.; van der Meijden, M. W.; Kellogg, R. M.; Deroover, G. Org. Process Res. Dev. 2010, 14, 908–911. (10) Noorduin, W. L.; van der Asdonk, P.; Meekes, H.; van Enckevort, W. J. P.; Kaptein, B.; Leeman, M.; Kellogg, R. M.; Vlieg, E. Angew. Chem. 2009, 121, 3328–3330. (11) Kaptein, B.; Noorduin, W. L.; Meekes, H.; van Enckevort, W. J. P.; Kellogg, R. M.; Vlieg, E. Angew. Chem., Int. Ed. 2008, 47, 7226–7229. (12) Kondepudi, D. K.; Kaufman, R. J.; Singh, N. Science 1990, 250, 975–976.

(13) Levilain, G.; Rougeot, C.; Guillen, F.; Plaquevent, J.-C.; Coquerel, G. Tetrahedron: Asymmetry 2009, 20, 2769–2771. (14) Leeman, M.; Noorduin, W. L.; Millemaggi, A.; Vlieg, E.; Meekes, H.; van Enckevort, W. J. P.; Kaptein, B.; Kellogg, R. M. CrystEngComm 2010, 12, 2051–2053. (15) Skrdla, P. J. Cryst. Growth Des. 2011, 11, 1957–1965. (16) McBride, J. M.; Tully, J. C. Nature 2008, 452, 161–162. (17) Noorduin, W. L.; Meekes, H.; Bode, A. A. C.; van Enckevort, W. J. P.; Kaptein, B.; Kellogg, R. M.; Vlieg, E. Cryst. Growth Des. 2008, 8, 1675–1681. (18) Mersmann, A. Crystallization Technology Handbook, 2nd ed.; Marcel Dekker: New York, 2001.  (19) Stahl, M.; Åslund, B.; Rasmuson, A. C. Ind. Eng. Chem. Res. 2004, 43, 6694–6702. (20) Uwaha, M. J. Phys. Soc. Jpn. 2004, 73, 2601–2603. (21) Uwaha, M. J. Phys. Soc. Jpn. 2008, 77, 083802. (22) Saito, Y.; Hyuga, H. J. Phys. Soc. Jpn. 2005, 74, 535–537. (23) Wattis, J. Origins Life Evol. Biospheres 2011, 41, 133–173. (24) Cartwright, J. H. E.; Piro, O.; Tuval, I. Phys. Rev. Lett. 2007, 98, 165501. (25) Uwaha, M.; Katsuno, H. J. Phys. Soc. Jpn. 2009, 78, 023601. (26) Uwaha, M. J. Cryst. Growth 2011, 318, 89
92. The 16th International Conference on Crystal Growth (ICCG16)/The 14th International Conference on Vapor Growth and Epitaxy (ICVGE14). (27) Katsuno, H.; Uwaha, M. J. Cryst. Growth 2009, 311, 4265–4269. (28) Saito, Y.; Hyuga, H. J. Phys. Soc. Jpn. 2008, 77, 113001. (29) Saito, Y.; Hyuga, H. J. Phys. Soc. Jpn. 2009, 78, 104001. (30) Saito, Y.; Hyuga, H. J. Phys. Soc. Jpn. 2010, 79, 083002. (31) Saito, Y.; Hyuga, H. J. Cryst. Growth 2011, 318, 93
98, The 16th International Conference on Crystal Growth (ICCG16)/The 14th International Conference on Vapor Growth and Epitaxy (ICVGE14). (32) Hatch, H. W.; Stillinger, F. H.; Debenedetti, P. G. J. Chem. Phys. 2010, 133, 224502. (33) Noorduin, W. L.; Kaptein, B.; Meekes, H.; van Enckevort, W. J. P.; Kellogg, R. M.; Vlieg, E. Angew. Chem., Int. Ed. 2009, 48, 4581–4583. (34) Lifshitz, I. M.; Slyozov, V. V. J. Phys. Chem. Solids 1961, 19, 35–50. (35) Wagner, C. Z. Elektrochem. 1961, 65, 581–591. (36) Tavare, N. S. AIChE J. 1987, 33, 152–156. (37) Valentas, K. J.; Amundson, N. R. Ind. Eng. Chem. Fundamentals 1966, 5, 533–542. (38) Ramkrishna, D. Population Balances: Theory and Applications to Particulate Systems in Engineering; Academic Press: San Diego, CA, 2000. (39) Lindenberg, C.; Sch€ oll, J.; Vicum, L.; Mazzotti, M.; Brozio, J. Cryst. Growth Des. 2008, 8, 224–237. (40) Valentas, K. J.; Bilous, O.; Amundson, N. R. Ind. Eng. Chem. Fundamentals 1966, 5, 271–279. (41) Randolph, A. D. Ind. Eng. Chem. Fundamentals 1969, 8, 58–63. (42) Ziff, R. M. J. Physics A 1991, 24, 2821–2828. (43) Izumi, T.; Klussmann, M.; Viedma, C.; Blackmond, D. G. How did Amino Acids Take a Left-Hand Turn? Poster presentation at 2010 Annual Conference of the British Association of Crystal Growth, Manchester, September, 5
7, 2010. (44) Noorduin, W. L.; Vlieg, E.; Kellogg, R. M.; Kaptein, B. Angew. Chem., Int. Ed. 2009, 48, 9600–9606. (45) Siegel, J. S. Chirality 1998, 10, 24–27. (46) Viedma, C. Cryst. Growth Des. 2007, 7, 553–556. (47) Viedma, C. Astrobiology 2007, 7, 312–319. (48) Kostoglou, M.; Dovas, S.; Karabelas, A. J. Chem. Eng. Sci. 1997, 52, 1285–1299.

4622

dx.doi.org/10.1021/cg2008599 |Cryst. Growth Des. 2011, 11, 4611–4622