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

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\sim C t^{\alpha}$. To understand the structure-property relationship, the knowledge of the coarsening exponent $\alpha$ 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 mobility term and show that the coarsening for off-critical mixtures is slower than the expected $t^{1/3}$-growth. Instead, we find $\alpha$ to be dependent on the mixture composition and thus from the morphology. Finally, we propose a model to describe the complete coarsening kinetics including the rate constant $C$.


I. 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 0 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][2][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 4 and later Lifshitz, Slyozov and Wagner 5,6 formalized a theory (commonly known as LSWtheory) for 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 Here, L(t 0 ) = L 0 is the characteristic length of the phase separated system at t 0 , α will be called in what follows the "coarsening exponent" and C the "coarsening rate constant". a) Electronic mail: j.harting@fz-juelich.de 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][12][13][14][15][16] , analytical models 1,17,18 and still an active area of research [19][20][21][22][23][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 Eq. (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][26][27][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 31 .
A widely used model for phase separation due to SD and subsequent coarsening is the Cahn-Hilliard equation 41,42 . 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 de-   Table I: Dependencies of the coarsening exponent α using a symmetric phase diagram and a constant mobility Λliterature excerpt noted by ϕ 0 . For polymeric solutions, the Gibbs free energy of mixing is often expressed using the Flory-Huggins equation 43 . 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 = 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 ϕ 0 , 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 = 0.5) constitutes the critical mixture. 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 I 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 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, Refs. 14 and 15 report a ϕ 0 -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 ϕ 0 -dependent and found 1/3 (25/75 and 50/50) and 1/4 (75/25). For this mobility and ϕ 0 = 0.5, Andrews et al. 22 report a time-dependent coarsening rate constant C = C(t).
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 44-47 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,48 .
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 (ϕ 0 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 kinetics of SD. Section II describes the used phase field framework in detail. In Sec. III, we investigate the early stages of SD and check the validity of the well-known scaling laws for the time until demixing t 0 , and the corresponding initial average domain size of the phase separated structure L 0 , over a broad range of parameters. In Sec. IV, we investigate the coarsening behavior during the late stages of the phase separation, focusing on a symmetric material system (N 1 = N 2 = 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 Sec. V.

II. 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 49,50 .
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 43 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 0 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 i v 0 , while χ is the Flory-Huggins interaction parameter. The first term in Eq. (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 an 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 1 = N 2 = 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 1 = 100, N 2 = 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 SD. The equation of the spinodal line χ s is given by 51 .
The composition of the phases after demixing are given by the binodal line (blue dashed line). In the region between the binodal and the spinodal line, the system is metastable and can demix by nucleation and growth of isolated droplets. The binodal concentrations can be computed by the equality of the chemical potential in both phases for both materials (common tangent construction 52 ). 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 41,42 : 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 Sec. I. 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" 53 and the "fast-mode theory" 54 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 In Eqns. (8) and (9), D 1 and D 2 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 55 and the one proposed by Liu, Bardow and Vlugt (LBV) 55 . For both of these models, the compositiondependent diffusion coefficients can be calculated with the help of the self-diffusion coefficients in the pure materials, with D i j 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 11 , D 12 , D 21 , D 22 ) need to be supplied and D 1 and D 2 are given as D 1 (ϕ) = f (D 11 , D 12 , ϕ) and D 2 (ϕ) = f (D 22 , D 21 , ϕ). For the binary case, the LBVassumption reads while the Vignes-model is given by The same holds true for D 2 with D 21 and D 22 . Furthermore, a linear weight of D 11 and D 12 can also be used, 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 1 = N 2 = 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 equation (6) is written in the split form 56 , 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 57-59 with an iterative solver (GMRES method together with a bloc Jacobi preconditioner 60 ). 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 61 , 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. 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.  From the linear analysis of the Cahn-Hilliard equation the time needed for the SD to take place t 0 , and the initial characteristic size L 0 of the separated phases can be calculated as 23,62

III. EARLY STAGE BEHAVIOR
and  Here, a is an arbitrary constant which depends on the definition taken for t 0 . In this work, t 0 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 0 and L 0 given by Eqns. (13) and (14) are compared to numerous simulations were we vary the parameters over a broad range accord-ing to Tab. II. Thereby, we also use sophisticated mobility functions given by Equation (8) and 9 in combination with composition-dependent self-diffusion coefficients (Eqns. (10) to (12)) which themselves differ by decades. The results, shown in Figs. 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 0 and L 0 are required in order to determine L(t) according to Eq. (1). In addition, the form of our equations for t 0 and L 0 allows to directly identify the impact of parameter variations. As illustrated in Figs. 2 and 3, t 0 and L 0 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.   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. We characterize the coarsening kinetics by fitting the characteristic length scale L(t) obtained from the simulations according to Eq. (15), with the growth law, Eq. (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.

IV. LATE STAGE COARSENING KINETICS
The morphological evolution for ϕ 0 = 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 −5 s), while the snapshots in the second and third column are at two respectively three decades later.
For the critical mixture, ϕ 0 = 0.5, we observe the typical highly interconnected, co-continuous morphology (Fig. 5(j)-(l)) which persists until the end of the simulation. For both off-critical compositions ϕ 0 = 0.25 and 0.4 (Figs. 5(a)-(c) and (d)-(f)), the main characteristic is a droplet-in-matrix structure. For ϕ 0 = 0.25, all particles or droplets are close to spherical during the whole coarsening process. In contrast, for ϕ 0 = 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 = 0.45-mixture in the first coarsening stage (Fig. 5(g)) contains mainly co-continuous or worm-like domains. During further coarsening, the co-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. Figure 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 = 0.25, 0.4 and 0.5, the L increases according to power laws up to the latest simulated times so that the data can be nicely fitted with Eq. (1). The obtained value for the coarsening exponent α is slightly sensitive to the time domain selected for the fit and 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, but also perform 10 successive simulations of the same system (ϕ 0 = 0.3, χ = 4). Altogether, we estimate the standard deviation to be roughly s = 0.0075.
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 Sec. I in Tab. I. We note that the typical co-continuous morphology is far away from the LSWassumption but we confirm the LSW-exponent of α = 1/3 for the critical composition. For the very off-critical system ϕ 0 = 0.25, we find the coarsening exponent to be α = 0.3. For the slightly off-critical composition ϕ 0 = 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. 37 , by Rogers and Desai 38 and by Brown and Chakrabarti 39 , although newer studies report α = 1/3 for the off-critical compositions ϕ 0 = 0.25 and 0.3 14,15 .
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 co-continuous morphology (ϕ 0 = 0.5, α = 1/3) while the coarsening is slightly slower for the spherical droplet structure (ϕ 0 = 0.25, α = 0.3). The ellipsoidalshaped droplet structure produces the slowest coarsening kinetics (ϕ 0 = 0.4, α=0.26). This hypothesis is supported by the more complex coarsening kinetics of the ϕ 0 = 0.45 blend: for ϕ 0 = 0.45, a good fit to Eq. (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 −4 s) the morphology is almost co-continuous or showing wormlike 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 co-continuous morphology (ϕ 0 = 0.5, Figs. 5(j)-(l)). At the end of the coarsening (from approximately 1 × 10 −4 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 = 0.4 blend (Fig. 5(f)) for which we obtain α = 0.264.
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 ϕ 0 varying from 0.1 to 0.4, and for ϕ 0 = 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 ϕ 0 . We find that the exponent reaches its maximum value α = 1/3 for ϕ 0 = 0.5, and its minimum value for ϕ 0 = 0.4 whereby α ≈ 0.26. For ϕ 0 smaller than 0.4, α increases again, reaching an asymptotic value around 0.3. The co-continuous morphology (ϕ 0 = 0.5) seems to evolve faster than the spherical droplet structure (ϕ 0 = 0.1 to 0.25), while the ellipsoidal-shaped droplet structure (ϕ 0 = 0.3 to 0.4) produces the smallest values of α. For ϕ 0 = 0.45 and for ϕ 0 = 0.475, a fit to Eq. (1) with a single exponent is not possible and we hence exclude these quench-ing conditions from Fig. 7. For these compositions, the morphology transitions from almost co-continuous to ellipsoidal droplets, associated with a significant decrease of the coarsening exponent α.
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 = 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 = 0.25, corresponding to morphologies made of spherical, isolated droplets, the time evolution of the energy therefore nicely matches the LSW-prediction.
In order to fully understand the coarsening kinetics given by Eq. (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 /µ 0 ,Ĝ nonloc = G nonloc /(µ 0 l 2 sc ) for the energies,l = l/l sc for the lengths,t = t/t sc for the times andΛ = Λ/D sc for the mobilities, with the scaling factors defined as µ 0 = RT /v 0 , l sc = κ/(2µ 0 ), and t sc = κ/(2µ 0 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 adimensionalized form, similar to Eq. (1): whereL = L/l sc is the dimensionless characteristic length and t = t/t sc ,t 0 = t 0 /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 con-stantĈ may vary asĈ = c 1 αΛσ ∆ϕ bino , where c is a fitting prefactor,Λ the dimensionless mobilityσ the dimensionless surface tension and ∆ϕ bino the volume fraction difference between the two separated phases: on the one hand, the dependence onΛ, σ 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 c 1 α , which leads to the proposed equation for C.
The third argument is the calculation of the surface tension of a diffuse interface using the Van der Waals equation σ = κ ( dϕ dx ) 2 dx, which leads to the scalingσ = σ / √ κ µ 0 . Using this in the equation proposed above forĈ and inserting it in Eq. (17), we finally obtain for the scaled coarsening rate constant Note that with these hypotheses, the unit of the constant is m 1/α s −1 as expected and that for α = 1/3 the LSW law C ∼ Λσ µ 0 ∆ϕ bino is recovered. In order to compare this equation with the C values obtained from the fit to Eq. (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 The comparison between the coarsening constants obtained from the fit of the time-dependent characteristic length scale on the one side, and Eq. (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 = 0.45 and ϕ 0 = 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 Eq. (18) with a value of c = 150, for all initial volume fractions ϕ 0 and quench depths.
At the end, all the parameters used in Eq. (1) (t 0 , L 0 , 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 Eqns. (13), (14), (19), (18) and the α values from Fig. 7.

V. 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 mobility. 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, 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 propose that the morphology of the phase-separated mixture might be the reason for these variations of the coarsening kinetics, the cocontinuous morphology leading to the highest values of α, while the smallest values are found for morphologies dominated by ellipsoidal droplets. In addition to the data for α, we proposed a model (Eq. (18)) to predict the coarsening rate constant C. All these results finally allow the prediction of the late-stage coarsening kinetics using Eq. (1) for symmetric blends under the assumption of a constant mobility Λ, without the need of additional simulations.
To shed more light on the observed ϕ 0 -dependence of the coarsening exponent α, we plan to quantitatively analyze the obtained morphologies regarding their interfacial shape distribution, mean curvature and other structure metrics according to Refs. 19,22 . This also includes the characterization using Minkowski descriptors as proposed by Manzanarez et al. 16,63 . 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.