Two-dimensional Cahn–Hilliard simulations for coarsening kinetics of spinodal decomposition in binary mixtures

Abstract

The evolution of the microstructure due to spinodal decomposition in phase separated mixtures has a strong impact on the final material properties. In the late stage of coarsening, the system is characterized by the growth of a single characteristic length scale L ∼ Ct^α. To understand the structure–property relationship, the knowledge of the coarsening exponent α and the coarsening rate constant C is mandatory. Since the existing literature is not entirely consistent, we perform phase field simulations based on the Cahn–Hilliard equation. We restrict ourselves to binary mixtures using a symmetric Flory–Huggins free energy and a constant composition-independent mobility term and show that the coarsening for off-critical mixtures is slower than the expected t^1/3-growth. Instead, we find α to be dependent on the mixture composition and associate this with the observed morphologies. Finally, we propose a model to describe the complete coarsening kinetics including the rate constant C.

---

1 Introduction

When a binary mixture like a polymer solution is quenched in the thermodynamically unstable region of its phase diagram, the system undergoes phase separation via spinodal decompostion (SD). The period until the phases separate and reach their equilibrium concentrations is called “early stage” of the demixing. The time needed for phase separation, denoted t₀ in the following, is material specific and can vary over decades. Then, the system further evolves towards the thermodynamic equilibrium by coarsening of the separated phases, whereby the energetic contribution due to the interfacial tension between the separated domains decreases progressively.^1–3 The characteristic length scale L(t) of the system therefore increases over time. Directly after the phase separation, during the so-called “intermediate stage”, the coarsening behavior is more sophisticated and there is no theory for it. The “late stage” of the evolution, at longer times, has been widely investigated and theoretically described. Ostwald⁴ and later Lifshitz, Slyozov and Wagner^5,6 formalized a theory (commonly known as LSW-theory) for diffusional coarsening of spherical precipitates for the limiting case of zero volume fraction. They predicted that the domains grow with time t as L ∼ t^1/3, the larger domains growing at the expense of the smaller ones. In general, such a power law can be written as

L(t)^1/α − L(t₀)^1/α = C(t − t₀). (1)

Here, L(t₀) = L₀ is the characteristic length of the phase separated system at t₀, α will be called in what follows the “coarsening exponent” and C the “coarsening rate constant”.

The properties of a material blend are greatly affected by its morphology. Coarsening of the morphology is therefore expected to have a significant impact on the material properties. Hence, coarsening is a long lasting subject of interest, mostly investigated together with SD in binary systems with experimental methods,^7–10 numerical simulations,^11–16 analytical models^1,17,18 and still an active area of research.^19–24 Being able to predict the actual time-dependent average domain size or characteristic length scale L(t) is crucial for understanding the morphology–property relationship. In particular, as indicated by eqn (1), the exact knowledge of the coarsening exponent α is therefore of vast importance.

A specific example is the formation of the photoactive layer in organic solar cells. Due to the polymeric nature of the donor materials, phase separation via SD and subsequent coarsening is likely to happen and is assumed to play an important role. During processing, both photo-active materials are dissolved in solution and deposited on a substrate. In the drying process and with the associated concentration increase in the solution, the mixture may become immiscible and undergo phase separation. Whether SD is beneficial^25–28 or undesired and a source of low efficiency of the solar cell^29,30 is still unclear. In any case, the size and morphology of the microstructure is crucial for both efficiency and stability of the photoactive layer.³¹

A widely used model for phase separation due to SD and subsequent coarsening is the Cahn–Hilliard equation.^48,49 It is a classical diffusion equation, but where the driving force for the system evolution is the gradient of the exchange chemical potential, which includes not only mixing terms but also surface tension terms in order to handle separated phases.

In a binary blend like the ones investigated in this article, the system is described by the time- and space-dependent volume fraction field φ of one of the two materials, while the other one is given by 1 − φ. The initial volume fraction is denoted by φ₀. For polymeric solutions, the Gibbs free energy of mixing is often expressed using the Flory–Huggins equation.⁵⁰ It allows to take into account different molar volumes, often denoted as degree of polymerization N_i of the i^th-mixture component, as well as different interaction parameters between each pair of materials. Equal degrees of polymerization result in a phase diagram which is symmetric about φ₀ = 0.5 as depicted in Fig. 1(a), denoted as a symmetric blend or symmetric material system. Differences in N_i lead to a strongly asymmetric phase diagram as illustrated in Fig. 1(b). Besides the molar volumes and the initial volume fraction or initial composition of the mixture φ₀, a binary blend is characterized by the interaction parameter χ which describes the degree of incompatibility of the materials. Focusing on the unstable region of the phase diagram, the degree of incompatibility will be called quench depth in this paper. Almost pure phases are characteristic for deep quenches. The mixture composition with the lowest value of χ which will lead to SD is defined as the critical composition, while other mixture compositions are denoted as off-critical. For a symmetric blend, Fig. 1(a), a mixture with equal volume fractions (denoted in this article as a 50/50-mixture or a mixture with initial composition of φ₀ = 0.5) constitutes the critical mixture.

Fig. 1 (a) Symmetric phase diagram for N₁ = N₂ = 1 and (b) asymmetric phase diagram for N₁ = 100 and N₂ = 1. The spinodal line is calculated according to eqn (5).

An essential property in the Cahn–Hilliard model is the mobility, which describes how fast the system reacts to the thermodynamic driving force and is characterized by the mobility coefficient Λ. Beyond the simplest assumption of a composition-independent constant value for Λ, there exist many other composition-dependent mobility models in the literature.^14,15,22

