Modeling size controlled nanoparticle precipitation with the co-solvency method by spinodal decomposition

Abstract

The co-solvency method is a method for the size controlled preparation of nanoparticles like polymersomes, where a poor co-solvent is mixed into a homogeneous copolymer solution to trigger precipitation of the polymer. The size of the resulting particles is determined by the rate of co-solvent addition. We use the Cahn–Hilliard equation with a Flory–Huggins free energy model to describe the precipitation of a polymer under changing solvent quality by applying a time dependent Flory–Huggins interaction parameter. The analysis focuses on the characteristic size R of polymer aggregates that form during the initial spinodal decomposition stage, and especially on how R depends on the rate s of solvent quality change. Both numerical results and a perturbation analysis predict a power law dependence R ∼ s^−⅙, which is in agreement with power laws for the final particle sizes that have been reported from experiments and molecular dynamics simulations. Hence, our model results suggest that the nanoparticle size in size-controlled precipitation is essentially determined during the spinodal decomposition stage.

---

1. Introduction

Because of their great potential in nano- and biotechnology, polymeric nanoparticles such as polymersomes have attracted growing interest during the last decades.^1,2 One possible application can be found in the field of drug delivery, where they serve as transport vehicles for medication.³ A crucial property of such transport vehicles is their size, as it does not only determine their loading capacity, but also the composition of their protein corona in blood, which affects the retention times in the circulatory system.⁴ Furthermore, the nanoparticle size plays a critical role in passive targeting of tumors, which is based on the enhanced permeability and retention effect.⁵

A method to prepare polymersomes of a particular size is the co-solvency method or flash nanoprecipitation: a poor co-solvent is mixed into an initially homogeneous solution of a good solvent and a block copolymer to induce particle formation via self-assembly of the polymer.^6,7 The co-solvency method can be implemented in different ways. A straightforward approach is to add co-solvent by drop injection to a polymer solution. In this method the size of the produced nanoparticles depends on the rate of co-solvent addition.⁸ Thiermann et al. fed the co-solvent and the copolymer solution into continuous flow multilamination micro mixers and observed that the size of the synthesized particles decreases with an increasing flow rate, which connects nanoparticle size to an easily adjustable parameter of the experimental setup.⁹ The micro mixer approach has several advantages as it yields narrower size distributions and can be done without additional steps like membrane extrusion to achieve acceptably low polydispersities.

A superficial comparison between the two realizations of the co-solvency method suggests that the size of the produced nanoparticles depends on a completely different quantity in both cases: the rate of co-solvent addition on the one hand and the flow rate on the other. However, due to the special design of multilamination micro mixers, an increase in flow rate decreases the mixing time of liquids that are fed into its inputs.¹⁰ Changing the flow rate indirectly changes the mixing rate. Hence, in both cases particle sizes are found to depend on the mixing rates of the co-solvent and the copolymer solution.

To gain a deeper understanding of the size controlled preparation of polymeric nanoparticles with the co-solvency method, one must analyze how different rates of co-solvent addition affect the particle formation and how the particle size depends on the rate of co-solvent addition. In this article we consider the earliest stage of particle formation, the spinodal decomposition of an oversaturated polymer solution. We present a simple phase field model that can be used to determine the size of aggregates in that stage. The model is an extension of the popular Cahn–Hilliard model for the dynamics of phase separation, where we represent the effect of mixing solvent and co-solvent in an effective manner by using time dependent interaction parameters. Using numerical simulations of an idealized mixing process, we show that the model can reproduce the dependence of particle size on mixing speed observed experimentally.

Our article is organized as follows: we present the model in Section 2 and describe the simulation method in Section 3. In Section 4, we explain the evaluation method, show our simulation results, compare them to a perturbation theory, and discuss scaling laws both have in common. These scaling laws are then compared to experimentally observed power laws in Section 5. We summarize and conclude in Section 6.

2. Theoretical model

Experimentally, three different components are involved in the co-solvency method of nanoparticle synthesis: a polymer, a good solvent and a poor or selective co-solvent. The precipitation of the polymer is triggered and influenced by the continuous addition of co-solvent into the polymer solution – i.e., by solvent mixing. We assume that solvent mixing is fast on the time scales of the polymer phase separation and that the main effect of solvent mixing is a change of ‘mean solvent quality’ from ‘good’ to ‘poor’. Thus, we incorporate solvent mixing by only taking into account the change in solvent quality: the three component system from the experiment is modeled by a two component system containing a polymer and only one solvent, which changes its quality over time. More specifically, at any given time we describe the momentary solvent mixture by one single effective solvent with an associated interaction parameter χ at a polymer–solvent contact. The addition of co-solvent into the solvent mixture is modeled by a temporal increase of χ.

