Three eras of micellization

Micellization is the precipitation of lipids from aqueous solution into aggregates with a broad distribution of aggregation number. Three eras of micellization are characterized in a simple kinetic model of Becker-D\"oring type. The model asigns the same constant energy to the $(k-1)$ monomer-monomer bonds in a linear chain of $k$ particles. The number of monomers decreases sharply and many clusters of small size are produced during the first era. During the second era, nucleii are increasing steadily in size until their distribution becomes a self-similar solution of the diffusion equation. Lastly, when the average size of the nucleii becomes comparable to its equilibrium value, a simple mean-field Fokker-Planck equation describes the final era until the equilibrium distribution is reached.


I. INTRODUCTION
Spontaneous self-assembly of small molecular aggregates in aqueous solutions forms association colloids or complex fluids [1]. Depending on their mean aggregation number, molecular volume and critical hydrocarbon chain length, lipids can pack into spherical or cylindrical micelles. The surfaces of these structures are formed by the hydrophilic heads of the monomer molecules, whose hydrophobic tails lie inside the aggregate. Equilibrium thermodynamics shows that rod-like cylindrical aggregates have a polydisperse distribution of sizes (micellization), whereas the sizes of spherical aggregates grow indefinitely (phase segregation) [1]. The latter process is similar to other examples of first order phase transitions [2] such as condensation of liquid droplets from a supersaturated vapor, colloidal crystallization [3] and the segregation by coarsening of binary alloys quenched into the miscibility gap [4,5,6]. Understanding the kinetics of nucleation and growth beyond the determination of the steady-state nucleation rate is a task of great importance and not yet completely accomplished. This is so despite the large literature on nucleation and growth [7], and several attempts at bridging the gap between nucleation and late-stage coarsening theories [8].
In this paper, we study asymptotically a simple discrete model of micellization kinetics of Becker-Döring type [7,8,9,10]. Starting from an initial condition of pure monomers, we expect the system to evolve to the well known polydisperse equilibrium distribution [1].
However, the nonequilibrium evolution is interesting per se and because the methodology employed here may be applicable to the kinetics of phase segregation. We find that the approach to equilibrium occurs in three well defined stages or eras. Starting from the initial state of pure monomers, the number of monomers decreases sharply and many clusters of small size are produced during the first era. During the second era, aggregates are increasing steadily in size until their distribution becomes a self-similar solution of the diffusion equation. Lastly, when the average size of the nucleii becomes comparable to its equilibrium value, a simple mean-field Fokker-Planck equation describes the final era until the equilibrium distribution is reached. Numerical solution of the model confirms all the theoretical predictions.
The rest of the paper is as follows. In Section II, we review the equilibrium properties of self-assembling aggregates and introduce discrete kinetic models of Becker-Döring type to describe them. Depending on the binding energy of the aggregate with k monomers (k cluster), micellization or phase segregation occurs. For rod-like aggregates, the binding energy of a k cluster (relative to isolated monomers in solution) is (k − 1) times the monomer-monomer bond energy, and an equilibrium size distribution exists (micellization).
For spherical aggregates, the binding energy includes a term proportional to the surface area of the aggregate and no equilibrium size distribution exists beyond a critical density. Then aggregates grow indefinitely and phase segregation occurs following the typical nucleation and growth kinetics. Section III presents a numerical simulation of micellization kinetics which clearly revels its three eras. The agenda of the asymptotic analysis is now clear, and is carried out in Section IV. The last section contains our conclusions and suggestions for experiments.

II. THERMODYNAMICS AND KINETIC MODELS
The model presented here is nucleation in a lattice. There are systems, such as proteins aggregating in a cubic phase of lipid bilayers, for which a lattice formulation is physically correct. In this paper, the main reasons for a lattice model are clarity, and the expectation that the dilute limit of the lattice model (in which there are many more binding sites, M, than particles, N) should closely resemble crystallization from a dilute solution. The latter is a classical problem in the kinetic theory of first order phase transitions [2]. We shall now review the equilibrium statistical mechanics of aggregates, distinguishing between micellization and phase segregation, and then introduce the kinetic models we study.