A lot of research has already been done on simulations of the Cahn–Hilliard equation in 2D using different setups. Table 1 is a selection that summarizes to the best of our knowledge the key references for coarsening kinetics using 2D numerical simulations of the Cahn–Hilliard equation under the assumption of a symmetric phase diagram and a constant composition-independent mobility Λ. It can be seen that the coarsening exponent α for the critical composition is found very robust to be 1/3, while the results for α for off-critical mixtures are inconsistent. Additionally, the coarsening kinetics under the assumption of composition-dependent mobility functions have been reported. Using a parabolic or double-degenerate mobility, ref. 14 and 15 report a φ₀-independent α equal to 1/4. For a one-sided mobility, the authors of ref. 14 obtain α = 1/3.3, while in ref. 15α is reported to be φ₀-dependent and found 1/3 (25/75 and 50/50) and 1/4 (75/25). For this mobility and φ₀ = 0.5, Andrews et al.²² report a time-dependent coarsening rate constant C = C(t).

Table 1 Dependencies of the coarsening exponent α using a symmetric phase diagram and a constant mobility Λ – literature excerpt

φ 0	Quench depth	Coarsening exponent α
0.5	Deep	1/3^14,15,22,32–43
0.45	Deep	0.33^44
0.4	Deep	1/4,^43 0.29,^44 0.30^45
0.4	Shallow	0.21^45
0.35	Deep	1/3^46
0.3	Deep	1/3,^14 0.32,^39 0.29^44
0.25	Deep	1/3^15
0.22	Deep	1/3^20
0.21	Deep	1/3^47
0.05	Deep	1/3^47

---

While symmetric or nearly symmetric phase diagrams are often encountered in metallic alloys, polymeric solutions are typically characterized by highly asymmetric phase diagrams. For these systems, some recent articles^51–54 reveal novel details for the microstructure formation in the early stages of SD. However, the literature on coarsening kinetics is scarce in comparison to the symmetric case.^16,55 Besides the presented studies on coarsening due to diffusion which is the topic of this article, we note that also hydrodynamic interactions influence the coarsening kinetics,^36,56 which is seen in binary fluid mixtures^57,58 and polymer solutions.^59–61

Despite the numerous already available theoretical studies listed above, we believe that a systematic and complete overview of the coarsening kinetics is still missing. A systematic investigation should basically differentiate between the effect of quenching conditions (φ₀ and χ) and the effect of the assumption for the mobility term Λ, both for symmetric as well as asymmetric material systems.

In this work, we perform large scale simulations of the Cahn–Hilliard equation in 2D and systematically investigate the coarsening kinetics of SD. Section 2 describes the used phase field framework in detail. In Section 3, we investigate the early stages of SD and check the validity of the well-known scaling laws for the time until demixing t₀, and the corresponding initial average domain size of the phase separated structure L₀, over a broad range of parameters. In Section 4, we investigate the coarsening behavior during the late stages of the phase separation, focusing on a symmetric material system (N₁ = N₂ = 1) and a constant, composition-independent mobility coefficient Λ. There, we investigate how the coarsening kinetics may vary with the quench depth and the blend composition. Finally, we give a short summary of our findings and conclude in Section 5.

2 Simulation model and implementation

The kinetic evolution of the system towards its thermodynamic equilibrium is simulated with the phase-field framework proposed in our previous papers.^{62, 63}

As a starting point, we describe the system with the help of its free energy functional which defines its thermodynamic properties. Without loss of generality, the free energy functional reads

where V is the system volume, G^loc is the local free energy density and G^nonloc the non-local contribution due to the field gradients. The local part of the free energy corresponds to the Flory–Huggins theory of mixing⁵⁰ and an additional purely numerical contribution

In the equation above, φ is the volume fraction of the first material, R is the gas constant, T the temperature, and v₀ the molar volume of the lattice site in the Flory–Huggins theory. The molar volume of the material i is given by v_i = N_iv₀, while χ is the Flory–Huggins interaction parameter. The first term in eqn (3) corresponds to the entropy of mixing, the second one to the enthalpic interactions and the last one is a purely numeric correction which supports numerical stability, especially for high values of Nχ for which the binodal concentrations can be very close to 0 and 1.

The non-local part of the free energy functional describes the contributions of the interfacial tension with

where κ is the surface tension parameter.

The phase diagram can be computed from the local free energy. As an example, the phase diagram for a symmetric material system with N₁ = N₂ = 1 is given in Fig. 1(a), while a phase diagram where one material has a much higher molar volume than the other one (N₁ = 100, N₂ = 1) is shown in Fig. 1(b). The spinodal line (red full line) defines the instability region of the blend, above which it is unstable and undergoes spontaneous phase separation due to SD which occurs simultaneously throughout the material. In the region between the binodal and the spinodal line, the system is metastable and demix via nucleation and growth. This phase separation mechanism starts locally when an energy barrier is overcome, typically by the presence of nuclei. The equation of the spinodal line χ_s is given by⁶⁴

The composition of the phases after demixing are given by the binodal line (blue dashed line). The binodal concentrations can be computed by the equality of the chemical potential in both phases for both materials (common tangent construction⁶⁵). This can be solved analytically in the symmetric case but only numerically in the asymmetric case.

The kinetic equation for the volume fraction φ, which is a conserved quantity, is based on the formalism initially proposed by Cahn and Hilliard:^48,49

This can be understood as a continuity equation, the fluid flux being proportional to the gradient of the exchange chemical potential density μ^ex, which is the driving force for the system evolution. Λ is the mobility coefficient related to the diffusional properties of both materials. The exchange chemical potential μ^ex is calculated from the free energy functional in the following way:

The first term corresponds to the chemical potential, whereas the second contribution takes into account the potential due to concentration gradients. The Cahn–Hilliard equation ensures that the system progressively minimizes its free energy relative to the volume fraction, i.e. it relaxes towards its thermodynamic equilibrium.