Close to a homogeneous ground state isothermal phase separation of incompressible binary mixtures in a fixed volume can be modeled by the Cahn–Hilliard equation. It describes the local evolution of a globally conserved dimensionless composition field u, for example the volume fraction of one component. The Cahn–Hilliard equation is a special case of the generalized diffusion equation

ζ represents the scale of the mobility and M(u) describes its dependence on the composition. v is the volume of a polymer or solvent segment. Here we will use the “degenerate mobility” M(u) = u·(1 − u), which is suitable for the description of composition currents in incompressible mixtures.^12,13 δF/δu is the functional derivative of the free energy functional and can be interpreted as a chemical potential. With the free energy functional F proposed by Cahn and Hilliard,¹¹(d is the spatial dimension) one obtainsHere, f(u) is the free energy per segment in a homogeneous system, ½λ(∇u)² represents surface contributions, and λ is the gradient energy parameter. Insertion of eqn (3) into eqn (1) yields the Cahn–Hilliard equationIn this article we specifically consider polymeric systems. Hence we choose ζ = D/k_BT and λ = R_g²·k_BT with the segment diffusion coefficient D = D_pN, the Boltzmann factor 1/k_BT and the radius of gyration R_g. D_p is the diffusion coefficient of a polymer chain composed of N segments. The expression for λ is an approximation to the gradient energy parameter for a binary homopolymer solvent mixture, which holds in the weak segregation limit where concentration gradients are weak.¹⁴ We also usewhich is free energy per segment from the Flory–Huggins solution theory,^15,16 where χ is the Flory–Huggins interaction parameter mentioned earlier.

To describe phase separation during solvent mixing, we assume that χ in eqn (5) explicitly depends on time, i.e. χ = χ(t), and increases from χ(0) = χ₀ to χ_max. We choose χ₀ to be the spinodal interaction parameter, which is the value of χ where a homogeneous system becomes unstable. It depends on the mean polymer volume fraction in the system,

as well as N and is defined by the condition , which yieldsχ_max > χ₀ is a constant and we consider a situation where χ(t) grows linearly aswith t_max = (χ_max − χ₀)/s.

In the following, all lengths will be given in units of , all times in units of t₀ = l₀²/(ζk_BT) = R_g²/D, and all energies in units of k_BT(R^d_g/v). Eqn (4) can then be rewritten as

and the derivative of the free energy becomes In experiments, polymersomes are typically formed from amphiphilic diblock copolymers and stabilized by the hydrophilic block. The model presented in this article neither incorporates the stabilization effects from copolymers nor is it able to describe an internal structure of polymer aggregates, because eqn (5) is restricted to homopolymers. Thus, the equilibrium state will always correspond to macroscopic phase separation. However, simulations of more detailed models have shown that nanoparticle self-assembly is initially dominated by the formation of unstructured droplets, and that the number of droplets after the initial stage largely determines the final number and size of particles.¹⁷ The system defined by eqn (5) represents the simplest model system that reproduces this early stage of particle assembly.

In this article we restrict our investigations to the very early stages of phase separation, where the first patterns in the concentration profile form and where the gradients in the composition field are still small, which motivates the application of eqn (4) with our choice of λ. In the context of different possible mechanisms that lead to the formation of structured copolymer–nanoparticles,^17,18 we focus on the spinodal decomposition before the first micelles appear. In these very early stages, the self-assembly should be driven mainly by the energetically unfavorable interaction between the co-solvent and the co-solvent-phobic block of the polymer, which leads to typical ‘Cahn–Hilliard-type’ spinodal decomposition patterns in the concentration profiles.¹⁹ The solvent-philic block of the polymer, which is often incompatible with the other one, is mainly responsible for internal structure formation in aggregates once they have formed. So if the internal structuring does not significantly change the size of an aggregate, the substitution of the copolymer by a homopolymer of its hydrophobic block might still yield approximate results for its size. We shall see below that eqn (5) is indeed sufficient to describe the relation between particle sizes and mixing rates in the early stages of mixing.

A very recent publication also shows that it is possible to produce homopolymer particles stabilized by surface charges.²⁰ Besides an experimental part it also contains molecular dynamics simulations, where solvent mixing is modeled by a time dependent strength of the repulsive force. They observe very similar scaling laws as we do. Another work from the same group applied time dependent repulsive forces to dissipative particle dynamics simulations for copolymers.²¹ Although they only slightly touched the issue of particle size dependence on mixing time their curves look also similar to ours and the ones from.²⁰ Thus, the model presented in this article reproduces important features observed in much more complex particle models. Its simple structure allows a perturbation treatment and we will see that typical scaling laws observed in our simulations, the molecular dynamics simulations and the experiments are already inherent in the perturbation theory, which might pave the way for semi analytical approaches. In addition, the phase field model allows to study slow mixing processes with characteristic mixing times in the range of milliseconds or more, whereas molecular dynamics simulations are limited to microseconds.²⁰