A. Equilibrium size distribution of aggregates
Let us assume that we have p k ≥ 0 clusters with k particles (in short, k clusters) so that Let e k be the energy of a k cluster. The total energy of the lattice system is where we have used the particle conservation (2.1). Except for a constant Ne 1 , the total energy is Now E is the total lattice energy measured with respect to a configuration in which all clusters are monomers, and ε k is the binding energy of the k cluster (notice the sign convention). We will obtain the equilibrium configuration by minimizing the free energy density with respect to the density of k clusters. To calculate the entropy, we proceed as follows. Let  5) and the entropy of the system is k B ln Ω. In the appropriate thermodynamic limit, N → ∞ with fixed densities ρ ≡ N/M (particles) and ρ k ≡ p k /M (k clusters), particle conservation and we can show that the entropy density is where ρ 1 = ρ − ∞ k=2 kρ k and r = 1 − ∞ k=1 ρ k . In the dilute limit, 1 − r = ∞ k=1 ρ k < ∞ k=1 kρ k = ρ ≪ 1, and therefore r ∼ 1, r ln r ∼ − ∞ k=1 ρ k , and Eq. (2.9) becomes (2.10) which corresponds to the Boltzmann counting. The equilibrium density of k clusters (k ≥ 2) can be found by differentiating this equation with respect to ρ k and equating the result to zero. Taking into consideration that ∂ρ 1 /∂ρ k = −k (k ≥ 2), we obtaiñ (the positive sign in the argument of the exponential is due to our definition of the binding energies). Eq. (2.11) can be rewritten as g k as a function of k can be interpreted as the activation energy of nucleation theory. The equilibrium density of monomers can be found by inserting Eq. (2.11) into Eq. (2.6) and solving the resulting self-consistent equation for ρ 1 in terms of the constant density ρ: (2.14) Whether this self-consistent equation has a solution depends on the value of ρ and on the model we adopt for the binding energy of a k cluster. Typical models are as follows. For rod-like aggregates, where αk B T is the monomer-monomer bonding energy [1]. For spherical aggregates,  Inserting Eq. (2.15) in Eq. (2.14) and using ∞ k=1 (2.17) This equation has the unique solution is the average cluster size in equilibrium. Notice that, for ρe α ≫ 1, k ∼ √ ρe α and For spherical aggregates, the self-consistency condition based on the approximation to ε k in Eq. (2.16) is Clearly this series converges provided ρ 1 e α < 1 and it diverges if ρ 1 e α > 1. The critical monomer concentration ρ 1 = e −α is called critical micelle concentration (CMC) [1]. Below CMC eq. (2.20) can be solved for ρ 1 and the aggregates eventually form micelles with an equilibrium size distribution, whereas phase segregation and indefinite aggregate growth results if more monomers are added above the CMC. For k ≫ 1, the free energy (2.13) is . For ϕ > 0, g k increases for small k, it has a maximum at the critical cluster size k c ≈ (σ/ϕ) 3 , and then it decays monotonically as k further increases.

B. Kinetic models
Let us now formulate the kinetic theory of aggregation in these systems. As in the Becker-Döring kinetic theory, we shall assume that a k cluster can grow or decay by capturing or shedding one monomer at a time. Theṅ or finally, Here D + ε k ≡ ε k+1 − ε k = −D + g k + k B T ln(1/ρ 1 ) and j k is the net rate of creation of a k + 1 cluster from a k cluster, given by the mass action law. We have made the detailed balance assumption to relate the kinetic coefficient for monomer aggregation to that of decay of a (k + 1) cluster, d k . Thenρ k given by Eq. scales describing aggregation kinetics range from microseconds to milliseconds [11,12,13] to be solved together with the conservation condition (2.6), namely, ∞ k=1 kρ k = ρ. At t = 0, we assume that ρ k = ρδ k1 . We shall consider the limit ρ ≫ e −α , in which the initial monomer concentration is much larger than the CMC. The parameters ρ and α are not really independent: if we rescale the cluster densities with ρ, so that and define a scaled time, to be solved with initial conditions Lastly, notice that we can straightforwardly derive two global identities from Eqs. (2.27) and (2.28): Here r c is the total density of clusters and, initially, r c (0) = 1.