The importance of the mobility Λ for the kinetic evolution of the blend, and especially its dependence on the composition variable φ has already been intensively investigated by many authors, as stated in Section 1. In Cahn's original assumption, the flux is simply proportional to the gradient of the generalized exchange chemical potential through a constant mobility coefficient. However, the mobility has to depend on the local mixture composition in order to ensure the incompressibility constraint together with the Gibbs–Duhem relationship. Several theories have been proposed to derive correct expressions for the coupled fluxes in multinary mixtures, among which the “slow mode theory”⁶⁶ and the “fast-mode theory”⁶⁷ are the most successful. Their names come from the fact that the mutual diffusion coefficient derived from the Cahn–Hilliard equation in a binary system is controlled by the slowest component in the “slow-mode theory”, while it is controlled by the fastest component in the “fast-mode theory”. For a binary system, the mobilities for the fast mode theory Λ_FM and for the slow mode theory Λ_SM respectively read

Λ_FM = φ(1 − φ)²N₁D₁ + φ²(1 − φ)N₂D₂, (8)

In eqn (8) and (9), D₁ and D₂ are the self-diffusion coefficients of both materials, which in general also depend themselves on the composition of the mixture. Various phenomenological models exist in the literature in order to describe the composition dependence of the diffusion coefficients, among others the equation initially proposed by Vignes and the one proposed by Liu, Bardow and Vlugt (LBV).⁶⁸ For both of these models, the composition-dependent diffusion coefficients can be calculated with the help of the self-diffusion coefficients in the pure materials, with D_ij being the self-diffusion coefficient of material i in a matrix of 100% material j. Hence, for a binary mixture, only four self-diffusion coefficients (D₁₁, D₁₂, D₂₁, D₂₂) need to be supplied and D₁ and D₂ are given as D₁(φ) = f(D₁₁, D₁₂, φ) and D₂(φ) = f(D₂₂, D₂₁, φ). For the binary case, the LBV-assumption reads

while the Vignes-model is given by

D₁ = D^φ₁₁D^(1−φ)₁₂. (11)

The same holds true for D₂ with D₂₁ and D₂₂. Furthermore, a linear weight of D₁₁ and D₁₂ can also be used,

D₁ = D₁₁φ + D₁₂(1 −φ). (12)

This renders the final expression of the mobility Λ complicated in the general case. However, if we consider the case of a symmetric material system (N₁ = N₂ = 1) and assume that the system can be described by one single, constant composition-independent self-diffusion coefficient D, both the fast-mode and the slow-mode model reduce to the same expression Λ = Dφ(1 − φ) which is often referred to as the “double degenerate” mobility in the literature. This mobility is used as a model for surface diffusion driven coarsening.^14,15 It emphasizes the diffusion in the interface regions between the bulk phases. Nevertheless, one needs to keep in mind that these simplifications are not suitable to simulate polymeric material systems which are characterized by highly dynamic asymmetries.

The Cahn–Hilliard eqn (6) is written in the split form,⁶⁹ discretized with finite volumes, and solved using an implicit backward Euler finite difference scheme. To this end, the local part of the free energy is linearized consistently. The code is parallelized and the linear system is solved using PETSc^70–72 with an iterative solver (GMRES method together with a bloc Jacobi preconditioner⁷³). We make use of the unconditional stability of the Euler backward scheme to achieve large time steps. Hereby, we use an adaptive time stepping strategy similar to the one described by Wodo and Ganapathysubramanian,⁷⁴ based on the number of iterations required for the convergence of the solver. We perform 2D simulations with periodic boundary conditions in both directions. The phase separation is initiated with the help of a small initial perturbation to the homogeneous volume fraction field. Note that we also could have used the Cahn–Hilliard–Cook equation,⁷⁵ but several works showed that thermal noise has an impact on the early stages but no to very little impact on the coarsening process.^{34,35,37,76,77} The grid resolution is adapted to each parameter set in order to discretize the interface between separated phases with at least five grid points. With the chosen parameters, the interface thickness lies typically between 5 nm and 50 nm so that the grid spacing lies between 1 nm and 8 nm. It has been verified that all simulations are numerically converged in both time and space.

3 Early stage behavior

From the linear analysis of the Cahn–Hilliard equation the time needed for the SD to take place t₀, and the initial characteristic size L₀ of the separated phases can be calculated as^{23, 78}andHere, a is an arbitrary constant which depends on the definition taken for t₀. In this work, t₀ is determined as the time required for the first phase to reach the equilibrium binodal composition within an error of 1%. As a measure for the average domain size or the characteristic length scale L(t), we calculate the 2D-structure factor from the Fourier transform of the volume fraction field. We obtain the probability distribution p(q, t) of q-vectors by integration over all directions. Taking the inverse of the mean value of q over this distribution gives the characteristic length scale as

Both theoretical predictions for t₀ and L₀ given by eqn (13) and (14) are compared to numerous simulations where we vary the parameters over a broad range according to Table 2. Thereby, we also use sophisticated mobility functions given by eqn (8) and 9 in combination with composition-dependent self-diffusion coefficients (eqn (10)–(12)) which themselves differ by decades. The results, shown in Fig. 2 and 3, show an excellent agreement with the predictions and validate our numerical implementation. These results are also useful for the evaluation of the late stage coarsening, since t₀ and L₀ are required in order to determine L(t) according to eqn (1). In addition, the form of our equations for t₀ and L₀ allows to directly identify the impact of parameter variations. As illustrated in Fig. 2 and 3, t₀ and L₀ may change by orders of magnitude. For a quantitative simulation of realistic material systems, this emphasizes the importance of the correct choice of the simulation parameters as well as the mobility assumption.

Table 2 Simulation parameters for investigation of the early stage behavior in Section 3

