A Thermodynamic Approach for the Prediction of Oiling Out Boundaries from Solubility Data

: Many pharmaceutical molecules, fine chemicals, and proteins exhibit liquid–liquid phase separation (LLPS, also known as oiling out) during solution crystallization. LLPS is of significant concern in crystallization process development, as oiling out can compromise the effectiveness of a crystallization and can lead to operational problems. A comprehensive methodology that allows a process scientist/engineer to characterize the various phase boundaries relevant to oiling out is currently lacking. In this work, we present a modeling framework useful in predicting the binodal, spinodal, and gelation boundaries starting from the solubility data of a solute that is prone to oiling out. We collate the necessary theoretical concepts from the literature and describe a unified approach to model the phase equilibria of solute–solvent systems from first principles. The modeling effort is validated using experimental data reported in the literature for various solute– solvent systems. The predictive methods presented in this work can be easily implemented and help a process engineer establish the design space for a crystallization process that is affected by liquid– liquid phase separation.


Introduction
Solution crystallization is an important unit operation typically practiced in the pharmaceutical and fine chemical industry to produce solid compounds with desired purity and physical attributes [1,2].The design of an industrial crystallization process is generally concerned with the control of the nucleation and growth of solid particles.However, for many organic active pharmaceutical ingredients (APIs) and biomolecules that form molecular crystals (as opposed to inorganic salts which form ionic crystals), this task may not be straightforward because these systems can exhibit unusual phase separation phenomena during solution crystallization.For many APIs and proteins, often the result of a crystallization attempt may not be the formation of solid particles, but that of metastable dense liquid phases and gels.These intermediate phases typically occur due to liquidliquid phase separation (LLPS-which in the process engineering parlance is known as "oiling out") of the solute, a phenomenon that complicates the process development effort.The occurrence of these intermediate dense liquid and gel phases from oiling out is usually considered undesirable, as these phases can lead to uncontrollable nucleation, nucleation of metastable crystal forms, impurity entrapment, and other operational difficulties [3][4][5][6][7][8].However, circumstances may exist in which a crystallization mediated through these phases may be needed to achieve solid particles of special attributes, e.g., clusters and spherical particles that easily filter downstream [9][10][11].Thus, a comprehensive understanding of the phase behavior of the solute is essential for a rational design and trouble-free operation of a crystallization process.
While the solid-liquid equilibrium boundary (solubility curve) is easily determined experimentally for many solute-solvent systems, a complete characterization of the liquid-liquid and liquid-gel boundaries remains a daunting task that requires a contextual understanding of the phenomena and laborious experimentation.Two different oiling out boundaries are commonly encountered in experiments-the binodal and the spinodal (see Section 2 below).When a process scientist observes an oiling-out phenomenon in a (typically cooling) crystallization experiment, most likely they are in the vicinity of the binodal at that process condition.This is because penetrating the spinodal envelope requires rapid change in process conditions (quenches) that are typically not achievable in cooling crystallization that involves solutions of significant volume.Thus, a portion of the binodal curve may be accessible to experimenters using techniques, such as focused beam reflectance measurement (FBRM) [7,12].However, a complete experimental characterization of the binodal (including the dense liquid branch) is rather difficult.Hence a strong need exists among the crystal growth community to model and predict the liquid-liquid phase boundaries for the solute of interest under process conditions.
The aspect that distinguishes real systems that form various condensed phases from an ideal gas is the existence of intermolecular interactions.For this reason, the study of phase behavior of molecules starts with the description of intermolecular interactions.A large body of literature exists in the condensed matter physics and colloidal science fields that describes methods for the prediction of various phase boundaries for molecules with well-defined interaction potentials.Theoretical approaches developed to this end include perturbation methods that describe the liquid structure [13][14][15][16][17][18][19][20], integral equation approaches [21,22], and the recent density functional methods [23][24][25], which are considered more rigorous.In general, the theoretical predictions of these models are validated against the results obtained from computer simulations (also known as machine calculations) [26][27][28][29][30], the results of which are regarded as "experimental data" in this context.
Given that the various theories of phase behavior of molecules described above mainly exist in the computational chemistry and condensed matter physics domains, unless a crystallization process engineer is well-versed with this literature, they are at a disadvantage with respect to conceiving a suitable predictive approach.Also, while these theoretical methods have evolved over time to various degrees of sophistication, as a rule, the computations involved in many of these approaches are extremely complicated and are tedious to be executed by practicing scientists and engineers in the context of industrial process development.Hence, in general, these methods remain obscure and foreign to industrial practitioners of crystallization and process modelers.Moreover, all these models require the intermolecular interaction parameters (such as the range and strength of interaction) to be pre-defined.This information is not readily available for systems of industrial interest.To address these gaps, here, we present a methodology to predict the various phase boundaries of solutes using simple fluid models.We collate the necessary theoretical concepts from various sources in the literature, and propose a modeling approach that is simple enough to be implemented readily by process scientists.Our method derives the necessary molecular interaction parameters from experimental solubility data and uses these parameters to predict the liquid-liquid and liquid-gel boundaries.Apart from the availability of an experimentally determined solubility data set, the approach we discuss here requires only the knowledge of the molecular weight of a compound and the crystal density.
In Section 2, we first outline the theoretical underpinnings of the proposed modeling approach.In Section 3, we consider a few solute-solvent systems discussed in the literature for which experimental data on the solubility and binodal are available.In Section 4, we model the solubility data for these systems, and from that predict the liquid-liquid phase boundaries (binodal and spinodal).This effort establishes the effectiveness of the method.Section 5 that follows discusses the context and limitations of the proposed approach.Finally, we summarize the work in Section 6.