III. NUMERICAL RESULTS
Numerical solution of the initial value problem given by Eqs.  reaches a maximum and then decreases to a constant value, as can be seen in Fig. 3.1(c). By the end of the initial stage at time τ = 10, the creation of smaller clusters, with 2 ≤ k ≤ 5, has slowed down greatly relative to the initial spurt for times 0 < τ < 2. Furthermore, the number of clusters with more than 5 monomers is negligible. At τ = 10, k ≈ 2.69, much smaller than the equilibrium value k ≈ √ ρe α = ǫ − 1 2 ≈ 46.9. To determine the time scales appropriate for exploring the subsequent kinetics, it is highly instructive to plot the average cluster size k as a function of time, based on the numerical solution. Fig. 3.5 is a log-log plot of k /e as a function of τ . It reveals an initial rapid growth of k to a "plateau value" close to e, roughly located in the interval 10 < τ < 100. In the subsequent growth after the plateau, large clusters with k ≫ 1 eventually appear.   distribution shape over the whole time span 20 < τ < 1.5 × 10 5 is not very great.
The self-similar stage is not the final chapter of the kinetics story either. By τ = 10 6 , the linear dependence of ln( k /e) with lnτ breaks down. In fact, at τ = 10 6 , k ≈ 31.1, which is comparable to the equilibrium value of 46.9 mentioned before. Evidently, there is a final stage of kinetics in which the size distribution asymptotes to its equilibrium form. Fig. 3.3 is the final era of cluster aggregation, continued from Fig. 3.2, in which snapshots of the size distribution are taken at τ increments of 0.2 × 10 6 , from 0.2 × 10 6 to 4 × 10 6 . Convergence to an exponential distribution with k equal to the equilibrium value of 46.9 is clear.

IV. ASYMPTOTIC THEORY OF MICELLIZATION
In this section, we shall interpret the numerical results shown in Section III by using singular perturbation methods; see Ref. 14 for a general description thereof.

A. Initial transient
Initially, r 1 (0) = 1 and there are no multiparticle aggregates. As we have seen in Section III, the numerical solution of the complete model shows that there is an initial transient stage during which dimers, trimers, etc. form at the expense of the monomers, and that r k ≈ 0 for sufficiently large k. Taking the ǫ → 0 limit of Eqs. (2.30) and (2.31) yields the following planar dynamical system: with ǫ = 0 becomes d(r k e s )/ds = r k−1 e s , which can be solved recursively to yield As τ → ∞, r k → (k − 1)e −1 /k!. Since r 6 (1) = 0.00255, after the initial transient stage there are insignificant numbers of aggregates with more than 5 monomers. In fact, the average aggregate cluster size is k = 1/r c = e, whereas at equilibrium, k ∼ √ ρe α ≫ 1. We therefore conclude that there must be successive transients on time scales much larger than t = O(ǫ).