There also exists a pinning effect of structures,²² which is caused by viscoelasticity in systems with asymmetric molecular dynamics, i.e. polymer solutions. It might also affect particle sizes and there are models that incorporate this effect.²³ However, the present model does not, because pinning did not occur in the experiments⁹ we aim to describe.

Since we focus on situations where the particle formation is a thermodynamically driven process (initiated by spinodal decomposition), which does not involve a thermally activated crossing of free energy barriers as in nucleation theory, we do not include thermal noise in our theoretical model, eqn (1). This corresponds to the limit v → 0 in eqn (1) and (2) (thermal noise would scale with ) and is also motivated by the fact that the relative thermal fluctuations are generally small in polymeric systems.

For simplicity, we call the first aggregates of well-defined shape that emerge during spinodal decomposition ‘particles’. How we exactly define these particles and how we determine their size is described in the discussion in Section 4.

3. Simulation method

Eqn (7)–(9) with M(u) = u(1 − u) constitute the model equations. The applied scheme is very similar to the one used by Zhu et al.,²⁴ which is a first order time accurate pseudo spectral method, and any Fourier transform was calculated by the FFTW library.²⁵ The domains are boxes [0, L_b)^d × [Doublestruck R]₀⁺ with d = 2, 3 and periodic boundary conditions. As initial conditions we use uniformly distributed random perturbations in the interval [u̅ − 0.001, u̅ + 0.001] generated by the Mersenne twister.²⁶ All numerical results are averages over 5 simulation runs with different initial perturbations and we performed simulations in both 2 and 3 dimensions to check the influence of dimensionality. As it will turn out, the particular dimension plays only a minor role, which allows investigations of main dependencies in 2D to speed up computation times.

The adjustable physical parameters in the model are N, s, χ_max and u̅. The slope s in eqn (7) parameterizes the rate of solvent quality change. It is varied to investigate the effect of different solvent mixing rates on phase separation, while the three remaining parameters are kept constant. In this article, we will mostly study a model with parameters set to u̅ = 0.1, N = 14, and χ_max = 2. Simulations for more realistic parameters are shown in Section 5. In 2D we used 400 × 400 grid points with a lattice constant 0.25, and set the time step to 0.005.

Using that lattice constant assures that the spatial resolution does not limit the smallest particles we encounter in our simulations (which is the particle size at χ_max). For constant interaction this problem could also be approached by a rescaling of the spatial coordinate that involves the quench depth Δχ.²⁷ However, this might introduce numerical artifacts that lead to unphysical pinning close to the spinodal. Even though these artifacts can be avoided by proper normalization,^28,29 we did not scale the spatial coordinate by Δχ because in our case it depends on time, meaning the system size would exhibit a temporal change if we kept our lattice size and number of grid points constant as is customary in simulations. We should also note that we did not encounter any pinning artifact in our simulations either. In 3D we used 64 × 64 × 64 grid points and a lattice constant of 1.

4. Results and discussion

4.1. Qualitative characterization of the demixing process

In all simulation runs the phase separation proceeds in a similar manner as in the case of constant interaction parameters χ. In the first stage, termed spinodal decomposition, a bicontinuous pattern emerges in the concentration profiles and initially grows and coarsens on a relatively fast time scale, until droplets with well-defined interfaces have formed (see examples in Fig. 2). In the second stage, called Ostwald ripening, the droplet pattern coarsens on a very slow time scale. As stated at the end of Section 2, we focus on the spinodal decomposition stage. This restriction requires us to identify the crossover time between spinodal decomposition and Ostwald ripening. To this end, we use a procedure proposed by Sofonea and Mecke, which is based on Minkowski measures.³⁰ Minkowski measures are a complete set of additive motion-invariant measures for unions of convex sets. Each measure assigns one real number to any polymer volume fraction profile depending on its morphology. Since the morphology of concentration profiles during phase separation depends on time, the Minkowski measures also do. One of these measures, from here on denoted by C, is the total boundary length of the union over all subsets in [0, L_b)^d where the polymer volume fraction u exceeds a predefined threshold u_th. During spinodal decomposition, polymer aggregates form on a fast time scale leading to a rapid temporal increase of C and during Ostwald ripening, polymer aggregates merge on a large time scale leading to a slow decrease of C. These two characteristic regimes can be seen in Fig. 1, which shows C as a function of time for two examples discussed below. The regimes are separated by a maximum of C and the corresponding time is called the transition time t_tr.³⁰ Therefore, spinodal decomposition dominates for t < t_tr and Ostwald ripening for t > t_tr. The remaining Minkowski measures yield equivalent estimates for the transition time.³⁰ To calculate the Minkowski measures we use the algorithm proposed by Mantz et al.³¹