Theory and Model Development
When quenched rapidly, a system that contains a dissolved solute can become unstable with respect to the existing single phase composition, and may undergo a liquid-liquid phase separation through a spinodal decomposition mechanism before it can nucleate crystals [31].The locus of the temperature and composition points that define the unstable region of the phase diagram is known as the spinodal curve.Upon spinodal decomposition, the system typically splits into two liquid phases, with one phase containing a high concentration of the solute and the other being lean in solute.These two liquid phases can co-exist in a local equilibrium.The locus of the temperature and compositions that define this local equilibrium is known as a coexistence or binodal curve.The binodal region completely envelopes the spinodal region.For molecules experiencing short-ranged interactions, the region of the phase diagram given by the liquid-liquid coexistence curve (binodal) is metastable with respect to the solid-liquid coexistence curve [22,28,32].In such a case, the two liquid phases resulting from LLPS are expected to produce crystals eventually.
The basis for the phase equilibrium calculations in which we are interested is the thermodynamic condition that at equilibrium the chemical potentials and pressures of the co-existing phases are equal.Thus, one needs to obtain expressions for the chemical potential of the solid and liquid phases to be able to determine the co-existence boundaries.To that end, we start by considering the solute in the crystallizing solution as a simple fluid and use a square-well potential to describe the intermolecular interactions.Several studies in the literature have established the effectiveness of the square-well interaction model in capturing the solution thermodynamics for systems that exhibit liquid-liquid phase separation [33][34][35][36].The square-well model defines the pair interaction energy, U(r), as a function of the distance, r, between molecules as: in which σ corresponds to the molecular diameter of the species in consideration, and ε and λ denote the strength and range of the interaction, respectively.Below, we first discuss the liquid (fluid) phase chemical potential and pressure for a square-well fluid.The chemical potential for a square-well solid is discussed in the subsection that follows.

Fluid Phase Thermodynamics
For describing the solution thermodynamics of a square-well fluid, we follow the approach by Gil-Villegas et al. [37] and Patel et al. [38] in which the molecular liquid phases were described using the statistical associating fluid theory for chain molecules with attractive potentials of variable range (SAFT-VR).Detailed derivation of the equation of state for SAFT-VR square well fluids is available in the original works [37,38].Here, we limit our discussion to the essential elements needed to understand the phase diagram calculations.
By considering a hard sphere fluid [39] as a reference fluid, the expression for the Helmholtz free energy, A, of N particles interacting via the square-well potential can be written using a second order perturbation theory expansion [18] as: in which k is the Boltzmann's constant and T is the absolute temperature of the system.In Equation ( 2), the first term on the right-hand side (RHS) of the equality sign (A ideal ) is the ideal gas Helmholtz free energy, the second term (A HS ) is the excess free energy of the hard-sphere (reference) fluid.Thus, the sum of the first two terms on the RHS gives the free energy of the reference fluid.The remaining two terms, A₁ and A₂, express the first order and second order perturbation corrections to the reference fluid to yield the Helmholtz free energy of the fluid of interest.

Free Energy of the Reference Fluid
The ideal contribution to the free energy, A ideal , is given by: where ρ = N/V is the molecular number density of the fluid, with V being the volume of the system and Λ being the thermal De Broglie wavelength.To describe the hard sphere (reference fluid) excess Helmholtz free energy (A HS ), the Carnahan-Staling expression is used [40,41]: ρ σ is the hard sphere packing (volume) fraction.

First Perturbation Term
The first perturbation term, A₁, corresponds to the mean-attractive energy of the fluid, and uses the radial distribution function (RDF) of the reference hard sphere fluid, g HS (x;η), to approximate the structure of the square-well fluid as [18]: in which x = (r/σ).The mean value theorem [42] can be used to simplify the right-hand side of the equation above to give [37]: in which ξ is the value of x, which satisfies the equality: As argued in [37], it is possible to simplify this expression further by assuming that the value of the hard-sphere RDF at a given ξ and η can be mapped on the contact value of the hard-sphere RDF at an effective packing fraction, eff , η defined by: The Carnahan-Starling expression for the contact value of the hard-sphere RDF [43] is then used to obtain the RHS of Equation ( 8) as: The resulting expression for the first-order perturbative contribution to the free energy is given by: Equation (11) , where the coefficients, i b , were given by fourth-order polynomials in (1/λ) as:

Second Perturbation Term
The second order correction term, A₂, in Equation ( 2) describes changes in the energy due to compression of the fluid arising from the attractive interaction.This term is approximated using the local compressibility approximation [14], which yields for the square-well potential [37]: where HS K is the isothermal compressibility of the hard-sphere reference fluid, estimated by the Percus-Yevick expression as [18]:

Chemical Potential and Pressure
Calculating phase equilibria requires expressions for the chemical potential, μL, and pressure, pL, of the liquid phases.These quantities can be obtained from the Helmholtz free energy using the standard thermodynamic relations: , The above relations can be written in a more convenient form using the product rule of differentiation and the definition of Helmholtz free energy as:

Solid Phase Thermodynamics
The required calculations are greatly simplified when one assumes the solid phase to be incompressible, in which case an expression for the pressure of the molecular solid is not needed.For the excess chemical potential of a square-well solid, μS, we use the simple expression due to Asherie et al. [33] and write: in which nS is the number of molecules whose centers lie within the attractive well of a given molecule.This number is determined by the structure of the crystalline lattice.As a simplification, we follow the approach of Asherie et al., and assume the solid to be of face-centered cubic (fcc) structure and hence consider nS = 12, i.e., the attractions in the solid are dominated by the molecules in the first coordination shell.This model for the chemical potential of a square-well solid has been found to represent the solid-liquid equilibria of several small molecules and colloidal solutions reasonably well [11,33,44,45].

Gelation Boundary
The experimental observations of gel phase formation and the resulting coating on the reactor walls during oiling out also motivate us to present a model to estimate this gelation boundary for a given solute.For square-well fluids, on the basis of mode-coupling theory arguments, Bergenholtz et al. [46] have derived a simple expression for the volume fraction, ηg, at which gelation occurs at a given temperature, T. This model relates ηg to the intermolecular interactions through: 1 e 1 1.42.
Thus, if one has the knowledge of ε and λ, one can determine the gelation region on a phase diagram by solving Equation ( 18) for ηg at various temperatures.

Model Implementation
In obtaining the free energy, A, of a molecular fluid, the calculation of A₁ through Equations ( 5) -( 8) is the most challenging and tedious part.Thanks to the approach conceived by Gil-Villegas [37], this difficulty is greatly circumvented.One can simply calculate A₁ through Equations ( 10) -( 12), and need not be concerned with obtaining the radial distribution function, g HS , of the reference hardsphere fluid and solving Equation ( 5) numerically.These simplifications make the following approach attractive and readily implementable with minimal computational time.

Obtaining Molecular Interaction Parameters
To calculate the phase boundaries for the solute-solvent system of interest, the molecular interaction (square-well) parameters, ε and λ, need to be known.When experimental data on solidliquid equilibrium (i.e., solubility data) are available, these parameters can be deduced as follows.First, one converts the solubility concentrations (typically available in the units of g/L or wt%) to the number densities, ρ, and from ρ to the volume fractions, η (=(π/6)ρσ³).Here, σ represents the equivalent spherical diameter of the solute molecule (see Section 3.6).One then proceeds to obtain the square-well parameters, ε and λ, through the regression of this η vs. T data set to the relation, μL = μS.Note that the other condition of phase equilibrium, the equality of pressures, is neglected because the solid phase is considered to be incompressible.
The strength of interaction, ε, itself is expected to be a function of temperature, which may be expressed through a simple linear relationship [45] as ε/k = α₀+ α₁ T. One may choose other forms for the temperature dependence of ε/k.However, one needs to ensure that the total number of adjustable parameters in the model does not increase unreasonably.With this linear expression, we have three parameters (α₀, α₁, and λ) that need to be fit to describe the solubility data.For this, first one provides initial guesses for α₀, α₁, and λ and calculates μS through Equation ( 17) for each T in the experimental data set.The equation of μL = μS needs to be solved iteratively by providing an appropriate starting value for ηm, the volume fraction, η, that represents the predicted solubility (solution of the equation) at the chosen T. This initial value of ηm enables one to compute the Helmholtz free energy of the reference fluid by evaluating Equations ( 3) and (4) for A ideal and A HS , respectively.The De Broglie wavelength, Λ, in Equation ( 3) is neglected, as this term effectively cancels out as an ideal contribution on both sides of the equality of μL = μS at constant T. One then has a choice of evaluating ηeff using either Equation (11) or (12).From ηeff, A₁ is obtained using Equation (10), and A₂ is obtained by numerically differentiating A₁ and using that derivative in Equation (13).Once the total Helmholtz free energy, A, is obtained for the guessed ηm, this expression is numerically differentiated to evaluate μL (Equation ( 16)).Thus, the relation of μL = μS is solved to yield the value of ηm for each temperature in the experimental data set.The best-fit model parameters, α₀, α₁, and λ, are those that minimize the objective function: i.e., that minimize the sum of squared errors between the predicted ηm and experimental η for the nexp data points considered.

Practical Considerations
A few practical aspects need to be noted toward implementing the calculation procedure discussed above.The relation of μL = μS may have up to three solutions (i.e, one, two, or three solutions) for ηm due to the cubic nature of the μL vs. η curve.The volume fraction, η, is obviously constrained between 0 and 1, but volume fractions greater than 0.5 are unphysical (for the liquid phase).Hence, it is desirable to constrain the numerical solution of μL = μS to 0 ≤ ηm ≤ 0.5.Finding the ηm that corresponds to the minimum of all the possible solutions in this range (which is a guaranteed solution) ensures that the solution obtained is physically meaningful and the iterative procedure for parameter estimation converges.Also, one needs to be mindful of the range of λ over which the chosen model for ηeff (Equation ( 11) or ( 12)) is applicable and constrain the search space for λ appropriately during the regression.The strength of interaction, ε, between the molecules is expected to be in a range of 1 to 5 kT.Hence, one may have to limit the search space for α₀ and α₁ accordingly during parameter estimation, so that the values of the resulting ε/k (= α₀ + α₁ T) fall in the range of 1 to 5 T, where T is the absolute temperature.The bounds on the possible parameter space and the resulting constrained error minimization, while help keep the result of the calculation physically relevant, make it hard to estimate the uncertainties in the parameters.