Parameter	Value per unit
T	300 K
ν 0	10^−3 m^3 mol^−1
ρ i	1000 kg m^−3
k	10^−5 J m^−3
N 1	1–100
N 2	1
χ	0.6–10
κ i	10^−8–10^−12 J m^−1
D ij	10^−9–10^−16 m^2 s^−1
Λ	See text

---

Fig. 2 Simulated time until phase separation t₀versus expected time from eqn (13). We determine the prefactor a to be 9.7.

Fig. 3 Simulated initial characteristic size L₀versus expected size from eqn (14).

Table 3 Simulation parameters for investigation of the late stage coarsening kinetics in Section 4

Parameter	Value per unit
T	300 K
ν 0	10^−3 m^3 mol^−1
ρ i	1000 kg m^−3
k	10^−5J m^−3
N 1 = N2	1
χ	See Fig. 4
κ i	2 × 10^−10J m^−1
D	10^−10m^2 s^−1
Λ	10^−10m^2 s^−1 (=D)

---

4 Late stage coarsening kinetics

In this section, we investigate the late stage coarsening kinetics of a binary symmetric material system for various quenching conditions. The mobility Λ is assumed to be constant and independent of φ. The used simulation parameters are compiled in Table 3. We characterize the coarsening kinetics by fitting the characteristic length scale L(t) obtained from the simulations according to eqn (15), with the growth law, eqn (1). It is our aim to identify the coarsening exponent α and the coarsening rate constant C. We perform 2D simulations with a grid size of 1024 × 1024 elements. The large grid size ensures that the late stage coarsening lasts for at least two decades of physical time with a significant number of phase-separated domains. This is necessary for a precise evaluation of α and C. The tested compositions and quench depths are presented in Fig. 4. We vary the composition φ at fixed interaction parameter χ, the interaction parameter χ at fixed composition φ, and both together at fixed quench depth χ − χ_s. Altogether, we investigate 60 different conditions covering a wide range of the unstable region.

The morphological evolution for φ₀ = 0.25, 0.4, 0.45 and 0.5 at constant quench depth χ − χ_s = 2 (asterisk markers in Fig. 4) is shown in Fig. 5. The snapshots represent the volume fraction of the second material, from 0 (blue) to 1 (yellow), at different times after phase separation: the snapshots in the first column are taken approximately after one decade of late-stage coarsening (4 × 10⁻⁵ s), while the snapshots in the second and third column are at two respectively three decades later.

Fig. 4 Tested quenching conditions in the spinodal region for the late stage coarsening kinetics. The results for the asterisk-marked conditions are presented in Fig. 5 and 6.

Fig. 5 Resulting morphologies for the asterisk-marked quench points in Fig. 4 for φ₀ = 0.25 ((a)–(c)), φ₀ = 0.40 ((d)–(f)), φ₀ = 0.45 ((g)–(i)) and φ₀ = 0.5 ((j)–(l)). The 1st column corresponds to 4 × 10⁻⁵ s, the 2nd column to 4 × 10⁻⁴ s and the 3rd one to 4 × 10⁻³ s. For φ₀ = 0.25 to 0.45, the majority concentration of the mixture constitutes the yellow matrix phase.

For the critical mixture, φ₀ = 0.5, we observe the typical highly interconnected, bi-continuous morphology (Fig. 5(j)–(l)) which persists until the end of the simulation. For both off-critical compositions φ₀ = 0.25 and 0.4 (Fig. 5(a)–(c) and (d)–(f)), the main characteristic is a droplet-in-matrix structure. For φ₀ = 0.25, all particles or droplets are close to spherical during the whole coarsening process. In contrast, for φ₀ = 0.4 some of the droplets have a sausage-like shape at the first stage of the coarsening (Fig. 5(d)). Still, at the end of the simulated coarsening, some droplets deviate from sphericity and have slightly ellipsoidal shapes. The microstructure for the φ₀ = 0.45-mixture in the first coarsening stage (Fig. 5(g)) contains mainly bi-continuous or worm-like domains. During further coarsening, the bi-continuous and worm-like structure elements gradually transform into an ellipsoidal-shaped dominated morphology (Fig. 5(i)). Note that even after three decades of late-stage coarsening (last column of Fig. 5), the simulation boxes are still much larger than the domain size, ensuring good statistics for the evaluation of the coarsening kinetics.

Fig. 6 shows the evolution of the characteristic length scale L(t) for the four morphologies presented in Fig. 5 during the late stage of the coarsening. For φ₀ = 0.25, 0.4 and 0.5, L increases according to power laws up to the latest simulated times so that the data can be nicely fitted with eqn (1). The obtained value for the coarsening exponent α is slightly sensitive to the time domain selected for the fit, because the fit is performed over a time period that is not well defined. Moreover, α can slightly differ between two different simulations performed with exactly the same parameters. In order to evaluate the precision on the measurement of α, we not only vary the time domain used for the fit (start and end time), but also perform 10 successive simulations of the same system (φ₀ = 0.3, χ = 4). Altogether, we estimate the standard deviation to be roughly s = 0.0075. On top, we cross-checked the obtained coarsening exponents α shown in Fig. 6 with the effective growth exponent as proposed by Huse³³ and get the same values.

As it can be seen from Fig. 6, the critical mixture coarsens with a coarsening exponent α = 1/3. This result is in line with the literature as stated in Section 1 in Table 1. We note that the typical bi-continuous morphology is far away from the LSW-assumption but we confirm the LSW-exponent of α = 1/3 for the critical composition. For the very off-critical system φ₀ = 0.25, we find the coarsening exponent to be α = 0.3. For the slightly off-critical composition φ₀ = 0.4, the coarsening exponent is as small as α = 0.26, significantly below the initially expected value of 1/3. This evidence for late-stage coarsening occurring slower than t^1/3 for off-critical mixture compositions has also been obtained by Garcke et al.,⁴³ by Rogers and Desai⁴⁴ and by Brown and Chakrabarti,⁴⁵ although newer studies report α = 1/3 for the off-critical compositions φ₀ = 0.25 and 0.3.^14,15