Fig. 1 Time series of the Minkowski measure C (see text) for slopes s = 5 × 10⁻⁵ (a) and s = 5 × 10⁻³ (b). The time when C reaches its maximum is defined as the transition time. The threshold value u_th is set to 0.3. The insets show a magnification of the regions in the dashed rectangles.

In the following, we first discuss exemplarily the effect of a time dependent interaction parameter on spinodal decomposition by comparing the results from s = 5 × 10⁻⁵ and s = 5 × 10⁻³. Fig. 1 shows the corresponding time series of C. The transition time obviously depends on s. Hence, spinodal decomposition happens faster for large values of s. Fig. 2 illustrates how different values of s affect the morphology of the polymer volume fraction profiles during spinodal decomposition. The upper panel (Fig. 2(a), (c) and (e)) and the lower panel (Fig. 2(b), (d) and (f)) show temporal evolutions of the same initial polymer volume fraction profile for different growth rates s. At t = 10 the volume fraction profiles look very similar (Fig. 2(a) and (b)). At t = 200, however, they deviate significantly from each other (Fig. 2(c) and (d)) and at the end of the spinodal decomposition there are significantly smaller droplets for the larger quench rate (Fig. 2(e) and (f)). Hence, we see that the time dependence in the interaction parameter does not only affect transition times but also the length scales of structures during spinodal decomposition. We can rationalize this observation as follows: with incasing χ, one reaches deeper into the miscibility gap and the characteristic wavelength of the most unstable mode decreases. If χ(t) increases very slowly, the initially unstable modes have time to grow and dominate also the later stages of demixing. If χ(t) increases more rapidly, modes with smaller wavelengths take over and determine the final structure. Indeed, Fig. 2 demonstrates that the characteristic length scale of patterns in the initial stage of demixing (Fig. 2(c) and (d)) is larger than the characteristic length scale of the final droplet pattern (Fig. 2(e) and (f)). We will analyze this effect at a more quantitative level further below in Section 4.3.

Fig. 2 Spatial distribution of polymer volume fraction u(x,y,t) (color coded) in the domain [0, L_b)² at different times t (t = 10, t = 200, and t = t_tr) during spinodal decomposition for s = 5 × 10⁻⁵ (upper panel: (a), (c), (e)) and s = 5 × 10⁻³ (lower panel: (b), (d), (f)). The transition time t_tr depends on s and marks the time at which the behavior crosses over from spinodal decomposition to Ostwald ripening. The color coding is different for every snapshot and chosen such that the smallest value is blue and the largest dark red.

In general, bicontinuous patterns are favored for a larger range of composition variables u̅ in 3D than in 2D, but for our set of parameters we observed droplets in both dimensions.

4.2. Definition of particles and quantitative determination of particle size

To examine the length scales in the volume fraction profiles more quantitatively, we evaluate the normalized radially averaged structure factor²⁴
where S_c is the absolute value from the radial average of
Here S is the discrete Fourier transform of composition correlations in d dimensions with wave vectors and the summation is carried out over all grid points r⃑, r⃑′. n is the number of grid points per site. The quantity S_c(k,t) is calculated by averaging S(k⃑,t) over the discs {k⃑:|k⃑| ∈ [k,k + 2π/L_b]}, and 〈·〉 denotes the mean over the grid. The maximum and the first moment of the structure factor,are usually used to quantify a characteristic inverse length scale. We define the polymer aggregates at transition time as ‘particles’ and estimate their radius with

l₁ := γ2π/k₁(t_tr) (11)

and

l_max := γ2π/k_max(t_tr) (12)

where k_max(t_tr) := arg max_kS(k,t_tr) and γ := 1/4. We use two estimators since both k₁ and k_max are reasonable choices to quantify a characteristic inverse length scale and we want to assess the difference between the two. The particle radius is thus taken to be one fourth of the characteristic wave length 2π/k₁(t_tr) or 2π/k_max(t_tr), respectively.

In addition, we used a standard image labeling algorithm to determine the total particle number n_p and for each particle we calculate its sphere equivalent radius by

with V and A being the volume or the area of particle i. As a measure for the mean particle size we use the mean radiusIt should be emphasized that we have to restrict to the very early stages of particle formation, where no sharp interfaces are present, if we use the Cahn–Hilliard equation. So actually we are interested in structures as they appear for example at t = 200 in Fig. 2(d) but since it is hard to define a clear measurement specification in the early stages we pick the particles at transition time as representatives because the structures from earlier times seem to imprint onto them.

4.3. Dependence of particle size and transition time on solvent mixing rate