Calculation of Other Phase Boundaries
Once the square-well interaction parameters (α₀, α₁, and λ, with ε/k = α₀ + α₁ T) are known, the co-existence curve (binodal) may be obtained by solving simultaneously I η are then passed on as initial guesses for the next temperature, which will be slightly less than the previous one (about 0.5 °C).This procedure results in quick convergence of the iterative numerical computation and will not result in trivial solutions.The spinodal curve can also be calculated using this procedure.In this case, we need for the two solutions of η at a given temperature.With known square-well interaction parameters, the calculation of the gel boundary from Equation ( 18) is straightforward for a given temperature, T.

Experimental Section -Methods
To demonstrate the effectiveness of the above approach toward modeling different phase boundaries, we considered experimental data available in the literature for various solute-solvent systems.We found that experimental studies that report a detailed characterization of the binodal and/or spinodal for solutes are rare.Among the few studies available in which sufficient data on the solubility and binodal are reported together, we found five systems for which the necessary information on the molecular weight of the compound is available.The data on three of these systems are used here to show the usefulness of the proposed modeling approach, and the other two systems are discussed in the context of the limitations of the approach.Below, we summarize the experimental methods used by the authors.For detailed discussion on the experimental aspects of each study, the reader is referred to the original publications.

Pyraclostrobin -Isopropanol/Cyclohexane
First, we consider experimental data on the solubility and oiling-out (binodal) curves reported by Li et al. [4] on pyraclostrobin (form II crystals) crystallizing in 10% (w/w) methanol -90% (w/w) cyclohexane mixture.Li et al. determined the solubility of these crystals in the temperature range of 5 to 35 °C using a laser obscuration technique that monitored the solutions for the dissolution of solids during gradual addition of solute at a constant temperature.The oiling-out data points at various concentrations were determined using the FBRM/PVM instruments.These authors obtained the oiling-out data points with multiple cooling rates.Here, we consider the LLPS data set generated with the slowest cooling rate reported (0.1 K/min), as this data set is expected to be close to the true equilibrium binodal exhibited by the system.

Compound Z -Methanol/Water
The next system we consider was recently reported by Bhamidi et al. [11].This compound (compound Z) was crystallized in 2% (w/w) methanol -98% (w/w) water.The solubility data were determined in the temperature range of 5 to 65 °C using the Crystal 16 apparatus at a 0.5 °C/min cooling rate.The oiling-out data were also determined using this apparatus separately at a 1.0 °C/min cooling rate.These relatively fast cooling rates in this case were necessary due to the possible decomposition of the compound in extended contact with water at elevated temperatures.An additional data point on the LLPS curve was determined separately using a cold stage microscopy technique, which agreed with the LLPS data trend obtained with Crystal 16.

Idebenone -Hexane/Methylene Chloride
A third solute system we consider here is the API idebenone crystallizing in a mixture of 7:1 (by volume) hexane/methylene chloride, as reported by Lu et al. [5,6].Lu et al. determined the solubility of idebenone in the above solvent mixture by a gravimetric method.The oiling-out points reported in [5] were obtained by cooling the solution of a pre-determined composition at 0.1 °C/min.The authors observed that when the extent of oiling out was small, normal crystallization occurred, presumably through oil droplets.However, at high degrees of oiling out, the droplets coalesced to form a thick gel phase, which did not crystallize readily.The authors also observed that for the cases of crystallization mediated through LLPS, impurity rejection was suboptimal.

C₃₅H₄₁Cl₂N₃O₂ -Ethanol/Water
Another system we considered was reported by Veesler and co-workers [47][48][49][50][51] for an API with a molecular formula of C₃₅H₄₁Cl₂N₃O₂ crystallizing in ethanol/water mixtures.The solubility data of form I and form II crystals for this compound were determined using a temperature bracketing technique that monitored the growth and dissolution of small particles [47].The binodal data were determined using static light scattering to monitor the occurrence of the second liquid phase in the solution [50], and also by turbidity and FBRM probes [49].These authors also report data on the spinodal boundary extracted from the light scattering intensities through appropriate extrapolations [50].The available data for this system are in the temperature range of 10 to 60 °C.

Vanillin -1-Propanol/Water
The phase behavior of vanillin crystallizing in a 1-propanol/water mixture was reported by Zhao et al. [52].These authors obtained the solubility data using a gravimetric method.Oiling out in their experiments was detected with the help of FBRM/PVM instruments.The oiling-out points reported in [52] were determined by cooling the solution of a pre-determined composition at 0.1 °C/min.Available data for this solute-solvent system are in the temperature range of 5 to 35 °C.