B. Intermediate transient
Examination of the exact equation (2.27) shows that when r 1 decreases to size O(ǫ) but r 2 , r 3 , . . . are of order 1, all terms in its right hand side are O(ǫ). This suggests rescaling r 1 = ǫR 1 , so that ρ 1 = e −α R 1 , and using the original time t = ǫτ . Eq. (2.27) becomes The global identities (2.30) and (2.31) become where now r c = ǫR 1 + ∞ k=2 r k ∼ ∞ k=2 r k , as ǫ → 0. In the limit ǫ → 0, R 1 − 1 = r 2 /r c and Eq. (4.8) becomes This is a closed system of equations for r 2 , r 3 , . . . , to be solved with the asymptotic values r k = (k − 1)e −1 /k! as initial conditions. It can be shown that the reduced versions of Eqs. (4.10),ṙ c = −(R 1 − 1)r c , and the conservation condition, ∞ k=2 kr k = 1, are upheld automatically by the solution of Eqs. (4.11), so that they are redundant for this stage.
The numerical solution of the reduced system of equations (4.11) for r k , k ≥ 2 closely approximates that of the full system of kinetic equations at this stage. It can be seen that more and more r k become different from zero as t increases and that r k − r k−1 becomes small. This strongly suggests that r k can be approximated by a continuum limit for long times. To find the continuum limit, we set A similar calculation for the total number of clusters yields r c ∼ δ ∞ 0 r(x, T ) dx, which suggests the definition (4.14) We now substitute Eq. (4.12) in Eq. (4.11) and use (4.14) instead of r c . The result is The right hand side of this expression is of order O(δ 2 ), so that the following distinguished limit is obtained if we set b = 2 and take δ → 0: For k = 2, Eq. (4.11) and the scaling (4.12) with a = b = 2 imply that r(0, T ) = 0. Therefore The numerical prefactor is chosen so that particle conservation, given by Eq. (4.13), holds.
It follows from Eq. (4.14) that R c = (πT ) − 1 2 . Hence the average aggregate size is In terms of the original variables k, t and r k , the previous expressions are as t → ∞. These two equations yield which resembles the behavior of the numerical solution of the full kinetic model as indicated in Fig. 3.4. Notice that the average cluster size k corresponding to the solution of Eqs.

C. Equilibrium transient
The large time limit of Eq. (4.19) does not match the equilibrium size distribution, which is r k ∼ ǫe −k √ ǫ in the same scaled units; see Section II. Thus the limit given by Eq. This is the same scaling as in Eq. (4.12) with a = b = 2 and δ = √ ǫ, and therefore we use here the same notation for the variables. With this scaling, the scaled particle conservation Similarly, The scaled version of the global identity (2.32) is Here r 1 = ǫR 1 = ǫ r(ǫ 1 2 , T ). It follows from Eq. (4.25) that The scaled kinetic equation (2.27) is .
We now substitute Eq. (4.26) in this expression, divide it by ǫ 3 and take the limit ǫ → 0.
The result is

Matching with the intermediate transient stage
We represent r(x, T ) as With prefactor 1/T , the particle conservation equation (2.6) and the total cluster density adopt the invariant forms . (4.33) Inserting this equation together with Eq. (4.30) in Eq. (4.28), we obtain Asymptotic similarity as T → 0 means that h(ζ, T ) in Eq. (4.30) has a limit H(ζ) as T → 0.

Trend towards equilibrium
The stationary solution of Eq. (4.28) with the condition (4.29) is r e = e −x x , and the particle conservation condition gives x 2 = 1, so that x = 1. Then the stationary solution of Eq. (4.28) is r e = e −x , which is the sought equilibrium solution. To show that r(x, T ) → r e (x) as T → ∞, we define the associated free energy If we now substitute Eq. (4.28), integrate by parts, and use r(0, T ) = r 0 (0) = 1 and ∞ 0 r dx = 1/ x , we obtain The right hand side of this equation is less or equal than zero because of the Cauchy-Schwarz inequality: Therefore, we have proven that the free energy is a Lyapunov functional. We can rewrite Eq. (4.38) in an equivalent form by definingr 0 = exp[−x x ], and using the identities:   k (t), which solves the approximate system of Equations (4.11) and r c = ∞ k=2 r k with the initial conditions r k (0) = (k−1)e −1 /k!, and (iii) r(x, T ), which solves the nonlinear Fokker-Planck equation (4.28) with the condition