Fig. 3 shows the simulation results for R, l₁, l_max and the transition time t_tr as a function of the mixing rate s. l₁, l_max, and R take slightly different values but the progression of their data points is the same. All simulation results decrease monotonically in s and show the same asymptotic behavior for large s.

Fig. 3 Particle size (a–c) and transition time with t_max (d) vs. solvent mixing rate s in double logarithmic representation for 2D. Error bars in (b–d) indicate the standard deviation over 5 simulation runs. In (a) the error bars indicate the mean standard deviation of droplet radii within a specific simulation run to indicate polydispersity. The statistical variance of R is similar to l₁ and l_max and is not shown. The dashed horizontal lines give simulation results for χ(t) = χ_max and the grey solid lines are regression lines (discussed in Section 4.3). (e–h) Show the same as (a–d) but for 3 dimensions.

We begin with discussing the asymptotic behavior. Since eqn (7) indicates that the value of the time dependent interaction parameter in the limit of infinitely fast solvent mixing is given by

the asymptotes are expected to correspond to the simulation results for an instantaneous quench with constant interaction parameter χ_max. To check this assumption the corresponding results are represented by the dashed horizontal lines in Fig. 3. The data points clearly converge to these lines, hence the asymptotes are consistent with the expectations. In our set of eqn (7)–(9), a constant interaction parameter is achieved by substituting eqn (7) by χ(t) = χ_max. Next we define an asymptotic regime and discuss which values of s belong to that regime. The growth of χ(t) is cut off when t becomes greater than t_max. We call the asymptotic regime the values of s for which the choice of the cutoff affects the simulation results at the transition time. This is clearly the case if t_tr exceeds t_max. Hence the conditiondefines the asymptotic regime. The function t_max(s) is represented by the dash-dotted line in Fig. 3(d) and (h). It crosses the simulation results for t_tr at s ≈ 2 × 10⁻⁴–5 × 10⁻⁴ in both 2D and 3D. The complement of the asymptotic regime is called non-asymptotic regime.

The numerical results for particle size in Fig. 3 show a remarkable similarity to the behavior of structure sizes that occur during continuous cooling of an alloy, which was investigated with a perturbation theory long time ago, including a typical scaling law with an exponent −1/6.³² This similarity does not come as a surprise due to the formal relation of the underlying models. We are going to verify the scaling laws in a semi analytical manner and check if the actual values of the data points agree with this theory and not only their qualitative progression. To this end, we first expand eqn (8) in u about the homogeneous ground state with a mean polymer volume fraction ū. After a Fourier transformation in space, we obtain the ordinary differential equations

where c_m⃑ are the Fourier coefficients of the perturbation (u(r⃑,t) − ū), and a(k,t) is given by
with and m_j ∈ [Doublestruck Z]. The solution to eqn (15) reads and we define
Inserting eqn (7) and (9) into a(k,t) givesWe use k* to estimate the particle radius within the linearized theory. The corresponding estimator is
which is defined analogously to eqn (11) and (12). The comparison between the numerical results and l* is depicted in Fig. 4. The data points are identical to the data points in Fig. 3(a)–(c) but they are plotted in a different representation, namely versus Δχ* := k*² instead of s. In this representation, l* becomeswhich is shown as a straight line in Fig. 4(a). The data points in the asymptotic regime in Fig. 3 collapse onto a single accumulation point at because while t_tr has a lower bound greater than zero (cf.Fig. 3(d) and (h)). Fig. 4 shows that the prediction for the particle size from the linearized theory, l*, approximates the numerical results with a relative deviation of less than 20%. Regarding our interpretation of l*, l_max should be the best approximation and it can be seen from Fig. 4(b) that its deviation is even less than 10%. Hence, we conclude that the perturbation theory yields a good approximation to the numerical results at transition time. The scaling l* ∝ Δχ*^−0.5 also reminds of the relation between particle size and quench depth for constant interaction parameters,^33,34 which gives Δχ* the interpretation of an effective constant quench depth.

Fig. 4 (a) Characteristic particle size vs. parameter Δχ*. The symbols correspond to simulation data shown in Fig. 3(a)–(c), the line represents the prediction of the linear approximation l*. (b) Relative deviation between simulation data and l*. The data for 3D is not shown but looks very similar.

To establish a relation between k* and s we plot t_tr against Δχ* in Fig. 5 and observe that t_tr ∝ Δχ*⁻², which is also reminiscent of a scaling behavior for constant quench depths. The proportionality can be used to formulate approximate scaling laws for the non-asymptotic regime in Fig. 3. Inserting t_tr ∝ Δχ*⁻² = k*⁻⁴ into the case for t_tr < t_max from eqn (16) yields

t_tr ∝ s^−⅔ (18)

for t_tr < t_max. Employing the proportionality (18) into eqn (16) and combining it with definition of l* leads to