Data Extraction and Processing
The original data reported in the above publications are available in the form of temperatureconcentration plots for various phase boundaries.These experimental data points were extracted from the published figures using Engauge Digitizer (v 10.4), a free utility [53].The density of the solvent mixture in each case was determined using standard density data for the solvents [54] and a mixing rule that considers ideal mixing (i.e., volume additivity holds).The volume fractions of the solutes were estimated using simple mass balance calculations with the known densities of the solvent (temperature dependent) and solute (temperature independent).For each solute studied, a consistent value for the molecular diameter (σ) as an equivalent spherical diameter was calculated by assuming the crystals to be of an fcc lattice with a packing fraction of 0.74, as suggested in [44].From this consideration, we obtained σ from the solid density and molecular weight using the relation

Results
The results of the calculations for the pyraclostrobin -isopropanol/cyclohexane system are shown in Figure 1.Table 1 lists the fit parameters that best describe the solubility data for this system obtained through nonlinear regression.Note that in Figure 1, only the solubility curve (curve 1 passing through the open circles) is a curve fit.The other curves were calculated (i.e., not fit) using the square-well interaction parameters obtained from the solubility fit (Table 1).The good agreement between the experimental oiling-out data points and the calculated curve 2 shows that the modeling approach described in Section 2 captures the essential features of the phase equilibria for this system.We also note that the various LLPS curves predicted for this system using Equation (11) (dotted lines) and Equation ( 12) (solids lines) for ηeff are close to each other and either of those curves can be used as a guideline in crystallization process development.

Table 1.
Best fit parameters for the pyraclostrobin -isopropanol/cyclohexanone system estimated through fitting the experimental data of Li et al. [4].Parameters obtained using both the models for the first perturbative term of the Helmholtz free energy (i.e., ηeff through Equations ( 11) and ( 12)), are listed.Also given are the estimated critical temperature (Tc) and critical concentration (Wc).The calculations were performed with σ = 0.896 nm.16) and ( 17).The square-well interaction parameters (listed in Table 1) obtained from this fit were used to further calculate the binodal (curve 2), spinodal (curve 3), and the gelation (curve 4) boundaries.The dotted lines are the curves obtained using Equation ( 11) and the solid lines are those obtained with Equation ( 12), respectively, for calculating ηeff.The dotted and solid lines for curve 4 are almost identical.
Figure 2 shows the results of the modeling effort for the compound Z -methanol/water system.Here, we also note good agreement between the LLPS line (the binodal calculated using the squarewell parameters regressed from the solubility data) and the experimental data.Table 2 gives the bestfit parameters for this case.From Table 2, we note that for this system, even though the estimated parameters are nearly equal in both the cases of using Equations ( 11) and ( 12) for ηeff, the critical temperatures predicted by the models are significantly different.The critical concentration, however, is predicted to be about the same.The open circles and diamonds indicate the experimental data points for solubility and oiling out (binodal), respectively, adapted from [11].Curve 1 represents the model (μL = μS) fit to solubility data using Equations ( 16) and (17).The square-well interaction parameters (listed in Table 2) obtained from this fit were used to further calculate the binodal (curve 2), spinodal (curve 3), and the gelation (curve 4) boundaries.The dotted lines are the curves obtained using Equation ( 11) and the solid lines The results of the calculations for the ibedenone -hexane/methylene chloride system are shown in Figure 3.The values of square-well parameters that best describe the solubility data for this case are given in Table 3.As seen from Figure 3, in this case also, the model captures the trends in the experimental data very well.Note that for this solute-solvent system, the solubility curve appears to penetrate the oiling-out region for temperatures greater than 50 °C and emerge on the other side.This aspect of the model predictions is discussed further in Section 5.  16) and (17).The square-well interaction parameters (listed in Table 3) obtained from this fit were used to further calculate the binodal (curve 2), spinodal (curve 3), and gelation (curve 4) boundaries.The dotted lines are the curves obtained using Equation (11) and the solid lines are those obtained with Equation ( 12), respectively, for calculating ηeff.The dotted and solid lines for curve 4 are almost identical.
Table 2. Best fit parameters for the compound Z -methanol/water system estimated through fitting the experimental data of Bhamidi et al. [11].Parameters obtained using both the models for the first perturbative term of the Helmholtz free energy (i.e., ηeff through Equations ( 11) and ( 12)) are listed.Also given are the estimated critical temperature (Tc) and critical concentration (Wc).The calculations were performed with σ = 0.664 nm.  3. Best fit parameters for the idebenone -hexane/ methylene chloride (7:1 by volume) system estimated through fitting the experimental data of Lu et al. [5].Parameters obtained using both the models for the first perturbative term of the Helmholtz free energy (i.e., ηeff through Equations ( 11) and ( 12)), are listed.Also given are the estimated critical temperature (Tc) and critical concentration (Wc).The calculations were performed with σ = 0.872 nm.Now, let us examine the performance of the model in capturing the phase behaviors of C₃₅H₄₁Cl₂N₃O₂ in the ethanol/water mixture, studied by Veelser and co-workers [47][48][49][50][51], and of vanillin in 1-propanol/water, reported by Zhao et al. [52].Figure 4 shows the results of the modeling effort for these systems.Table 4 presents the values of the molecular interaction parameters obtained from regressing the solubility data for these systems.One notes from the plots in Figure 4 that while the theoretical framework discussed describes the solid-liquid equilibrium (SLE) aspects of the systems very well, the resulting molecular interaction parameters fail to capture the experimentally observed LLE boundaries.This aspect is further discussed below.Table 4. Best fit parameters for the C₃₅H₄₁Cl₂N₃O₂ -ethanol/water system estimated through fitting the experimental data of Veesler et al. [48], and for the vanillin -1-propanol/water system estimated through fitting the experimental data of Zhao et al. [52].Parameters obtained using both the models for the first perturbative term of the Helmholtz free energy (i.e., ηeff through Equations ( 11) and ( 12)), are listed.Also given are the estimated critical temperature (Tc) and critical concentration (Wc).The calculations were performed with σ = 1.044 nm for C₃₅H₄₁Cl₂N₃O₂ and 0.697 nm for vanillin.In (a), the triangles are the spinodal points reported in [48] obtained through extrapolation of the light scattering intensity data [50].Curve 1 represents the model (μL = μS) fit to solubility data using Equations ( 16) and ( 17).The square-well interaction parameters (listed in Table 4) obtained from these fits were used to further calculate the binodal (curve 2) and spinodal (curve 3) boundaries.The dotted