Fig. 6 Characteristic length scales L as a function of time t for φ₀= 0.5 (blue plus signs), φ₀ = 0.45 (ochre crosses), φ₀ = 0.4 (magenta squares) and φ₀ = 0.25 (red points) corresponding to Fig. 5 and power law reference lines t^α. Note that we multiplied L(t) for φ₀ = 0.50, 0.40 and 0.25 by arbitrary constants for a better visualization.

The comparison of the morphologies (Fig. 5) and of the coarsening kinetics (Fig. 6) suggests an influence of the morphology on the coarsening exponent. The highest value of α is found for the bi-continuous morphology (φ₀ = 0.5, α = 1/3) while the coarsening is slightly slower for the spherical droplet structure (φ₀ = 0.25, α = 0.3). The ellipsoidal-shaped droplet structure produces the slowest coarsening kinetics (φ₀ = 0.4, α = 0.26). This hypothesis is supported by the more complex coarsening kinetics of the φ₀ = 0.45 blend: for φ₀ = 0.45, a good fit to eqn (1) is not possible. We find that the evolution of L(t) is characterized by a transition between two different coarsening exponents. At the beginning of the late-stage coarsening (until approximately 1 × 10⁻⁴ s) the morphology is almost bi-continuous or showing worm-like droplets (Fig. 5(g)), and the structure coarsens with α = 0.314, which is very close to the value of 1/3 obtained for the fully bi-continuous morphology (φ₀ = 0.5, Fig. 5(j)–(l)). At the end of the coarsening (from approximately 1 × 10⁻⁴ s to the end of the simulation), the morphology is more and more dominated by ellipsoidal-shaped domains (Fig. 5(i)) and the structure finally coarsens with α = 0.264. This is similar to the φ₀ = 0.4 blend (Fig. 5(f)) for which we obtain α = 0.259.

The complete overview of the coarsening exponents α obtained for the quenching conditions marked in Fig. 4 is given in Fig. 7. In all simulations with φ₀ varying from 0.1 to 0.4, and for φ₀ = 0.5, the morphology coarsens following a power law. On the one hand, at fixed composition, we are not able to observe a significant dependence (within two times the standard deviation) of the coarsening exponent on the quench depth. The red crosses in Fig. 7 are the mean values of the obtained coarsening exponent α for a specific φ for all the tested quench depths. On the other hand, as already discussed, there is a strong dependence of α on φ₀. We find that the exponent reaches its maximum value α = 1/3 for φ₀ = 0.5, and its minimum value for φ₀ = 0.4 whereby α ≈ 0.26. For φ₀ smaller than 0.4, α increases again, reaching an asymptotic value around 0.3. The bi-continuous morphology (φ₀ = 0.5) seems to evolve faster than the spherical droplet structure (φ₀ = 0.1 to 0.25), while the ellipsoidal-shaped droplet structure (φ₀ = 0.3 to 0.4) produces the smallest values of α. For φ₀ = 0.45 and for φ₀ = 0.475, a fit to eqn (1) with a single exponent is not possible and we hence exclude these quenching conditions from Fig. 7. For these compositions, the morphology transitions from almost bi-continuous to ellipsoidal droplets, associated with a significant decrease of the coarsening exponent α.

Fig. 7 Summary of the obtained coarsening exponents α as a function of the initial mixture composition φ₀ for different quench depths. The red stars indicate the mean values of the obtained exponents.

In addition, we also evaluate α based on the free energy decrease of the system: similar to the domain sizes, the energy is expected to decrease following a power law. We therefore evaluate and fit the decrease of the interfacial energy for all simulations reported in Fig. 7 and find again composition-dependent coarsening exponents. Interestingly, the α-values obtained with this method are systematically about 0.02 higher than the α-values obtained from the evaluation using the structure factor. For φ₀ = 0.1, we obtain a mean value for α of 0.32 which is closer to the LSW-prediction of t^1/3 than the value we found using the characteristic length scale method. For compositions up to φ₀ = 0.25, corresponding to morphologies made of spherical, isolated droplets, the time evolution of the energy therefore nicely matches the LSW-prediction. To the best of our knowledge, the question of growth exponents for off-critical compositions has not been resolved, even in recent studies.^14,15,20 These papers claim α = 1/3 for off-critical compositions below φ₀ = 0.3 and they evaluate α using the decay of the interfacial energy. For such compositions and with this evaluation method we find α-values around 0.32. Therefore, we have the feeling that our results are not contradictory to previous results, but provide more precise insights.

In order to fully understand the coarsening kinetics given by eqn (1), we finally need to focus on how the constant C depends on the model parameters. In particular, since the coarsening exponent α does not always have the same value, we also expect C to be α-dependent. In the following, we propose an equation for C and check its validity with the help of three arguments.

The first argument is a scaling argument. In our phase-field framework, the equations are written and solved in a dimensionless form using Ĝ^loc = G^loc/μ₀, Ĝ^nonloc = G^nonloc/(μ₀l_sc²) for the energies, l̂ = l/l_sc for the lengths, t̂ = t/t_sc for the times and [capital Lambda, Greek, circumflex] = Λ/D_sc for the mobilities, with the scaling factors defined as μ₀ = RT/v₀, , and t_sc = κ/(2μ₀D_sc). D_sc has the unit of a diffusion coefficient and, since we consider a simple model with constant mobility, can be chosen for instance as D_sc = Λ. Naturally, the equation for the coarsening kinetics can be written in the non-dimensionalized form, similar to eqn (1):