l* ∝ s^−⅙ (19)

for t_tr < t_max or s < max{s: t_tr < t_max}. The solid lines in Fig. 3(a)–(c) and (e)–(g) are regression lines to the corresponding data points. Their equations are shown in the diagrams and their exponents deviate about 10% and less from −1/6. Hence, the semi analytical approach verifies the predictions from the perturbation theory. The deviation of the exponents in Fig. 3(d) and (h) from −2/3 in eqn (18) is 3.5% and less. Therefore, the scaling behavior of the numerical data in the non-asymptotic regime comes very close to the predicted scaling behavior from eqn (18) and (19). These two equations are independent from ū and N. The parameter χ_max affects t_max and thus the extent of the non-asymptotic regime, but not the simulation results within that regime. Since the only independent parameters in the model other than s are ū, N, χ_max, the scaling laws in the non-asymptotic regime seem to be an universal feature – at least provided that different choices of ū and N do not destroy the analogy between linearly time dependent and constant interaction parameters. Even though they are not shown in the current publication we performed simulation runs for different parameters and the scaling laws were always observed with the same exponents within an error of 20% and less. Fig. 4 implies that even the actual values of the particles sizes correspond to the perturbation theory.

Fig. 5 Transition times t_tr from Fig. 3(d) plotted vs. Δχ*. The solid line is a regression line. The simulation results for 3D show a very similar scaling.

5. Reference to experiments

From a practical point of view, the most relevant part in Section 4 is the scaling law l* ∝ s^−⅙. Batch experiments with drop injection of selective solvent⁸ report that the mean vesicle or micelle radius depends on the rate of co-solvent addition according to a power law with an exponent of approximately −0.13. The drop-wise co-solvent addition at a constant rate could imply the applicability of a linear time dependence of the interaction parameter allowing a direct comparison between −0.13 and −1/6, which is a good agreement. In experiments where nanoparticles are produced continuously inside micro mixers⁹ it was also observed that the mean particle radius depends on the flow rate according to a power law with an exponent of −0.11, −0.13, or −0.17 depending on the mixer. For a comparison with the micro mixer approach, however, s has to be translated into a flow rate ν. Usually, the mixing time (corresponding to t_max in our model) is inversely proportional to the Reynolds number and thus, to ν.¹⁰ So linear interpolations of the temporal co-solvent volume fraction evolutions in such a mixer show slopes proportional to ν. This leads to scaling laws l* ∝ ν^−⅙, which is also in good agreement with the experiments.

To make a more quantitative comparison we calculate mixing times τ for different flow rates in the Cater Pillar Micro Mixer with an analytical approach³⁵ and assume s = (χ_max − χ₀)t₀/τ, with the time scale t₀ = R_g²/D defined in the beginning. This leads to

where R_g² and D have SI units and ν is given in ml min⁻¹ like in the experiments. The polymer PB₁₃₀PEO₆₆ possesses a molar mass of M ≈ 10 kg mol⁻¹.⁹ Unfortunately, the density for the copolymer was not measured but the homopolymer densities are ρ_PB = 0.96 kg l⁻¹ and ρ_PEO = 1.2 kg l⁻¹, so we estimated the copolymer density by their mean value, ρ ≈ 1.08 kg l⁻¹. The polymer content in the dilute solution was about c = 4 (g polymer) (l solvent)⁻¹. Basic algebra leads to a mean polymer volume fraction of ū = α/(1 + α) with α = c/ρ. Using the values above we have ū = 0.004. Both the molar mass and the density of THF is comparable to the molar mass and the density of the monomers PB and PEO, resulting in similar molar volumes. Thus we estimated N = 190 and set D to the diffusion coefficient of THF in water, which is about 10⁻⁹ m² s⁻¹. We substituted the PEO part by PB and estimated R_g of the resulting homopolymer from its molar mass M by a relation³⁶ which is valid for PB in THF and gives R_g ≈ 10 nm.

Simulations were performed with ū = 0.004, N = 190 and χ_max = 16. It should be noted that we Taylor expanded ln(u) in eqn (9) up to 10th order around ū/10 to avoid numerical difficulties caused by large N. Using R_g and D as mentioned above we converted R to the nanometer scale and calculated flow rates with eqn (20). The results are shown in Fig. 6. The open symbols are data from the experiments for symmetric flow conditions and the black dots represent our simulation results. CPMM, SIMM, SFIMM denote specific types of micro mixers and A and B refer to different end groups attached to the polymer. The SFIMM and SIMM^37–39 are pictured for the sake of completeness – strictly speaking ν is the corresponding flow rate in the CPMM. It can be seen that the model is able to reproduce both the scaling law and typical length and time scales of the experiments but it predicts roughly two times smaller particles. This could either be due to the rather rough approximations for D and R_g, the application of an implicit solvent model,⁴⁰ or the restriction to homopolymers. The final particle size is also influenced by ‘technical’ issues like the choice of u_th, so strict quantitative comparisons should be taken with care. It also should be emphasized that Fig. 6 shows simulation results for homopolymers and experimental data for copolymers, i.e. components of very different composition. Comparing the experimental data for PB₁₃₀PEO₆₆–H (sample A) and PB₁₃₀PEO₆₆–CO–CH₂–CH₂–COOH (sample B) in the CPMM it can be seen that the composition of the polymer chain significantly shifts the data.