(a) (b)
Replace with the high resolution figure submitted lines are the curves obtained using Equation (11) and the solid lines are those obtained with Equation ( 12), respectively, for calculating ηeff.

Discussion
Several experimental studies in the literature have focused on the influence of the rate of cooling on the observed oiling-out boundary [4,5,7].The fact that these studies observed the oiling-out temperatures to be dependent on the cooling rate suggests that oiling out in these experiments was due to the penetration of a binodal.This is because phase separation between a binodal and a spinodal occurs through nucleation, which is an activated process similar to solid crystal nucleation.The phase separation that occurs beyond a spinodal is spontaneous and is usually uncontrollable.The binodal observed during a cooling crystallization may hence have a cooling rate dependence (similar to a metastable zone width).However, the true binodal represents a (local) thermodynamic equilibrium, and hence is not influenced by the rate of generation of supersaturation.The theoretical framework described in the present work provides one with estimates of the true binodal and spinodal.

Relevance of the Spinodal and Gelation Lines to Process Development
When LLPS occurs in a crystallizing system, the resulting oil phase may stay suspended to ultimately produce crystals, or the phase may form thick gels that coat the vessel internals and lead to operational problems [5,6].The predictive approach discussed in this work offers guidelines for a process development scientist in this regard.
As mentioned earlier, the oiling-out boundary observed during cooling crystallization usually represents a binodal because penetration of the spinodal envelope requires significantly faster cool down of solutions than is achievable experimentally.However, such deep quenches are possible during reactive crystallization, in which a fast chemical reaction can rapidly increase the concentration of the solute in a solvent at a given temperature, such that the solution becomes unstable to a homogeneous composition [11].In this case, a spinodal becomes relevant to the process design, as this boundary sets the limits on the extent of the chemical reaction possible before one loses control over the crystallization process.
During oiling out, the solute partitions into two liquid phases, one that has a high concentration of the solute and the other being relatively lean in solute.The relative amounts of each phase are governed by the "lever rule" [55], i.e., the relative amounts of these phases are inversely proportional to the difference between the concentration of solute in that particular phase and the concentration of the solute before phase separation.Thus, given the shape of the binodal (and spinodal), it is evident that when LLPS occurs at low temperatures and low solute concentrations, only a small amount of the oil phase forms as a dispersed phase relative to the continuous lean phase.At high solute concentrations and low temperatures, LLPS will result in a large amount of the solute-rich phase, which may tend to become the continuous phase, trapping the lean liquid in pockets.At high temperatures and high concentrations, the relative amounts of the dense and lean phases are comparable.Whether a coating on the reactor walls will result from the dense phase formed depends on the solute concentration.If the solute concentration in the dense liquid phase exceeds the estimated gelation concentration, one may expect to encounter gel phase formation and resultant operational issues during crystallization.
Thus, a contextual understanding of the spinodal and gelation boundaries greatly helps a process scientist/engineer in developing a crystallization protocol for solutes prone to oiling out.