L̂(t)^1/α − L̂(t₀)^1/α = Ĉ(t̂ − t̂₀), (16)

where L̂ = L/l_sc is the dimensionless characteristic length and t̂ = t/t_sc, t̂₀ = t₀/t_sc are the dimensionless characteristic times. Using the definition of t_sc and l_sc, this leads to

As a second argument, we propose that the adimensional constant Ĉ may vary as , where c is a fitting prefactor, [capital Lambda, Greek, circumflex] the dimensionless mobility, [small sigma, Greek, circumflex] the dimensionless surface tension and Δφ_bino the volume fraction difference between the two separated phases: on the one hand, the dependence on [capital Lambda, Greek, circumflex], [small sigma, Greek, circumflex] and Δφ_bino is simply inspired by the classical models based on the LSW theory. On the other hand, Ĉ has also to depend on α. Following the idea developed in the previous section that the variation of the coarsening exponent α may only be dependent on the morphology of the blend and neither on the mobility, nor on the quench depth and/or surface tension, we assume that the α dependency only applies to the proportionality factor as , which leads to the proposed equation for Ĉ.

The third argument is the calculation of the surface tension of a diffuse interface using the van der Waals equation , which leads to the scaling . Using this in the equation proposed above for Ĉ and inserting it in eqn (17), we finally obtain for the scaled coarsening rate constant

Note that with these hypotheses, the unit of the constant is m^1/α s⁻¹ as expected and that for α = 1/3 the LSW law is recovered.

In order to compare this equation with the C-values obtained from the fit to eqn (1), we not only need the binodal composition which can be picked up from the phase diagram, but also the value of the surface tension. To this end, we calculate the one dimensional equilibrium interface profiles between two separated domains, for various parameter sets, varying all the parameters of the Flory–Huggins equation, and calculate the associated surface tension with the van der Waals formula. We are able to fit the surface tension nicely (see Fig. 8) for any binary system with

Fig. 8 Surface tension obtained from the simulations versus values calculated with eqn (19).

The comparison between the coarsening constants obtained from the fit of the time-dependent characteristic length scale on the one side, and eqn (18) on the other side, is shown in Fig. 9. Once again, the data for which the late stage cannot be fitted with a power law (namely φ₀ = 0.45 and φ₀ = 0.475) are not taken into account, but apart from that, all simulations reported in Fig. 7 are also used for Fig. 9. Despite of some significant deviations, the measured coarsening rate constant can be very nicely predicted with eqn (18) with a value of c = 5.31, for all initial volume fractions φ₀ and quench depths.

Fig. 9 Coarsening rate constants C obtained from the simulations versus values calculated with eqn (18).

At the end, all the parameters used in eqn (1) (t₀, L₀, C, α) are known for systems investigated here (symmetric blends with constant mobilities in two dimensions). This means that the coarsening kinetics can be predicted without the need of any further simulation using eqn (13), (14), (19), (18) and the α values from Fig. 7.

5 Conclusion

In this paper, we simulated the spinodal decomposition and subsequent coarsening of immiscible binary blends in 2D, using conserved Cahn–Hilliard dynamics together with the Flory–Huggins free energy of mixing.

First, we simulated the early stages of SD for a broad range of parameters, including asymmetric blends and highly asymmetric, composition-dependent mobility functions. We could check that the time required for phase separation and the initial characteristic size matches the theoretical predictions.

Second, we investigated the late stage coarsening dynamics of a symmetric blend with a constant composition-independent mobility. In fact, at the beginning of our study, considering early papers and the LSW-theory, we expected to find α = 1/3 for every quench point in the unstable region of the phase diagram. But a closer look at the literature including recent papers showed that the situation is more complex than that. We therefore decided to conduct exhaustive investigations: the quench depth and blend composition were varied systematically. As expected, we found that the growth of the characteristic length scale follows a power law. We found no systematic and clear dependence of α on the quench depth. However, surprisingly, the coarsening exponent α was found to be composition-dependent, starting from α = 0.30 for strongly off-critical mixtures, with a minimum of about 0.26 for 40 : 60 blends and reaching a maximum of 1/3 for the critical mixture.

We hypothesize that the morphology of the phase-separated mixture might be the reason for these variations of the coarsening kinetics: When the asymmetry in composition decreases (i.e. we get closer to 50/50 starting from 10/90), we expect more and more droplets deviating from sphericity and we hypothesize that these changes of the morphology lead to lower α-values, up to φ₀ = 0.4. On the other hand the bi-continuous structure (symmetric blend) seems to be associated with the highest α-value (1/3). For φ₀ values in between (φ₀ = 0.45 and 0.475) the structure is at the beginning almost bi-continuous and accordingly α relatively high (close to the value for the bi-continuous network). But with time the structure breaks down in very aspheric droplets. This is correlated with a strong decrease of α (close to the minimum value measured for φ₀ = 0.4).

In addition to the data for α, we proposed a model (eqn (18)) to predict the coarsening rate constant C. All these results finally allow the prediction of the late-stage coarsening kinetics using eqn (1) for symmetric blends under the assumption of a constant mobility Λ, without the need of additional simulations.

We think that our simulations of large systems combined with long simulated times provide very precise results as compared to previous works. We therefore hope that our results are of valuable interest for a broad audience to re-stimulate discussion across the community. For the future, we definitely see the need to provide a quantitative structural measurement or characterization of the obtained morphologies. This would be highly beneficial to check our hypothesis of a morphology-dependent coarsening exponent α. We would like to emphasize that for now, this hyphothesis only relies on a qualitative observation of the morphologies. Therefore, we plan to quantitatively analyze the obtained morphologies regarding their interfacial shape distribution, mean curvature and other structure metrics according to ref. 19 and 22. This also includes the characterization using Minkowski descriptors as proposed by Manzanarez et al.^16,79 Furthermore, we plan to investigate asymmetric material systems in near future, as well as systems with more sophisticated mobility functions in order to extend, confirm and enrich the results on the coarsening behavior of binary blends.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