Fig. 6 Comparison between our model and experiments from Thiermann et al.⁹ CPMM, SIMM, SFIMM denote mixer types and A and B refer to different end groups of the polymer (see text for details).

6. Summary and outlook

We have described nanoparticle precipitation by spinodal decomposition.

The simulations reproduce power laws as well as typical length scales for the size of vesicles and micelles from experiments.^8,9 These scaling laws are also in par with analytical investigations of spinodal decomposition during continuous cooling³² and our results also agree with more complex particle models for homopolymer precipitation,²⁰ where similar exponents were observed (≈−0.17). Thus, the main result of the present article is that the thermodynamic notion of spinodal decomposition is a promising frame to study size controlled flash nanoprecipitation. Compared to particle models, field theories and especially phase field models require less computation time and grant access to time scales corresponding to mixing times in experiments. Computation time also benefits from the fact that the scaling laws can be investigated in 2D, since 2D and 3D simulations show the same behavior, which allows relatively efficient explorations of parameter spaces. Due to their simple structure even an analytical treatment in the frame of a perturbation theory might be possible.

Scaling laws l ∝ α^−⅙ were also found in a recent publication, which considered the structuring of polymer solutions in the spinodal area upon solvent evaporation,⁴¹ where l is a typical structure size and α a constant evaporation rate. The authors added α as a source term in a Cahn–Hilliard–Cook equation. Within a typical “Flory–Huggins”-phase diagram with axes ū and χ, they advance into the spinodal area in the ū-direction, while we move in the χ-direction. The fact that both processes yield the same scaling behavior suggests that the scaling should just depend on the distance to the spinodal line independent of the direction in the ū–χ-plane.

As far as the comparison between homopolymers and copolymers in Fig. 6 is concerned, a possible interpretation of the similar particle size behavior could be that the co-solvent addition controls the size of the vesicles mainly by determining the size of their micellar predecessors (cf. mechanisms I and II^17,18) and that ‘population balance effects’ like flow induced collision-coagulation and break-up of aggregates in the micro channels might play a minor role, if any. Thus we have also identified one possible mechanism that determines the nanoparticle size in micromixers.

The similar behavior of homopolymer and copolymer particle size might also imply that the principal effect behind size controlled nanoparticle precipitation could be independent of the actual polymer architecture.

In the future, we plan to couple solvent mixing to more sophisticated free energy models,^19,42 which are able to describe copolymers and the vesicle formation process, in order to capture the nanoparticle self-assembly also in the later stages of the aggregation process. Furthermore, it would be interesting to compare simulations for three component systems to our effective two component system and to analyze explicitly how the phase separation process depends on the time-dependent solvent composition. In our study, we have focused on liquid–liquid phase separation, where crystallization and solidification effects can be neglected. Recent experiments on semi-crystalline copolymers⁴³ have shown that the effect of solvent exchange (in this case, solvent evaporation) on the dynamics is very different if demixing interferes with solidification. For example, the characteristic length scales of the resulting structures no longer depend on the solvent evaporation rate, and the experiments can be described within a model based on homogeneous nucleation theory. In the future, it will also be interesting to consider the competition of liquid–liquid phase separation and solidification in more detail.

Acknowledgements

The author acknowledged funding from EFRE.

References