Limitations of the Current Modeling Approach
For the ibedenone -hexane/methylene chloride system (Figure 3), we observed that the solubility curve predicted by the model extrapolates into the LLPS region and crosses the liquidliquid equilibrium (LLE) region at temperatures greater than 50 °C.Whether this predicted trend is physical may be debated.The observation by Lu et al. [5] that at these temperatures, the oiled-out phases never produced crystals indicates that the predicted phase behavior may be real.Lu et al.'s report suggests that the system exhibits a stable liquid-liquid equilibrium (LLE) at high temperatures, as observed for the cases of vanillin -water [56], ethyl-2-ethoxy-3-(4-hydroxyphenyl)propanoate (EEHP) -cyclohexane [57], and L-menthol-water [3].One also notes that the solvents involved in Lu et al.'s work-hexane and methylene chloride-are low-boilers, and hence a significant amount of solvent would have evaporated (in an open system) or present in the vapor phase (in a closed system) at temperatures >50 °C.Such high temperatures are not practically relevant in a typical process development exercise.
For C₃₅H₄₁Cl₂N₃O₂ and vanillin (Figure 4), the predicted liquid-liquid phase boundaries do not agree with the experimental data.Several reasons may be contemplated for these results.Stable LLE, as seen in Figure 4b for vanillin, typically occurs when the solute molecules interact with significantly long-range attractions [22,28,32].Such long-range molecular interactions may not be captured by a square-well potential sufficiently accurately.The square-well potential provides a reasonable representation of the short-range intermolecular interactions for which the LLE is metastable with respect to SLE.Clearly, intermolecular interactions are more complex than captured by a square-well model.
Another issue is that the theoretical concepts that underlie Equation ( 17) are different from, and are less rigorous than, those used to describe the thermodynamics of the fluid phase.For example, the simple expression used for the solid chemical potential does not account for volumetric changes in the solid due to changes in temperature.It also fails to adequately describe chemical potential differences between various polymorphs, which arise from differences in the lattice arrangement of solute molecules in polymorphic crystals.The simplifying assumptions behind the description of solid chemical potential can lead to inaccuracies in the calculation of phase equilibria.If one intends to account for these shortcomings and improve the accuracy of the estimated solid chemical potential, one needs to bring in the radial distribution function of the solid phase into the calculation (see [20]), at which point the procedure becomes mathematically complicated and tedious, and loses its attractiveness.
Also, we note that the model treats the solute as a simple fluid acting alone and ignores the presence of the solvent.Thus, in essence, the ε and λ we obtain by regressing the solubility data represent an "effective" strength and range of interaction with solvent molecules in the background.This aspect may become significant when the LLPS results in a change in solvent composition in the resulting liquid phases.For a mixed solvent system, if one solvent selectively partitions into one of the liquid phases upon LLPS, the nature of the molecular interactions cannot remain the same in both the phases due to a change in the solvent backdrop experienced by the solute molecules.In other words, one cannot describe the chemical potential of the two liquid phases using single values for ε and λ.Indeed, such selective solvent partitioning has been reported for these systems through the ternary equilibria determined for C₃₅H₄₁Cl₂N₃O₂ [50] and for vanillin [58].This solvent partitioning effect may be significant for solvents that readily associate, such as alcohol(s) and water.Thus, the LLPS model discussed here may have a limited use in systems that contain a mixture of solvents (that associate readily), with each solvent being in considerable proportion.For solvent mixtures in which the component solvents are organic and/or non-associative, or if one component solvent is present in a small portion, the modeling approach may work better.This thought is supported by Figures 1-3, in which the model captured the physics of the systems very well despite a mixture of solvents being present.