The authors acknowledge financial support by the German Research Foundation DFG (project HA4382/14-1). The data that support the findings of this study are available from the corresponding author upon reasonable request.

Notes and references

1. F. Otto, C. Seis and D. Slepcev, Commun. Math. Sci., 2013, 11, 441–464.
2. A. J. Bray, Adv. Phys., 2002, 51, 481–587.
3. L. Ratke and P. W. Voorhees, Growth and Coarsening, Springer, 2002.
4. W. Ostwald, Z. Phys. Chem., 1897, 22, 289–330.
5. I. M. Lifshitz and V. V. Slyozov, J. Phys. Chem. Solids, 1961, 59, 263–268.
6. C. Wagner, Z. Elektrochem., 1961, 65, 581–591.
7. P. D. Graham, A. J. Pervan and A. J. McHugh, Macromolecules, 1997, 30, 1651–1655.
8. C. K. Haas and J. M. Torkelson, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 55, 3191–3201.
9. K. S. McGuire, A. Laxminarayan and D. R. Lloyd, Polymer, 1995, 36, 4951–4960.
10. S. Piçarra, E. J. N. Pereira, E. N. Bodunov and J. M. G. Martinho, Macromolecules, 2002, 35, 6397–6403.
11. D. Sappelt and J. Jäckle, Europhys. Lett., 1997, 37, 13–18.
12. R. Ahluwalia, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 59, 263–268.
13. J. Harting, G. Giupponi and P. V. Coveney, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 041504.
14. G. Sheng, T. Wang, Q. Du, K. G. Wang, Z. K. Liu and L. Q. Chen, Commun. Comput. Phys., 2010, 8, 249–264.
15. S. Dai and Q. Du, J. Comput. Phys., 2016, 310, 85–108.
16. H. Manzanarez, J. Mericq, P. Guenoun, J. Chikina and D. Bouyer, Chem. Eng. Sci., 2017, 173, 411–427.
17. R. V. Kohn and F. Otto, Commun. Math. Phys., 2002, 229, 375–395.
18. A. Novick-Cohen, in Handbook of Differential Equations: Evolutionary Equations, ed. C. M. Dafermos and M. Pokorny, North-Holland, 2008, vol. 4, pp. 201–228.
19. C.-L. Park, J. Gibbs, P. Voorhees and K. Thornton, Acta Mater., 2017, 132, 13–24.
20. Y. Sun, W. B. Andrews, K. Thornton and P. W. Voorhees, Sci. Rep., 2018, 8, 17940.
21. A. Banc, J. Pincemaille, S. Costanzo, E. Chauveau, M.-S. Appavou, M.-H. Morel, P. Menut and L. Ramos, Soft Matter, 2019, 15, 6160–6170.
22. W. B. Andrews, P. W. Voorhees and K. Thornton, Comput. Mater. Sci., 2020, 173, 109418.
23. J. S. Higgins and J. T. Cabral, Macromolecules, 2020, 53, 4137–4140.
24. P. K. Inguva, P. J. Walker, H. W. Yew, K. Zhu, A. J. Haslam and O. K. Matar, Soft Matter, 2021, 17, 5645–5665.
25. S. Kouijzer, J. J. Michels, M. van den Berg, V. S. Gevaerts, M. Turbiez, M. M. Wienk and R. A. J. Janssen, J. Am. Chem. Soc., 2013, 135, 12057–12067.
26. V. Negi, O. Wodo, J. J. van Franeker, R. A. J. Janssen and P. A. Bobbert, ACS Appl. Energy Mater., 2018, 1, 725–735.
27. L. Ye, H. Hu, M. Ghasemi, T. Wang, B. A. Collins, J.-H. Kim, K. Jiang, J. H. Carpenter, H. Li, Z. Li, T. McAfee, J. Zhao, X. Chen, J. L. Y. Lai, T. Ma, J.-L. Bredas, H. Yan and H. Ade, Nat. Mater., 2018, 17, 253–260.
28. L. Ye, S. Li, X. Liu, S. Zhang, M. Ghasemi, Y. Xiong, J. Hou and H. Ade, Joule, 2019, 3, 443–458.
29. N. Li, J. D. Perea, T. Kassar, M. Richter, T. Heumueller, G. J. Matt, Y. Hou, N. S. Güldal, H. Chen, S. Chen, S. Langner, M. Berlinghof, T. Unruh and C. J. Brabec, Nat. Commun., 2017, 8, 14541.
30. A. Wadsworth, Z. Hamid, J. Kosco, N. Gasparini and I. McCulloch, Adv. Mater., 2020, 32, 2001763.
31. A. Armin, W. Li, O. J. Sandberg, Z. Xiao, L. Ding, J. Nelson, D. Neher, K. Vandewal, S. Shoaee, T. Wang, H. Ade, T. Heumüller, C. Brabec and P. Meredith, Adv. Energy Mater., 2021, 11, 2003570.
32. K. Binder and D. Stauffer, Phys. Rev. Lett., 1974, 33, 1006–1009.
33. D. A. Huse, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 34, 7845–7850.
34. S. Puri, Phase Transitions, 1989, 16, 407–411.
35. S. Puri and Y. Oono, J. Phys. A: Math. Gen., 1988, 21, L755–L762.
36. K. Binder and P. Fratzl, Spinodal Decomposition, John Wiley and Sons, Ltd, 2001, ch. 6, pp. 409–480.
37. S. Puri and Y. Oono, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 1542–1565.
38. A. M. Lacasta, A. Hernáandez-Machado, J. M. Sancho and R. Toral, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 45, 5276–5281.
39. A. M. Lacasta, J. M. Sancho, A. Hernández-Machado and R. Toral, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 48, 6854–6857.
40. T. M. Rogers, K. R. Elder and R. C. Desai, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 9638–9649.
41. S. Puri, A. J. Bray and J. L. Lebowitz, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 56, 758–765.
42. J. Zhu, L.-Q. Chen, J. Shen and V. Tikare, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 60, 3564–3572.
43. H. Garcke, B. Niethammer, M. Rumpf and U. Weikard, Acta Mater., 2003, 51, 2823–2830.
44. T. M. Rogers and R. C. Desai, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 39, 11956–11964.
45. G. Brown and A. Chakrabarti, J. Chem. Phys., 1993, 98, 2451–2458.
46. S. Puri, Phys. Lett. A, 1988, 134, 205–210.
47. R. Toral, A. Chakrabarti and J. D. Gunton, Phys. A, 1995, 213, 41–49.
48. J. W. Cahn and J. E. Hilliard, J. Chem. Phys., 1958, 28, 258–267.
49. J. W. Cahn, Acta Metall., 1961, 9, 795–801.
50. P. J. Flory, J. Chem. Phys., 1942, 10, 51–61.
51. Y. Mino, T. Ishigami, Y. Kagawa and H. Matsuyama, J. Membr. Sci., 2015, 483, 104–111.
52. G. Zhang, T. Yang, S. Yang and Y. Wang, Phys. Rev. E, 2017, 96, 032501.
53. F. Wang, P. Altschuh, L. Ratke, H. Zhang, M. Selzer and B. Nestler, Adv. Mater., 2019, 31, 1806733.
54. F. Wang, L. Ratke, H. Zhang, P. Altschuh and B. Nestler, J. Sol-Gel Sci. Technol., 2020, 94, 356–374.
55. B. F. Barton, P. D. Graham and A. J. McHugh, Macromolecules, 1998, 31, 1672–1679.
56. E. D. Siggia, Phys. Rev. A: At., Mol., Opt. Phys., 1979, 20, 595–605.
57. N. González-Segredo, M. Nekovee and P. V. Coveney, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 046304.
58. P. Guenoun, R. Gastaud, F. Perrot and D. Beysens, Phys. Rev. A: At., Mol., Opt. Phys., 1987, 36, 4876–4890.
59. A. E. Bailey, W. C. K. Poon, R. J. Christianson, A. B. Schofield, U. Gasser, V. Prasad, S. Manley, P. N. Segre, L. Cipelletti, W. V. Meyer, M. P. Doherty, S. Sankaran, A. L. Jankovsky, W. L. Shiley, J. P. Bowen, J. C. Eggers, C. Kurta, T. Lorik, P. N. Pusey and D. A. Weitz, Phys. Rev. Lett., 2007, 99, 205701.
60. S.-W. Song and J. M. Torkelson, J. Membr. Sci., 1995, 98, 209–222.
61. F. S. Bates and P. Wiltzius, J. Chem. Phys., 1989, 91, 3258–3274.
62. O. J. J. Ronsin and J. Harting, Energy Technol., 2020, 8, 1901468.
63. O. J. J. Ronsin, D. Jang, H.-J. Egelhaaf, C. J. Brabec and J. Harting, Phys. Chem. Chem. Phys., 2020, 22, 6638–6652.
64. M. Rubinstein and R. H. Colby, Polymer Physics, Oxford University Press, 2003.
65. A. D. Pelton, Phase Diagrams and Thermodynamic Modeling of Solutions, Elsevier, 2019.
66. P. G. de Gennes, J. Chem. Phys., 1980, 72, 4756–4763.
67. E. J. Kramer, P. Green and C. J. Palmstrøm, Polymer, 1984, 25, 473–480.
68. C. Peters, L. Wolff, T. J. H. Vlugt and A. Bardow, Experimental Thermodynamics Volume X: Non-equilibrium Thermodynamics with Applications, The Royal Society of Chemistry, 2016, pp. 78–104.
69. C. M. Elliott, D. A. French and F. A. Milner, Numer. Math., 1989, 54, 575–590.
70. S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, D. Karpeyev, D. Kaushik, M. G. Knepley, D. A. May, L. C. McInnes, R. T. Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang and H. Zhang, PETSc Web page, https://www.mcs.anl.gov/petsc, 2021, https://www.mcs.anl.gov/petsc.
71. S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, D. Karpeyev, D. Kaushik, M. G. Knepley, D. A. May, L. C. McInnes, R. T. Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang and H. Zhang, PETSc Users Manual, Argonne National Laboratory Technical Report ANL-95/11 - Revision 3.15, 2021.
72. S. Balay, W. D. Gropp, L. C. McInnes and B. F. Smith, Modern Software Tools in Scientific Computing, 1997, pp. 163–202.
73. Y. Saad and M. H. Schultz, SIAM J. Sci. Comput., 1986, 7, 856–869.
74. O. Wodo and B. Ganapathysubramanian, Comput. Mater. Sci., 2012, 55, 113–126.
75. H. Cook, Acta Metall., 1970, 18, 297–306.
76. X. Zheng, C. Yang, X.-C. Cai and D. Keyes, J. Comput. Phys., 2015, 285, 55–70.
77. S. Pfeifer, B. S. S. Pokuri, O. Wodo and B. Ganapathysubramanian, in Uncertainty Quantification in Multiscale Materials Modeling, ed. Y. Wang and D. L. McDowell, Woodhead Publishing, 2020, pp. 301–327.
78. J. Cabral and J. Higgins, Prog. Polym. Sci., 2018, 81, 1–21.
79. H. Manzanarez, J. Mericq, P. Guenoun and D. Bouyer, J. Membr. Sci., 2021, 620, 118941.

---