1. L. Zhang, F. X. Gu, J. M. Chan, A. Z. Wang, R. S. Langer and O. C. Farokhzad, Clin. Parmacol. Ther., 2008, 83, 761–769.
2. M. E. Gindy and R. K. Prud'homme, Expert Opin. Drug Delivery, 2009, 6, 865–878.
3. D. E. Discher and F. Ahmed, Annu. Rev. Biomed. Eng., 2006, 8, 323–341.
4. S. Tenzer, D. Docter, S. Rosfa, A. Wlodarski, J. Kuharev, A. Rekik, S. K. Knauer, C. Bantz, T. Nawroth, C. Bier, J. Sirirattanapan, W. Mann, L. Treuel, R. Zellner, M. Maskos, H. Schild and R. H. Stauber, ACS Nano, 2011, 5, 7155–7167.
5. H. Maeda, J. Wu, T. Sawa, Y. Matsumura and K. Hori, J. Controlled Release, 2000, 65, 271–284.
6. D. E. Discher and A. Eisenberg, Science, 2002, 297, 967–973.
7. C. Zhang, V. J. Pansare, R. K. Prud'homme and R. D. Priestley, Soft Matter, 2012, 8, 86–93.
8. W. Müller, Hydrophobe und hodrophile Beladung polymerer Vesikel, Dissertation, Johannes Gutenberg University Mainz, German, 2009.
9. R. Thiermann, W. Mueller, A. Montesinos-Castellanos, D. Metzke, P. Löb, V. Hessel and M. Maskos, Polymer, 2012, 53, 2205–2210.
10. L. Falk and J.-M. Commenge, Chem. Eng. Sci., 2010, 65, 405–411.
11. J. W. Cahn and J. E. Hilliard, J. Chem. Phys., 1958, 28, 258–267.
12. J. W. Cahn and J. E. Taylor, Acta. Metall. Mater., 1994, 42, 1045–1063.
13. J. W. Barrett, J. F. Blowey and H. Garcke, ESAIM: Math. Modell. Numer. Anal., 2001, 35, 713–748.
14. M. V. Ariyapadi and E. B. Naumann, J. Polym. Sci., Part B: Polym. Phys., 1990, 28, 2395–2409.
15. P. J. Flory, J. Chem. Phys., 1942, 10, 51–61.
16. M. L. Huggins, Ann. N. Y. Acad. Sci., 1942, 43, 1–32.
17. X. He and F. Schmid, Phys. Rev. Lett., 2008, 100, 137802.
18. T. Uneyama, J. Chem. Phys., 2007, 126, 114902.
19. X. He and F. Schmid, Macromolecules, 2006, 39, 2654–2662.
20. A. Nikoubashman, V. E. Lee, C. Sosa, R. K. Prud'homme, R. D. Priestley and A. Z. Panagiotopoulos, ACS Nano, 2016, 10, 1425–1433.
21. J. R. Spaeth, I. G. Kevrekidis and A. Z. Panagiotopoulos, J. Chem. Phys., 2011, 135, 184903.
22. H. Tanaka, Phys. Rev. Lett., 1993, 71, 3158–3161.
23. D. Zhou, P. Zhang and E. Weinan, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 73, 061801.
24. J. Zhu, L.-Q. Chen, J. Shen and V. Tikare, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 60, 3564–3572.
25. M. Frigo and S. G. Johnson, The Fastest Fourier Transform in the West, MIT, Cambridge, MA, 1997.
26. M. Matsumoto and T. Nishimura, ACM Trans. Model. Comput. Simul., 1998, 8, 3–30.
27. M. A. Kotnis and M. Muthukumar, Macromolecules, 1992, 25, 1716–1724.
28. M. Fialkowski and R. Holyst, J. Chem. Phys., 2004, 120, 5802–5808.
29. M. Fialkowski and R. Holyst, Macromol. Theory Simul., 2008, 17, 263–273.
30. V. Sofonea and K. R. Mecke, EPJ B, 1999, 8, 99–112.
31. H. Mantz, K. Jacobs and K. Mecke, J. Stat. Mech.: Theory Exp., 2008, 12, P12015.
32. E. L. Huston, J. W. Cahn and J. E. Hilliard, Acta Metall., 1966, 14, 1053–1062.
33. C. P. Grant, Commun. Part. Diff. Eq., 1993, 18, 453–490.
34. C. M. Elliot, Internat. Ser. Numer. Math., 1989, 88, 35–71.
35. F. Schönfeld, K. S. Drese, S. Hardt, V. Hessel and C. Hofmann, NSTI-Nanotech, 2004, vol. 1, pp. 378–381.
36. C. I. D. Bica, W. Burchard and R. Stadler, Macromol. Chem. Phys., 1996, 197, 3407–3425.
37. K. S. Drese, Chem. Eng. J., 2004, 101, 403–407.
38. V. Hessel, S. Hardt, H. Löwe and F. Schönfeld, AIChE J., 2003, 49, 566–577.
39. S. Hardt and F. Schönfeld, AIChE J., 2003, 49, 578–584.
40. J. R. Spaeth, I. G. Kevrekidis and A. Z. Panagiotopoulos, J. Chem. Phys., 2011, 134, 164902.
41. C. Schaefer, P. van der Schoot and J. Michels, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 022602.
42. T. Ohta and A. Ito, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 52, 5250–5260.
43. J. J. van Franeker, G. H. L. Heintges, C. Schaefer, G. Portale, W. Li, M. M. Wienk, P. van der Schoot and R. A. J. Janssen, J. Am. Chem. Soc., 2015, 137, 11783–11794.

---