Concluding Remarks
While the phenomenon of oiling out in solution crystallization is well known, most studies that focus on the development of crystallization processes approach this problem somewhat empirically.Typically, the temperature at which a homogeneous solution with fixed composition oils out during cooling is determined (sometimes as a function of the cooling rate) experimentally.The thermodynamic relevance of this experimental LLPS (oiling out) boundary in relation to the possible binodal and spinodal curves is not usually emphasized and, in some cases, is completely missed.This work aims to demystify the oiling-out phenomenon often found daunting by some process engineers by placing it in its thermodynamic context.This work also introduces a broad class of literature on phase diagram modeling from first principles to a reader who may be unfamiliar with the approach.
The solubility and binodal boundaries of compounds are typically modeled through activity coefficient methods, in which the general solubility relation (van't Hoff equation) is used in conjunction with activity coefficient models (e.g., NRTL, PC-SAFT models) [56,59].In addition to the solubility (SLE) data, this approach requires the knowledge of a few physical properties of the solute and solvent systems, such as the melting point and enthalpy of fusion of the solute, the difference in the heat capacity of the solute in liquid and solid phases at the melting point, etc.This information typically is not available readily and one is required to perform additional solid-state characterization for the solute crystals.Moreover, while these activity coefficient models can predict a liquid-liquid co-existence boundary, prediction of the spinodal and gelation boundaries using these methods is difficult.The method we presented here adopts a more fundamental approach to the characterization of LLPS boundaries than considered by the (semi-empirical) activity coefficient models.It uses only the solubility data (which can be easily obtained experimentally) to extract the molecular interaction parameters necessary to estimate the other phase boundaries.
The method discussed here falls under the general category of "equation of state" (EOS) models, in which expressions for the Helmholtz free energy and compressibility factor (Equations ( 2) and ( 16)) are derived starting from suitable expressions for the intermolecular potential.The ensuing procedure of equating the chemical potentials and pressures of the phases in coexistence to calculate the phase boundaries follows standard thermodynamic concepts.The real challenge in executing these calculations often arises from the theoretical formalism used in the definition of various terms in Equation (2).The procedure described here using a square-well interaction model can also be readily adapted to other molecular potentials, such as a Yukawa interaction model or a Lennard-Jones model (for example, see [19] and [20]).Note, however, that we are interested in an LLPS that is metastable with respect to the crystalline phase.Such behavior is not captured by all molecular potentials.For example, a Lennard-Jones fluid does not exhibit a metastable LLPS [28].In the end, all these molecular potential models are abstractions of convenience to represent the real interactions between molecules.An appropriate choice for the pair potential is the one that can capture the observed thermodynamic features of a real system of interest and makes the resulting calculations tenable.
Clearly, the choice of the intermolecular potential and the perturbation scheme influence the predicted SLE and LLE boundaries.The various theoretical formalisms developed over the decades have attempted to improve the accuracy of model predictions by incorporating increasingly sophisticated theoretical arguments.Many such advanced theories are still being actively developed.In the context of industrial process development, one is more interested in a readily usable and reasonably accurate modeling approach than a hard-to-execute, but more sophisticated, calculation procedure.Our choices, of the square-well potential to represent the intermolecular interactions (Equation ( 1)), the SAFT-VR model for the Helmholtz free energy (Equations ( 10)-( 14)), and the simple equation for the chemical potential of a square-well solid (Equation ( 17)) are all influenced by this philosophy.Similarly, the choice of a linear relationship between ε/k and T is also arbitrary and is based on the values of ε/k for ibuprofen crystallizing in ethanol/water mixtures tabulated in [45].These choices influence the predicted phase boundaries, and the accuracy of any prediction with these or other choices can only be established through extensive experimentation.
Note that the goal of this paper is not to present an accurate and original theoretical model for the prediction of oiling-out boundaries.Our objective is to provide an industrial scientist the thermodynamic context of oiling out and related phenomena and a practical approach to estimate the phase boundaries.To the best of our knowledge, a unified approach based on thermodynamic first principles that allows a process scientist to predict the binodal, spinodal, and gelation boundaries for solutes has not been presented before.Through this work, we have attempted to bridge this gap.We readily recognize the limitations of the theoretical concepts that underlie the model equations discussed.Given the simplicity of the modeling approach, one is advised to obtain a few oiling-out data points experimentally when relevant and employ the calculation procedure judiciously to obtain working models for various phase boundaries for the systems of interest.
temperature, T, to obtain the volume fractions, I II and η η , of the solute in the two liquid phases I and II at (local) equilibrium.This solution procedure can be tricky because of the existence of a trivial solution at I II = η η .This complexity is tackled as follows.First, we determine the critical temperature, Tc, and critical volume fraction, ηc, at which the binodal and spinodal meet.At this upper consolute temperature (UCST) and composition, the derivatives of the chemical potential, w.r.t, to the volume fraction must vanish.From this condition, we obtain Tc and ηc by solving  of α₀, α₁, and λ regressed from the solubility data.Once Tc and ηc are obtained, we reduce the temperature slightly and calculate I , δη is a small value, typically of the order of 0.001.The freshly obtained I II and η ρs is the density of the crystal, MW is the molecular weight of the solute, and Nav is the Avogadro number.We implemented all the phase diagram calculations in Wolfram Mathematica (v11.3,Wolfram Research, Champaign, IL, U.S.A).While the curve fits were obtained by considering the errors in solubility in terms of volume fractions, the modeling results are shown below with concentrations expressed in wt% for convenience.

Figure 1 .
Figure 1.Experimental data and model predictions for the pyraclostrobin -isopropanol/cyclohexane system.The open circles and diamonds indicate the experimental data points for solubility and oiling out (binodal), respectively, adapted from [4].Curve 1 represents the model (μL = μS) fit to solubility data using Equations (16) and (17).The square-well interaction parameters (listed in Table1) obtained from this fit were used to further calculate the binodal (curve 2), spinodal (curve 3), and the gelation (curve 4) boundaries.The dotted lines are the curves obtained using Equation (11) and the solid lines are those obtained with Equation (12), respectively, for calculating ηeff.The dotted and solid lines for curve 4 are almost identical.

Figure 2 .
Figure 2. Experimental data and model predictions for the compound Z -methanol/water system showing (a) the complete phase diagram, and (b) magnified view of the experimental data region.The open circles and diamonds indicate the experimental data points for solubility and oiling out (binodal), respectively, adapted from[11].Curve 1 represents the model (μL = μS) fit to solubility data using Equations (16) and(17).The square-well interaction parameters (listed in Table2) obtained from this fit were used to further calculate the binodal (curve 2), spinodal (curve 3), and the gelation (curve 4) boundaries.The dotted lines are the curves obtained using Equation(11) and the solid lines Replace with the high resolution figure submittedReplace with the high resolution figure submitted are those obtained with Equation (12), respectively, for calculating ηeff.The dotted and solid lines for curve 4 are almost identical.

Figure 3 .
Figure 3. Experimental data and model predictions for the idebenone -hexane/methylene chloride (7:1 by volume) system showing (a) the complete phase diagram, and (b) magnified view of the experimental data region.The open circles and diamonds indicate the experimental data points for solubility and oiling out (binodal), respectively, adapted from [5].Curve 1 represents the model (μL = μS) fit to solubility data using Equations (16) and(17).The square-well interaction parameters (listed in Table3) obtained from this fit were used to further calculate the binodal (curve 2), spinodal (curve 3), and gelation (curve 4) boundaries.The dotted lines are the curves obtained using Equation(11) and the solid lines are those obtained with Equation (12), respectively, for calculating ηeff.The dotted and solid lines for curve 4 are almost identical.

Figure 4 .
Figure 4. Experimental data and model predictions for (a) C₃₅H₄₁Cl₂N₃O₂ (form I) -ethanol/water system[48] and (b) for vanillin in 1-propanol/water[52].In both (a) and (b), the open circles and diamonds indicate the experimental data points for solubility and oiling out (binodal), respectively.In (a), the triangles are the spinodal points reported in[48] obtained through extrapolation of the light scattering intensity data[50].Curve 1 represents the model (μL = μS) fit to solubility data using Equations (16) and(17).The square-well interaction parameters (listed in Table4) obtained from these fits were used to further calculate the binodal (curve 2) and spinodal (curve 3) boundaries.The dotted Replace with the high resolution figuresubmitted