Continuous theory of phase separation:The Cahn-Hilliard equation

Abstract

In this paper, we present a detailed deduction of the Cahn-Hilliard equation. It is a non-linear partial differential equation that continuously describes the dynamics of spinodal decomposition, i.e., spontaneous phase separation, of a binary mixture. From such an equation we can model systems that undergo phase separation toward the equilibrium, i.e., passive systems, and systems that are kept out of the equilibrium by the presence of either external forces or chemical reactions, i.e., active systems. To accomplish our main goal, we present the main concepts of non-equilibrium thermodynamics and from that, we obtain the Cahn-Hilliard equation by defining a general equation for the transport of matter and calculating the free energy for a homogeneous and non-homogeneous mixture of real fluids. By the end, we perform simulations of the Cahn-Hilliard equation and present two examples of passive and active systems that (quasi) suppress Ostwald ripening.

Keywords:
Non-equilibrium thermodynamics; liquid-liquid phase separation; passive and active systems

1. Introduction

Phase separation refers to the transient, non-equilibrium process by which a homogeneous mixture spontaneously separates into a multiphase mixture as a consequence of appropriate changes in thermodynamic variables, such as temperature, chemical concentration, pressure, etc. [¹]. Our daily life provides examples of this phenomenon, e.g., the condensation in the shower due to the supersaturation of the air with gaseous water and the formation of oil droplets in water in a vinaigrette [²–⁴]. These well-demixed mixtures have distinct functions and occur in different environments, ranging from industrial processes, e.g., medicine, food, energy, cosmetics, and agrochemicals; [⁵, ⁶] to nature, e.g., giving rise to colors in birds and forming biomolecular condensates in cells [⁷,⁸,⁹,¹⁰]. This last example, i.e., the biomolecular condensates, represents a novel class of physical systems and, because of that, it has been the main subject of scientific investigation by many physicists over the past decades.

Biomolecular condensates are membrane-less micro-compartments, or organelles, in the interior of eukaryotic and procaryote cells [¹¹, ¹²]. These structures concentrate biomolecules, creating a distinct chemical setting from their surroundings and aiding spatiotemporal regulation of biological reactions [¹, ¹¹–¹⁴,¹⁵]. They participate in diverse intracellular processes, such as RNA metabolism, ribosome biogenesis, DNA damage response, signal transduction, nucleation of microtubules, production of spliceosomes, etc [¹¹]. A well-characterized example of biomolecular condensates is the P granules first observed in the 1980s and found in Caenorhabditis elegans embryos [¹³, ¹⁶, ¹⁷]. Recent studies with P granules have shown they present liquid-like properties, i.e., they exchange material with the cytoplasm, exhibit a well-defined contact angle, present a round morphology that deforms in the presence of flows, and have an analogous viscosity to runny honey [¹, ¹³]. These facts suggest that biomolecular condensates are liquid droplets formed by phase separation [³]. However, differently from oil-in-water droplets, i.e., passive droplets, biomolecular condensates exist in a non-equilibrium environment of living cells, so that their material can chemically react and interact with external forces [³, ¹³]. In such conditions, these droplets present different stability properties than the passive ones, configuring a novel class of physical systems, named active droplets/emulsions [³].

Active emulsions represent a subset of a broad class called active matter. Active matter is a class of natural or artificial systems driven out of the thermodynamic equilibrium by a constant supply of energy at the single particle level, which can convert the energy consumed into a form of mechanical work [¹⁸]. Differently from passive systems, which evolve toward the state that corresponds to the free energy minimum, active systems exhibit exciting new dynamical properties including collective behavior, non-equilibrium order/disorder transitions, pattern formation on mesoscopic scales, unusual mechanical and rheological properties, among others [¹⁹]. In living entities, active behavior emerges in the cytoskeleton biopolymers, bacterial suspensions, and terrestrial, aerial, and aquatic flocks [¹⁹]. On the other hand, in nonliving systems active matter arises in layers of vibrated granular rods, self-propelled particles (Janus particles), arrays of coupled synthetic chemical cells, and reconfigurable, self-propelled sheets [¹⁹,²⁰,²¹,²²]. In particular, biomolecular condensate is an example of active matter originating from the interplay between phase separation and chemical reactions [², ³]. It has been shown that such condensates may perform self-propulsion, suppress Ostwald ripening, and stimulate growth and droplet division [¹⁵, ²³, ²⁴]. These intriguing behaviors can be described through physical models based on non-equilibrium thermodynamics, polymer, and soft matter physics [¹³].

Based on the crucial role of phase separation for the well-functioning of non-living and living systems, well-exemplified by the biomolecular condensates, we present in this paper the physical derivation of the Cahn-Hilliard (CH) equation,

$$\frac{\partial \varphi (r,t)}{\partial t}=\nabla \cdot (M\nabla \overline{\mu }),$$

where ϕ(r, t) is the local volume fraction of a chemical, M is the mobility coefficient, $\overline{\mu }=\delta F[\varphi ]/\delta \varphi$ is the local exchange chemical potential, and F[ϕ] is the free energy functional. The CH equation provides a continuous description of the phase separation for binary mixtures and it is the starting point to address more complex problems, e.g., active emulsions. To accomplish that, we begin introducing in section ² the main concepts of the theory of non-equilibrium thermodynamics. From such a theoretical ground, we obtain in section ³ a general dynamical equation for the transport of matter in an incompressible, isothermal, and non-reactive mixture. To proceed with our main task, we seek to calculate the exchange chemical potential through the determination of the free energy density for non-ideal fluids from an Ising-like Hamiltonian in Section ⁴. In section ⁵, we provide corrections to the free energy defined in section ⁴ to include the energetic effects of the presence of internal interfaces, resulting in a free energy functional. Finally, in section ⁶, we obtain the exchange chemical potential, determine the CH equation, carry out simulations of passive phase separation, and discuss two possible ways to (quasi) suppress Ostwald ripening from a passive and active system.

2. Theoretical Foundation of Non-equilibrium Thermodynamics

The classical thermodynamics (CT), as it was formulated in the nineteenth century, has its foundation in the concept of “equilibrium”. In general, such a theory assumes that a system evolves from one state to another reversibly, as a succession of equilibrium states [²⁵]. From that approach, it is possible to compute changes in thermodynamic quantities but not their evolution during the process [²⁵, ²⁷]. This is because at the idealized equilibrium scenario there is no dissipation, the physical properties are time-invariant, and homogeneous throughout the system [²⁵, ²⁸]. However, most of the natural phenomena take place out of equilibrium, i.e., the situation where the system dissipates energy and inhomogeneities induce internal flows [²⁹]. This fact required an extension of the CT to handle more realistic problems, resulting in the Non-equilibrium Thermodynamics theory.

Non-equilibrium Thermodynamics (NT) aims to describe irreversible processes such as transport phenomena and rates of chemical reactions [²⁵]. It is built upon the first and second laws of thermodynamics, i.e., the conservation of energy and the entropy laws, respectively [²⁸, ²⁹]. These fundamental laws are formulated for irreversible processes based on the Local Equilibrium Hypothesis [²⁵, ²⁸, ²⁹]. According to such a hypothesis, “the local and instantaneous relations between thermodynamic quantities in a system out of equilibrium are the same as for a uniform system in equilibrium” [²⁵]. In other words, it assumes that a system can be viewed as composed of subsystems, which are large enough for microscopic fluctuations to be negligible, but small enough for equilibrium to be a good approximation [²⁵, ³⁰]. The local equilibrium assumption results in two immediate consequences: (1) All physical-chemical variables defined in CT are valid in NT, but they are replaced by their respective densities, which are continuous functions of space coordinates and time; (2) The equations that relate state variables in equilibrium are locally valid for non-equilibrium processes, e.g., Gibbs’ relation [²⁵, ²⁸, ²⁹]:

(1)$$TdS=dU+pdV-{\displaystyle \sum _{k=1}^n{\mu }_{k}d{N}_k.}$$

It is worth mentioning that the Local Equilibrium Hypothesis is valid for systems where the Maxwellian distribution is maintained, otherwise, generalizations should be made [²⁵].

Considering the information exposed previously, we can provide a local formulation for the laws of thermodynamics [²⁸]. Let us initially address the conservation law. To describe irreversible processes, such as mass diffusion, thermal conduction, and chemical reactions, it is convenient to define macroscopic conservation laws of mass, energy, and momentum in a differential form [²⁵, ²⁶, ²⁸, ²⁹]. Following that, assuming a multi-component system with the presence of conserved external forces and chemical reactions, the conservation of these physical quantities is locally rewritten as a balance equation (also known as continuity equation), given in the following general form [²⁵, ²⁶, ²⁸],

(2)$$\frac{dz(r,t)}{\partial t}=-\nabla \cdot J_z+{\sigma }_z.$$

A deduction of the equation above can be found in Appendix A.

In Eq. (2), z is the density of a generic physical quantity (z per unit volume), J_z is the flux term expressing the exchange of z between the system and surroundings, and σ_z is the internal production of z. Note that, to obtain an explicit expression for z we shall identify J_z and σ_z, which are particular to each problem. This is done with the aid of the second law.

In the NT formalism, the changes in entropy are expressed as the sum of two parts [²⁵, ²⁶, ²⁸, ²⁹]:

(3)$$dS=d_{e}S+d_{i}S.$$

Where d_eS accounts for the entropy change due to reversible processes and d_iS is the entropy change of irreversible processes. To be a valid local formulation Eq. (3) must fulfill the requirements of the second law. Following Clausius’s statement, such a law asserts that an isolated system can only evolve from an initial state to a final state if the variation of entropy is greater or equal to zero [²⁷, ³¹, ³²]. Therefore, as long as d_es = 0 for an isolated system, the variation of entropy of irreversible processes must be [²⁵, ²⁸, ²⁹]

(4)$$d_{i}S\ge 0.$$

Notice that, such inequality holds for any process, whereas d_eS may be positive, zero, or negative [²⁵, ²⁶, ²⁸]. Upon the assumption that in non-equilibrium conditions the densities of extensive variables are continuous functions of space and time, entropy shall be given as S = ∫ sdV, with s the entropy density [²⁸]. Following that statement, we may rewrite Eq. (3) as the rate of variation of entropy [²⁵, ²⁸, ²⁹]

(5)$$\frac{dS}{dt}=\frac{d_{e}S}{dt}+\frac{d_{i}S}{dt}.$$

Basically, the equation above tells that the entropy of a material body with fixed mass, and volume and bounded by a surface is given as the sum of the exchange with the surroundings, i.e., $\frac{d_{e}S}{dt}$, and the internal production, i.e., $\frac{d_{i}S}{dt}$ [²⁵]. Owing to the physical interpretation of such quantities, it is advantageous to write them in terms of the entropy flux J_s, i.e., the entropy crossing the boundary surface per unit of area and unit of time, and the rate of entropy production σ_s, i.e., the entropy produced per unit of volume and unit of time inside the system [²⁵, ²⁸, ²⁹]

(6)$$\frac{d_{e}S}{dt}=-{\displaystyle {\int }_A(J_s\cdot n)dAand}\frac{d_{i}S}{dt}={\displaystyle {\int }_V{\sigma }_{s}dV.}$$

One can realize that Eqs. (6) give a precise general form to compute the entropy of reversible and irreversible processes. Moreover, using Eqs. (6) in Eq. (5) and applying Gauss’ theorem, we also obtain a balance equation, as Eq. (2), for entropy density [²⁸].

To complete our description of the second law and to be able to obtain explicit governing equations of irreversible processes, we have to relate d_iS with the various irreversible phenomena that may take place inside the system [²⁸]. To this end, further consideration of the rate of entropy production, σ_s, must be done. Such a term is expressed as a sum of products of thermodynamic fluxes, J_i, and thermodynamic forces, X_i[²⁵, ²⁸, ²⁹, ³³, ³⁴]:

(7)$${\sigma }_s={\displaystyle \sum _i{J}_i{X}_i.}$$

We can determine Eq. (7) by comparing the entropy balance equation with the resulting entropy evolution equation obtained from substituting Eq. (2) into Gibbs’s relation, Eq. (1), see Appendix B [²⁵, ²⁸]. We highlight that X_i are not forces in the mechanical sense, they are associated to the gradients of intensive variables [²⁵].

Moreover, it is also known, from empirical observations of a large class of irreversible processes, that fluxes are linear functions of forces. Such a relationship is given by Eq. (8) and is named phenomenological relations[²⁵, ²⁸, ²⁹]:

(8)$$J_i={\displaystyle \sum _k{L}_{ik}{X}_k.}$$

Where L_ik are the phenomenological coefficients and they were introduced by Lars Onsager [³³, ³⁴]. Introducing Eq. (8) into Eq. (7), one gets a quadratic expression in the thermodynamic forces, i.e., σ_s = ∑_ik L_ikX_iX_k. Therefore, considering σ_s > 0, as a consequence of the second law, the entropy production must be positive definite. To hold such a property, the phenomenological coefficients shall satisfy the following conditions: L_ii ≥ 0 and $L_{ii}{L}_{kk}\ge \frac{1}{4}{(L_{ik}+L_{ki})}^2$[²⁵, ²⁶, ²⁸, ²⁹, ³³, ³⁴].

It is important to call the reader’s attention to a particularity of Eq. (7). In principle, one is allowed to couple any flux to any force from Eq. (7). However, material symmetry, also known as the Curie symmetry principle, restricts some of them. Such a principle asserts that for an isotropic system, thermodynamic fluxes and forces of different tensorial characters do not couple. This prevents macroscopic causes from having more elements of symmetry than the effects they generate [²⁵, ²⁸, ³⁵, ³⁶]. The consequences of material symmetry are very important to the understanding of Soret, Dufour, Seebeck, Peltier, and Thomson effects, for example [²⁵, ²⁸].

Now we are capable of obtaining explicit evolution equations from the combination of Gibbs’ relation and the non-equilibrium equations. This is carefully done in the next section.

3. Transport of Matter

From the NT theory, we begin with the physical derivation of the CH equation. The first step is to determine a general dynamical equation for describing the transport of matter in a binary, incompressible, and non-reactive mixture. To achieve such a task, we consider a closed and isothermal system with a total volume V. From these considerations, let c_α(r, t) and c_β(r, t) denote the local chemical concentrations, i.e., c_k = N_k/V with N_k the number of particles of the k-component, with k = α, β. Following the assumptions imposed on the system, the chemical concentrations obey c_α + c_β = C, with C ∈ ℝ⁺, thus $\frac{\partial c_{\alpha }}{\partial t}=-\frac{\partial c_{\beta }}{\partial t}$. With that, we can write from Eq. (2) the continuity equations of α and β:

(9)$$\frac{\partial c_k}{\partial t}=-\nabla \cdot J_k,k=\alpha ,\beta .$$

Note that, ${\sigma }_{s_k}=0$ because α and β do not chemically react and the flux J_k can be split in advected and diffusive parts [²],

(10)$$J_k=c_{k}v+j_k.$$

Where v is the volume flow velocity, which obeys the incompressibility condition, i.e., ∇ · v = 0. To ensure the matter is conserved we have that ∑_k m_kj_k = 0, with m_k the molecular mass of k-component [², ³]. This implies that j_α + j_β = 0, see the Appendix C. Therefore, using the information presented previously and restricting ourselves to a reference frame where v = 0 [²⁵], we can write the time evolution of the chemical concentration of component α as

(11)$$\frac{\partial c}{\partial t}=-\nabla \cdot j,$$

where we have abbreviated c = c_α for simplicity.

Now, our goal is to derive an expression for j. From Gibbs’ relation for the present problem, we can write the entropy as [²⁵]

(12)$$T\frac{ds}{dt}=-({\mu }_{\alpha }-{\mu }_{\beta })\frac{dc}{dt}.$$

However, we also know from Eq. (2) that

(13)$${\sigma }_s=\frac{ds}{dt}+\nabla \cdot J_s\ge 0.$$

So, using Eqs. (11) and (12) into Eq. (13) we have

(14)$$T{\sigma }_s=\overline{\mu }\nabla \cdot j+\nabla \cdot J_s\ge 0.$$

In Eq. (14) $\overline{\mu }={\mu }_{\alpha }-{\mu }_{\beta }$ is the exchange chemical potential [³]. Continuing, rearranging the last equation we obtain,

(15)$${\sigma }_s=-j\cdot \nabla \left(\frac{\overline{\mu }}{T}\right)+\nabla \cdot \left[J_s+j\left(\frac{\overline{\mu }}{T}\right)\right].$$

Since σ_s represents the entropy production inside the system, such an expression cannot have a flux term as $\nabla \cdot [{\mathbf{J}}_s+\mathbf{j}(\overline{\mu }/T)]$, which describes the rate of exchange with the neighborhood. Therefore, this term must be zero [²⁵], requiring

(16)$$J_s=-j\left(\frac{\overline{\mu }}{T}\right).$$

As a consequence of that, we have the rate of entropy production given as

(17)$${\sigma }_s=-j\cdot \nabla \left(\frac{\overline{\mu }}{T}\right).$$

From Eq. (7), we know the σ_s is given as a product of a thermodynamic flux j and a thermodynamic force $X=-\nabla \overline{\mu }/T$. Moreover, it follows from Eq. (8) that the flux in the linear response regime is proportional to the thermodynamic force, implying

(18)$$j=-\text{λ}(c)\nabla \left(\frac{\overline{\mu }}{T}\right).$$

This is an important result. We conclude from the equation above that the flow of matter in an isothermal system with the absence of advection is driven by the gradient of the chemical potential [²⁶]. Here, λ(c) denotes a mobility coefficient and it is positive to ensure σ_s ≥ 0 [³, ²⁵, ²⁶, ²⁸]. Now, using Eq. (18) into Eq. (11), we get the general dynamic equation for the transport of matter

(19)$$\frac{\partial c}{\partial t}=\nabla \cdot (\mathrm{\Lambda }(c)\nabla \overline{\mu }),$$

where Λ(c) = λ(c)/T [³, ²⁵]. As we can realize from Eq. (19), to determine a dynamical equation for the transport of matter to a specific problem, we just have to find out $\overline{\mu }$. For example, let us assume the mixture is ideal, then the chemical potential can be written as μ(p, T, c) = μ(p, T) + k_BT ln(c), where k_B is the Boltzmann constant [³¹, ³², ³⁷]. Substituting the chemical potential in Eq. (19) and taking Λ(c) = Λ₀c we obtain the famous diffusion equation [²⁵]:

(20)$$\frac{\partial c}{\partial t}=D{\nabla }^{2}c.$$

With D = Λ₀k_BT. The solution of such an equation describes the mixing process of two fluids; however, it is incapable of representing the process of phase separation. The reason is that the model is based on ideal fluids, which by definition lack internal structure. In other words, the particle-particle interaction is neglected. In this scenario, the chemical potential is only composed of entropic contributions, which push the system to the mixed state. In conclusion, ideal fluids do not separate. Therefore, to observe phase separation or, equivalently, to obtain the CH equation, we have to determine an exchange chemical potential for a mixture of real fluids.

4. Free Energy of Homogeneous Systems

As we argued in the last section, Eq. (19) dictates the spatiotemporal dynamics of a binary mixture of incompressible and non-reactive fluids in a closed and isothermal system. We also showed that, if the fluids are ideal, i.e., the interparticle interactions are equal to zero, we obtain the diffusion equation, Eq. (20). However, here we are interested in the dynamics of phase separation, which can only happen to real fluids. Therefore, the next step in the derivation of the CH equation is to determine the free energy of a binary homogeneous mixture of real fluids in a system with N, V, and T constants. We proceed with that in this section.

We begin by assuming the system can be discretized in a lattice of M sites, in a way that N_α + N_β = M. From this assumption, we may write the total volume as V = νM, where ν is the molecular volume, and the chemical concentrations as c_k = N_k/V = N_k/νM = ϕ_k/ν, where ϕ_k is the volume fraction of the k-component. Moreover, for simplicity, we take ν_α = ν_β = ν, in other words, the components α and β occupy the same volume in the lattice. Here, we use the binary variable σ_i to represent the occupancy of a particular site [³]. (Please, note that σ_i is different from the entropy production, σ_s, mentioned in the last sections. We keep this notation because it is usual in the treatment of Ising-like Hamiltonians.) More precisely, if σ_i = 1 then the i-site is occupied by α, and if σ_i = 0 it is occupied by β. From that construction, the energy of a specific configuration is given by the following Hamiltonian [², ³]

(21)$$H(\{\sigma \})=\frac{1}{2}{\displaystyle \sum _{〈i,j〉}({\in }_{\alpha \alpha }{\sigma }_i{\sigma }_j+{\in }_{\beta \beta }(1-{\sigma }_i)(1-{\sigma }_j)+{\in }_{\alpha \beta }[{\sigma }_i(1-{\sigma }_j)+{\sigma }_j(1-{\sigma }_i)]).}$$

Notice that, the Hamiltonian only accounts for short-range interactions such that the summation in Eq. (21) is over the nearest pairs ⟨i, j⟩ on the lattice and the factor $\frac{1}{2}$ avoids double counting [³]. The terms ϵ_αα, ϵ_ββ, and ϵ_αβ are the interaction parameters of α − α, β − β, and α − β bonds, respectively; such that, if ϵ_αα < 0, for example, the total energy of the system is reduced when molecules of α are neighbors, making these arrangements more probable [³, ³⁸]. We highlight that these parameters represent physical interactions of different nature, e.g., dipolar and van der Waals interactions, screened electrostatic interactions between charged molecular groups, or entropy-driven hydrophobic interactions [³].

Continuing, we want to obtain the exchange chemical potential ($\overline{\mu }$), which for a (N, V, T) system is given by

(22)$$\overline{\mu }={\left.\frac{\partial F}{\partial N_{\alpha }}\right|}_{T,V}={\left.-\frac{\partial F}{\partial N_{\beta }}\right|}_{T,V}={\left.\frac{\partial f}{\partial c}\right|}_T=v{\left.\frac{\partial f}{\partial \varphi }\right|}_T,$$

where F is the free energy of Helmholtz [³¹, ³², ³⁷]. It is known from the CT that in the conditions established here, this free energy is expressed as

(23)$$F(N_{\alpha {,}_{\beta }},V,T)=-k_{B}T\ln Z(N_{\alpha {,}_{\beta }},V,T).$$

In the equation above, Z is the canonical partition function and it can be written as follows [³, ³¹, ³⁷]

(24)$$Z={\displaystyle \sum _{\Omega }exp\left[-\frac{H({\sigma }_1,\dots ,{\sigma }_M)}{k_{B}T}\right],}$$

where Ω is the set of all possible configurations, assuming the molecules of the same type are indistinguishable from each other [³]. Therefore, to determine $\overline{\mu }$ we shall evaluate the free energy F and the partition function Z. However, defining the exact solution for Z and, consequently, F is a difficult task in most cases. Thus, to overcome this problem, we use a mean-field approximation. Essentially, such an approach is based on the reduction of an N-body interaction system to a 1-body system with an average or effective interaction field [³⁹]. It results from that the probability of any site being occupied by a molecule α (⟨σ_i⟩) is equal to the average of molecules of α in the entire system, i.e., ⟨σ_i⟩ = N_α/M = ϕ. Moreover, following such approximation the spatial correlations are neglected, in a way that, ⟨σ_iσ_j⟩ ≈ ⟨σ_i⟩⟨σ_j⟩ = ϕ² [³, ³⁷]. Making use of the mean-field theory, we can conclude that the internal energy can be written as

(25)$$E(\varphi )\approx \mathrm{\langle }H\mathrm{\rangle }=\frac{zM}{2}[{\mathrm{\in }}_{\alpha \alpha }{\varphi }^2+{\mathrm{\in }}_{\beta \beta }{(1-\varphi )}^2+2{\mathrm{\in }}_{\alpha \beta }\varphi (1-\varphi )],$$

where z is the coordination number per lattice site, and zM/2 is the number of nearest neighbors [³]. From the same perspective, the partition function is given by

(26)$$Z\approx \left|\text{Ω}\right|exp\left[-\frac{E(\varphi )}{k_{B}T}\right].$$

In the equation above, $\left|\text{Ω}\right|=\frac{M!}{N_{\alpha }!N_{\beta }!}$ is the number of all possible arrangements on the lattice. From Eqs. (25) and (26), and using Stirling’s approximation, i.e., ln N! ≈ N ln N − N, we can calculate the following total free energy density:

(27)$$\mathfrak{F}(\varphi )\approx \frac{z}{2v}[{\mathrm{\in }}_{\alpha \alpha }{\varphi }^2+{\mathrm{\in }}_{\beta \beta }{(1-\varphi )}^2+2{\mathrm{\in }}_{\alpha \beta }\varphi (1-\varphi )]+\frac{k_{B}T}{v}[\varphi \ln \varphi +(1-\varphi )\ln (1-\varphi )].$$

Bear in mind that, $\mathfrak{F}(\varphi )=F/Mv$, once V = Mν. Moreover, we may realize that the total free energy density is expressed as $\mathfrak{F}(\varphi )=f_0(\varphi )+f(\varphi )$, where f₀(ϕ) is the free energy density of the non-interactive system, i.e., f₀(ϕ) = ϕf(1)+(1 − ϕ)f(0), and f(ϕ) is the free energy density of mixing, which we can explicitly write as

(28)$$f(\varphi )=\frac{k_{B}T}{v}[\varphi \ln \varphi +(1-\varphi )\ln (1-\varphi )+\chi \varphi (1-\varphi )].$$

In the Eq. (28), χ = (z/2k_BT)(2ϵ_αβ − ϵ_αα − ϵ_ββ) is the interaction parameter [³, ⁴⁰, ⁴¹]. We can go a step further and realize that it is possible to determine the entropy of mixing using S = k_B ln |Ω| = −k_BM [ϕ ln ϕ + (1 − ϕ) ln(1 − ϕ)]. Therefore, from this last result, we conclude that f(ϕ) is composed of an entropic term related to the demixing process, and an energetic term, i.e., χϕ(1 − ϕ), representing the interactions between the molecules [³].

Now, before we proceed to the study of the dynamics driven by f(ϕ), let us analyze the physical characteristics of such a thermodynamic potential in terms of the stability of the equilibrium state. By definition, an equilibrium state is locally stable if it corresponds to a minimum of the free energy. In other words, the stability condition requires the free energy to be a locally convex function [³¹, ⁴²]. Such a criterion can be written in terms of the second derivative of the f(ϕ) (assuming the free energy is a continuous and smooth function), i.e., d²f(ϕ)/dϕ² ＞ 0.

From the definition of stability of the equilibrium state, let us analyze f(ϕ) as a function of ϕ for different χ presented in Figure 1. It is evident from this Figure that when χ is lower than a critical value χ_c, which indicates that the energetic interaction between particles of the same composition is weak, the free energy is a convex function of ϕ and, therefore, the mixed equilibrium state is favored. However, as χ is increased, meaning stronger energetic interaction between particles of the same composition, a local maximum appears in the free energy from χ ＞ χ_c, forming a double-well potential. For the range of ϕ that the local maximum lies on, the mixed state is unstable. Therefore, to minimize the free energy of the equilibrium state, the system spontaneously separates into two phases with compositions ϕ₁ and ϕ₂, which are the two minima represented in Figure 2A. We have to point out that the coexistence of both phases in the equilibrium requires that the exchange chemical potential and the osmotic pressure between the phases are equilibrated, i.e., ${\overline{\mu }}_1={\overline{\mu }}_2$ and Π₁ = Π₂ [², ³, ⁴²]. These requirements are satisfied by a common tangent to the points (ϕ₁, f(ϕ₁)) and (ϕ₂, f(ϕ₂)). This last observation means we can determine the equilibrium compositions of each phase geometrically, through the “common tangent rule” [³, ⁴²]. Such a procedure is also known as Maxwell’s tangent construction [³].

Figura 1
The free energy density of mixing for an incompressible binary fluid as a function of ϕ for different χ. As one might realize, for χ ≥ 2, two local minima arise, indicating the mixed state is unstable.

Figura 2
A)A representation of the free energy density of mixing for χ ≥ 2, indicating the stable, metastable, and unstable regions. B) Phase diagram of ϕ × χ. The spinodal and binodal curves are shown in a red dashed line and a black thick line, respectively.

We can proceed and look in detail f(ϕ) with χ ＞ χ_c. Under such conditions, we may divide the free energy into two regions: (1) the miscibility gap, i.e., ϕ₁ ＜ ϕ ＜ ϕ₂, where the system minimizes the free energy by separating into two phases; (2) the stable region, i.e., ϕ ＜ ϕ₁ and ϕ ＞ ϕ₂, where the homogeneous state is maintained as presented in Figure 2A [⁴²]. We may also categorize two different domains inside the miscibility gap following the stability characteristics of the equilibrium state, i.e., the unstable and the metastable domains. These domains are separated by inflection points. For compositions in the unstable region, the mixed state spontaneously separates in the presence of small fluctuations, through a mechanism known as spinodal decomposition [³, ⁴²]. On the other hand, the metastable region is locally stable for small composition fluctuations, i.e., the system does not separate spontaneously. However, the presence of large fluctuations can induce demixing, which in this regime occurs by a mechanism called nucleation and growth [³, ⁴²].

We can precisely define the unstable and metastable regions using the criteria for stability. Initially, we can determine the necessary conditions for two distinct phases to coexist, i.e., df(ϕ)/dϕ = 0. Such an equation can be solved for the interaction parameter χ, see Figure 2B. These conditions define the gray region in this Figure, which is delimited by the coexistence curve or binodal curve. The reader must realize this curve delimits necessary, but not sufficient, conditions for phase separation. To be sufficient, the state must be unstable, which requires d²f(ϕ)/dϕ² ＜ 0. Solving this last equation for χ we obtain the spinodal curve, represented in red dashed lines in Figure 2B. Therefore, the system will separate into two phases for any pair of (ϕ, χ) chosen inside the spinodal curve. The region between the binodal and spinodal curves is metastable. Out of these regions, the system keeps mixed [⁴²].

Considering what we have exposed until now, we could proceed with the study of the phase dynamics associated with the free energy of mixing, given by Eq.(28). To do that, we simply have to calculate the exchange chemical potential and use it in Eq. (19), which would result in the CH equation. However, the free energy density determined above neglects the energy associated with the formation of interfaces between two regions with different compositions, in other words, the resulting CH equation is incomplete. As a consequence, the simulation of this version of the model will only result in short-range structures. Therefore, to accurately describe the process of phase separation and observe the emergence of long-range structures, e.g., droplets, we have to include the effects of the interface in the physical model [⁴³].

5. Free Energy of a Non-Homogeneous System

As we pointed out at the end of the last section, to effectively describe the process of phase separation the interfacial energy must be accounted for. Therefore, this section aims to include such energetic contributions to the physical model. This is the final step to obtain the CH equation.

In a system formed by two phases, the composition is spatially non-uniform. Following this observation, it is reasonable to suppose that the free energy density, f, should depend both on the local composition and the composition of the immediate environment [⁴³]. Therefore, we will express f as a function of the local composition and the local derivatives, i.e., f = f(ϕ,∇ϕ,∇²ϕ, . . .). We can expand this last expression in a Taylor series about $\overrightarrow{\varphi }=({\varphi }_0,0,0,\dots )$, resulting in [³⁸, ⁴³]

(29)$$f(\varphi ,{\varphi }_1,{\varphi }_2,{\varphi }_3,{\varphi }_{11},{\varphi }_{22},{\varphi }_{33},{\varphi }_{12},{\varphi }_{23},{\varphi }_{31},\dots )$$$$=f{(\varphi )}_0+{\displaystyle \sum _i^3\frac{\partial f({\varphi }_0)}{\partial {\varphi }_i}{\varphi }_i+}{\displaystyle \sum _{i,j=1}^3\frac{\partial f({\varphi }_0)}{\partial {\varphi }_{ij}}{\varphi }_{ij}}$$$$\mathrm{}+\frac{1}{2}{\displaystyle \sum _{i,j=1}^3\frac{{\partial }^{2}f({\varphi }_0)}{\partial {\varphi }_i\partial {\varphi }_j}{\varphi }_i{\varphi }_j+}\frac{1}{2}{\displaystyle \sum _{i,j,k=1}^3\frac{{\partial }^{2}f({\varphi }_0)}{\partial {\varphi }_k\partial {\varphi }_{ij}}{\varphi }_k{\varphi }_{ij}}$$$$+\frac{1}{2}{\displaystyle \sum _{i,j,k,l=1}^3\frac{{\partial }^{2}f({\varphi }_0)}{\partial {\varphi }_{ij}\partial {\varphi }_{kl}}{\varphi }_{ij}{\varphi }_{kl}+\mathcal{O}{(\varphi )}^3.}$$

Here, we ignore terms higher than second-order derivatives. Note that, in Eq. (29) ${\varphi }_i=\frac{\partial \varphi }{\partial x_i}$ and ${\varphi }_{ij}=\frac{{\partial }^2\varphi }{\partial x_i\partial x_j}$. Considering the system is isotropic, then the free energy must be invariant to the symmetry operations of reflection (x_i → −x_i) and rotation about a fourfold axis (x_i → x_i) [³⁸, ⁴³]. Hence,

(30)$$\frac{\partial f({\varphi }_0)}{\partial {\varphi }_i}=0,\frac{\partial f({\varphi }_0)}{\partial {\varphi }_{ii}}={\kappa }_1,\frac{{\partial }^{2}f({\varphi }_0)}{\partial {\varphi }_i^2}={\kappa }_2,i=1,2,3$$

(31)$$\frac{\partial f({\varphi }_0)}{\partial {\varphi }_{ij}}=\frac{{\partial }^{2}f({\varphi }_0)}{\partial {\varphi }_i\partial {\varphi }_j}=0,i\ne j.$$

Using the assumptions made above, the Eq. (29) can be rewritten as

(32)$$f(\varphi ,{\varphi }_1,{\varphi }_2,{\varphi }_3,{\varphi }_{11},{\varphi }_{22},{\varphi }_{33},{\varphi }_{12},{\varphi }_{23},{\varphi }_{31},\dots )=f{(\varphi )}_0+{\kappa }_1{\nabla }^2\varphi +\frac{{\kappa }_2}{2}{\left|\nabla \varphi \right|}^2.$$

We can determine the total free energy integrating Eq. (32) over the system’s volume,

(33)$$F[\varphi ]={\displaystyle \int dV\mathrm{}}\left[f{(\varphi )}_0+{\kappa }_1{\nabla }^2\varphi +\frac{{\kappa }_2}{2}{\left|\nabla \varphi \right|}^2\right].$$

Applying the divergence theorem to the second term of the right-hand side of the last equation and assuming the term ∇ϕ . n is zero at the boundary, we obtain the following expression [³⁸, ⁴³],

(34)$$F[\varphi ]={\displaystyle \int dV\mathrm{}}\left[f{(\varphi )}_0+\kappa {\left|\nabla \varphi \right|}^2\right].$$

Where, $\kappa =\frac{{\kappa }_2}{2}-\frac{\partial {\kappa }_1}{\partial \varphi }$, as we show in Appendix D. The functional given in Eq. (34) represents the total free energy of the system and it is the central one in our treatment. Observe that such an equation accounts for the free energy of mixing of the homogeneous system, i.e., f(ϕ)₀ = f(ϕ) given by Eq. (28), and the gradient energy related to the formation of interfaces that is a function of the local composition, i.e., |∇ϕ|².

Now we are able to calculate the exchange chemical potential and obtain the final equation to describe the dynamics of phase separation of a binary mixture, i.e., the CH equation.

6. The Cahn-Hilliard Equation

In this section, we use the information obtained previously to determine the CH equation. After that, we perform numerical simulations and analyze the outcomes.

The exchange chemical potential is given by:

(35)$$\overline{\mu }=\frac{\delta F[\varphi ]}{\delta \varphi }=\frac{\partial g(\varphi ,\nabla \varphi )}{\partial \varphi }-\nabla \cdot \frac{\partial g(\varphi ,\nabla \varphi )}{\partial \nabla \varphi }.$$

Where, g(ϕ,∇ϕ) = f(ϕ) + κ|∇ϕ|² and $\frac{\delta F[\varphi ]}{\delta \varphi }$ is a functional derivative (more information about it can be found in the references) [³, ⁴³,⁴⁴,⁴⁵]. Before calculating the chemical potential, it is convenient and a common practice to replace f(ϕ) with a fourth-order polynomial that keeps the form of a double-well potential, see Figure 3, such as [³, ³⁸]

(36)$$f(\varphi )=-\frac{b}{2}{(\varphi -{\varphi }_c)}^2+\frac{a}{4}{(\varphi -{\varphi }_c)}^4.$$

Figura 3
Fourth-order polynomial free energy, Eq. (36), with a = b = 1.0 and ϕ_c = 1.0, which results in ϕ ∈ (0, Φ) with Φ = 2.0.

The replacement of f(ϕ) is just a mathematical procedure to avoid numerical divergence. One can notice that the Laplacian of the exchange chemical potential calculated from Eq. (28) contains terms that can lead to singularity problems.

The last equation is obtained by expanding Eq. (28) around the critical volume fraction ϕ_c and assuming, for simplicity, the special case that the first and third-order components of the expansion are equal to zero [³]. Furthermore, the parameters a and b are related to entropic and energy terms [³]. Using Eq. (36) in Eq. (34) we have the Landau-Ginzburg free energy [⁴⁶]:

(37)$$F\left[\varphi \right]={\displaystyle \int dV\mathrm{}}\left[-\frac{b}{2}{(\varphi -{\varphi }_c)}^2+\frac{a}{4}{(\varphi -{\varphi }_c)}^4+\kappa {\left|\nabla \varphi \right|}^2\right].$$

Finally, calculating $\overline{\mu }$ from Eq. (37) and substituting it into the Eq. (19), we obtain the Cahn-Hilliard equation,

(38)$$\frac{\partial \varphi }{\partial t}=M{\nabla }^2(-b(\varphi -{\varphi }_c)+a{(\varphi -{\varphi }_c)}^3+\tilde{\kappa }{\nabla }^2\varphi ).$$

In which we have assumed the mobility (Λ(ϕ)ν) is a constant and equal to M, and $\tilde{\kappa }=2\kappa$.

The reader can immediately recognize that Eq. (38) is a nonlinear partial differential equation with second and fourth-order derivatives in the space coordinates. In general, such a class of equations does not have analytical solutions, requiring the usage of numerical methods. There is a great number of works in the literature presenting efficient numerical approaches based on finite difference [⁴⁷,⁴⁸,⁴⁹], finite element [⁵⁰,⁵¹,⁵²], Fourier-spectral [⁵³,⁵⁴,⁵⁵], cell dynamics system (CDS) [⁵⁶, ⁵⁷] methods to deal with the CH equation.

Here, we have employed a simple scheme of explicit finite differences (EFD), which can be easily implemented in Fortran and the simulations can be performed on a laptop. Summarily, the EFD method discretizes both spatial domain and time interval into a grid of integer numbers. In this construction, the derivatives are approximated by finite difference formulas, and the differential equation is transformed into a system of algebraic equations, i.e., one equation for each point of the grid [⁵⁸].

This procedure can be applied to Eq. (38) to obtain its algebraic version. To do that, we have considered a system with size L_x × L_y, where L_x and L_y are the system’s length in the x and y-direction. We can discretize it into a grid of N_x × N_y points, where N_x = L_x/Δx, N_y = L_y/Δy, and Δ_x,y are the distances between two consecutive points in both directions, i.e., the spatial steps. The same can be done with the time interval T, i.e., we can discretize it into S points, where S = T/Δt and Δt is the time step. In this discrete space $\varphi (x,y,t)\approx {\varphi }_{i,j}^n$, where n ∈ S, i ∈ N_x, and j ∈ N_y [⁵⁸]. Therefore, the algebraic version of the CH equation (38) is given as,

(39)$$\frac{\partial \varphi }{\partial t}=M\left[-b\left(\frac{{\partial }^2\varphi }{\partial x^2}+\frac{{\partial }^2\varphi }{\partial y^2}\right)+3a({\varphi }_c-\varphi )\left(-2\left(\frac{\partial {\varphi }^2}{\partial x}+\frac{\partial {\varphi }^2}{\partial y}\right)+({\varphi }_c-\varphi )\left(\frac{{\partial }^2\varphi }{\partial x^2}+\frac{{\partial }^2\varphi }{\partial y^2}\right)\right)+\tilde{\kappa }\left(\frac{{\partial }^4\varphi }{\partial x^4}+2\frac{{\partial }^4\varphi }{\partial x^2\partial y^2}+\frac{{\partial }^4\varphi }{\partial y^4}\right)\right].$$

Where the partial derivatives are approximated by finite differences:

(40)$$\frac{\partial \varphi }{\partial t}\approx \frac{{\varphi }_{i,j}^{n+1}-{\varphi }_{i,j}^n}{\Delta t}$$

(41)$$\frac{\partial \varphi }{\partial x}\approx \frac{{\varphi }_{i+1,j}^n-{\varphi }_{i-1,j}^n}{2\Delta x}$$

(42)$$\frac{\partial \varphi }{\partial y}\approx \frac{{\varphi }_{i,j+1}^n-{\varphi }_{i,j-1}^n}{2\Delta y}$$

(43)$$\frac{{\partial }^2\varphi }{\partial x^2}\approx \frac{{\varphi }_{i+1,j}^n-2{\varphi }_{i,j}^n+{\varphi }_{i-1,j}^n}{\Delta x^2}$$

(44)$$\frac{{\partial }^2\varphi }{\partial y^2}\approx \frac{{\varphi }_{i,j+1}^n-2{\varphi }_{i,j}^n+{\varphi }_{i,j-1}^n}{\Delta y^2}$$

(45)$$\frac{{\partial }^4\varphi }{\partial x^4}\approx \frac{{\varphi }_{i+2,j}^n-4{\varphi }_{i+1,j}^n+6{\varphi }_{i,j}^n-4{\varphi }_{i-1,j}^n+{\varphi }_{i-2,j}^n}{\Delta x^4}$$

(46)$$\frac{{\partial }^4\varphi }{\partial y^4}\approx \frac{{\varphi }_{i,j+2}^n-4{\varphi }_{i,j+1}^n+6{\varphi }_{i,j}^n-4{\varphi }_{i,j-1}^n+{\varphi }_{i,j-2}^n}{\Delta y^4}$$

(47)$$\frac{{\partial }^4\varphi }{\partial x^2\partial y^2}\approx \left(\frac{1}{16\Delta x^2\Delta y^2}\right)({\varphi }_{i+2,j+2}^n-2{\varphi }_{i,j+2}^n+{\varphi }_{i-2,j+2}^n-2{\varphi }_{i+2,j}^n+4{\varphi }_{i,j}^n-2{\varphi }_{i-2,j}^n+{\varphi }_{i+2,j-2}^n-2{\varphi }_{i,j-2}^n+{\varphi }_{i-2,j-2}^n).$$

One can realize that the approximated solution of the partial differential equation in a discrete point of the grid (i, j) for the time n + 1 is computed by solving an algebraic equation using information from previous time steps n. It is also important to emphasize that EFD has some limitations in terms of numerical stability. So, to avoid numerical divergence, the time and spatial steps must satisfy the Neumann criterion [⁵⁸].

The CH equation is solved by assuming L_x = L_y = 100, periodic boundary conditions, a = b = 1.0, $\tilde{\kappa }=0.5$, Δx = Δy = 0.5, Δt = 10⁻³, and different initial conditions ϕ₀ (which are displayed in the Figures’ captions). Moreover, all the simulations are carried out with ϕ_c = 1, which implies ϕ ∈ (0, Φ) with Φ = 2, as presented in Figure 3. Figure 4 shows snapshots of simulations performed with two different compositions.

Figura 4
A), B), C), and D) show the outcomes of the numerical simulation of the CH equation for ϕ₀ = 0.8 and M = 1 at t = 50, 250, 500, and 1000 time units, respectively. E), F), G), and H) show the outcomes of the numerical simulation of the CH equation for ϕ₀ = 1.35 and M = 1 at t = 50, 250, 500, and 1000 time units, respectively.

As we can note from the simulations, the topology of the fluid domains depends crucially on the volume fractions. We observe that after the spinodal decomposition, spherical droplets rapidly form and as the simulation evolves the mean droplet radius increases. Such a process is called droplet coarsening and it can proceed through two different mechanisms: Coalescence and Ostwald ripening [³, ⁵⁹]. The coalescence is driven by the presence of thermal noise. It induces Brownian motion, which increases the number of droplet-droplet collisions and promotes the formation of larger droplets. On the other hand, the driving force of Ostwald ripening is the Laplace pressure difference between droplets of different sizes. It creates a material transport flow from small to large droplets across the intervening continuous phase, maximizing the area of bigger droplets and minimizing, consequently, the energy [³, ⁵⁹]. It is worth realizing that both mechanisms occur through diffusion, but the nature of the diffusing species is quite different [⁵⁹]. Therefore, considering we neglected momentum transport due to thermal fluctuations, the coarsening mechanism observed in the simulations presented in Figure 4 is Ostwald ripening.

It is convenient in many situations to avoid or, at least to slow down, droplet coarsening. This can be achieved in passive and active systems.

6.1. Passive system

A passive system composed of a binary mixture of immiscible liquids in the presence of surfactants can slow down droplet coarsening.

Surfactants are amphiphilic organic compounds, i.e., molecules composed of both hydrophobic (tail) and hydrophilic (head) groups. These structural characteristics enable such molecules to interact with both polar and non-polar molecules simultaneously. Therefore, in a mixture of immiscible liquids, e.g., water and oil, surfactants are adsorbed at interfaces between coexisting phases, forming a “wall”. A direct consequence of the presence of surfactant molecules in such positions is the reduction of the interfacial energy, resulting in the stabilization of the droplets [⁵⁹].

Following that, the coalescence process may be reduced to negligible levels by the presence of surfactants [⁵⁹]. This is possible because the approximation of two droplets is prevented due to the repulsive force generated by the steric interactions between two surfactant tails. Such effects can be enhanced by considering charged surfactants, which create a Coulombic “barrier” moving droplets away from each other [⁵⁹]. Ostwald ripening can also be decreased in the presence of surfactants, once it is triggered by energetic differences between droplets; however, it can not be completely prevented [⁵⁹]. This happens because in nonideal situations there is a saturation point of adsorbed surfactant molecules at the interface of emulsions. At this point, there is the formation of structures called micelles and the interfacial energy cannot be further reduced [⁵⁹].

We may insert the effects of surfactants in the phase separation model by considering the resultant polarization vector field. The free energy accounting such effects is given below,

(48)$$F\left[\varphi ,\mathbf{P}\right]={\displaystyle \int dV\mathrm{}}\left[f{(\varphi )}_0+\kappa {\left|\nabla \varphi \right|}^2+\frac{{\left|\mathbf{p}\right|}^2}{2\chi }+\mathrm{\in }\mathbf{p}\cdot \nabla \varphi \right].$$

In the Eq. (48), p is the resultant polarization vector field, $\chi =〈{\left|\mathbf{p}\right|}^2〉\propto c_s$ is the osmotic compressibility, and c_s is the global concentration of surfactant molecules [⁵⁹]. The included terms represent the reduction in the free energy due to the presence of the surfactants in the solution, i.e., $\frac{{\left|\mathbf{p}\right|}^2}{2\chi }$ and their adsorption in the interface between the two phases, i.e., ϵp · ∇ϕ. Figure 5 shows the time evolution of a binary mixture of immiscible fluids in the presence of surfactants. It is clear from comparing Figures 4 and 5 that the presence of surfactants slows down Ostwald ripening; however, it does not suppress completely the coarsening, i.e., it quasi-suppresses Ostwald ripening.

Figura 5
Numerical simulation of the binary mixture of immiscible fluids and surfactant, described by Eq. (48), at t = A) 50, B) 250, C) 500, and D) 1000 time units. The conditions are ϕ₀ = 1.0, M = 1, p = (0, 0), χ = 11, and ϵ = 0.1. In this case, the simulations were performed with the xmds2 software package (version 2.2.2), which employs an adaptive fourth-fifth order Runge-Kutta and spectral methods for computing the spatial derivatives [⁶⁰].

6.2. Active system

An active system formed by the combination of chemical reactions and phase separation can suppress Ostwald ripening [²³].

Let us consider the scenario in which the components of the binary mixture, i.e., α and β, can convert into each other through a simple first-order reaction as represented below,

$$\alpha \underset{{\text{k}}_b}{\overset{{\text{k}}_f}{\rightleftharpoons }}\beta$$

Where k_f and k_b are the rate constants. The chemical reaction presented above indicates that the droplet material, α, chemically transforms into the background material, β, i.e., $\alpha \stackrel{{\text{k}}_f}{\to }\beta$. The same process may take place in the reverse direction, i.e., $\beta \stackrel{{\text{k}}_b}{\to }\alpha$. Notice that, if we assume the system can only react through these reversible steps, it would eventually reach the equilibrium state. Therefore, to keep the system out of equilibrium we add a source of energy, S, and a waste, W, in the following manner $\alpha +\text{S}\underset{{\text{k}}_b}{\overset{{\text{k}}_f}{\rightleftharpoons }}\beta +\text{W}$ [³]. This is motivated by biological systems, where S may represent the supply of chemical fuel such as adenosine triphosphate (ATP) [²³]. As we know, the chemical equilibrium is reached when Δμ = μ_α + μ_F − μ_β − μ_W = 0. So, as long as we keep μ_S − μ_W ≠ 0, the system will be away from thermodynamic equilibrium. This is possible because we would be able to control the chemical concentration of S and W externally.

We obtain the dynamic equation by applying the law of mass action and combining it with the CH equation, resulting in

(49)$$\frac{\partial {\varphi }_{\alpha }}{\partial t}=-k_f{\varphi }_{\alpha }{\varphi }_S+k_b{\varphi }_{\beta }{\varphi }_W+M{\nabla }^2\overline{\mu }.$$

As one can realize, the number of α and β molecules is constant, so that, ϕ_β = Φ − ϕ_α. And, for simplicity, we suppose that S and W are diluted and constant over time and space, so that, k_fϕ_αϕ_S = K_fϕ_α and K_bϕ_βϕ_W = K_bϕ_β. Using this considerations, and taking ϕ_α = ϕ, we have

(50)$$\frac{\partial \varphi }{\partial t}=K_f\varphi +K_b\left(\mathrm{\Phi }-\varphi \right)+M{\nabla }^2\overline{\mu }.$$

The outcomes of the numerical simulation of Eq. (50) are shown in Figure 6.

Figura 6
Numerical simulation of a binary mixture of immiscible fluids that can convert into each other through a first-order chemical reaction. The dynamics are described by Eq. (50) at t = A) 50, B) 250, C) 500, and D) 1000 time units. The conditions are ϕ₀ = 1.0, K_f = 9 × 10⁻³, K_b = 2 × 10⁻³ and M = 2.

It is evident from the results of the simulation that Ostwald ripening is suppressed in the presence of first-order chemical reactions. In fact, its direction is inverted, i.e., larger droplets shrink to form smaller droplets. Such a phenomenon is possible because when the system is driven out of the equilibrium by chemical reactions, it might evolve toward dynamically stable steady states. In this case, one such state is composed of stable monodisperse droplets. It is worth mentioning that these states can only be reached by the correct choice of physical parameters [²³].

7. Conclusion

In this paper, we have provided a detailed deduction of the Cahn-Hilliard equation, starting from the basic concepts of non-equilibrium thermodynamics to obtain a general equation for the transport of matter, Eq. (19), going through the general ideas of statistical mechanics to obtain the free energy density of mixing for a homogeneous system, Eq. (28), exploring the variational principle to extend such free energy and include the energetic contribution of interfaces, Eq. (34), obtaining an approximation for the free energy for a homogeneous system to avoid numerical problems, Eq. (36), calculating the exchange chemical potential from the Landau-Ginzburg free energy, Eq. (37), using the functional derivative, Eq. (35), and coming up by the end with the Cahn-Hilliard equation, by using Eq. (35) into Eq. (19). We then performed simple simulations of the Cahn-Hilliard equation, describing the coarsening processes and offering two alternatives to suppress, or at least decrease them: a passive system composed of surfactants, and an active system that combines chemical reactions with phase separation.

We expect that this material may offer the basic theoretical background on non-equilibrium thermodynamics, enabling the reader to deal with different irreversible phenomena close to the equilibrium in a continuous approximation. It also may serve as the starting point for studies on processes far from equilibrium, where non-linear chemical kinetics can induce the emergence of dissipative structures, such as pattern formation and chaos.

Acknowledgments

L. Silva-Dias thanks M.Sc. Barbara Bohlen, Dr. Asdrubal Lozada, Dra. Sthefany dos Santos Sena, and the referee for their helpful comments. L. Silva-Dias also thanks FAPESP, grant 2019/23205-3 and 2022/00257-0, and Coordenção de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), finance code 001.

Supplementary material

The following online material is available for this article:

A Fortran code with the explicit finite differences method to solve the Cahn-Hilliard equation:

Appendix A – General Balance Equation.

Appendix B – Entropy Production.

Appendix C – Conservation of Mass.

Appendix D – Boundary Terms.

References

• ^[1]
  J. Berry, C. Brangwynne and M. Haataja, Reports On Progress In Physics 81, 046601 (2018).
• ^[2]
  R. Seyboldt, The dynamics of chemically active droplets Masters Dissertation, Technische Universität Dresden, Dresden (2019).
• ^[3]
  C. Weber, D. Zwicker, F. Jülicher and C. Lee, Reports On Progress In Physics 82, 064601 (2019).
• ^[4]
  Y. Shin and C. Brangwynne, Science 357, eaaf4382 (2017).
• ^[5]
  D. Lohse and X. Zhang, Nature Reviews Physics 2, 426 (2020).
• ^[6]
  M. Cates and E. Tjhung, Journal Of Fluid Mechanics 836, P1 (2018).
• ^[7]
  E. Dufresne, H. Noh, V. Saranathan, S. Mochrie, H. Cao and R. Prum, Soft Matter 5, 1792 (2009).
• ^[8]
  A. Parnell, A. Washington, O. Mykhaylyk, C. Hill, A. Bianco, S. Burg, A. Dennison, M. Snape, A. Cadby and A. Smith, Scientific Reports 5, 18317 (2015).
• ^[9]
  S. Burg and A. Parnell, Journal Of Physics: Condensed Matter 30, 413001 (2018).
• ^[10]
  S. Mao, M. Chakraverti-Wuerthwein, H. Gaudio and A. Košmrlj, Physical Review Letters 125, 218003 (2020).
• ^[11]
  S. Banani, H. Lee, A. Hyman and M. Rosen, Nature Reviews Molecular Cell Biology 18, 285 (2017).
• ^[12]
  J. Kirschbaum and D. Zwicker, Journal Of The Royal Society Interface 18, 20210255 (2021).
• ^[13]
  A. Hyman, C. Weber and F. Jülicher, Annual Review Of Cell And Developmental Biology 30, 39 (2014).
• ^[14]
  M. Feric, N. Vaidya, T. Harmon, D. Mitrea, L. Zhu, T. Richardson, R. Kriwacki, R. Pappu and C. Brangwynne, Cell 165, 1686 (2016).
• ^[15]
  L. Demarchi, A. Goychuk, I. Maryshev and E. Frey, Physical Review Letters 130, 128401 (2023).
• ^[16]
  N. Wolf, J. Priess and D. Hirsh, Development 73, 297 (1983).
• ^[17]
  S. Strome and W. Wood, Cell 35, 15 (1983).
• ^[18]
  A. Doostmohammadi, J. Ignés-Mullol, J. Yeomans and F. Sagués, Nature Communications 9, 3246 (2018).
• ^[19]
  M. Marchetti, J. Joanny, S. Ramaswamy, T. Liverpool, J. Prost, M. Rao and R. Simha, Reviews Of Modern Physics 85, 1143 (2013).
• ^[20]
  A. Walther and A. Muller, Chemical Reviews 113, 5194 (2013).
• ^[21]
  R. Manna, A. Laskar, O. Shklyaev and A. Balazs, Nature Reviews Physics 4, 125 (2022).
• ^[22]
  L. Silva-Dias and A. Lopez-Castillo, The Journal Of Physical Chemistry Letters 13, 296 (2022).
• ^[23]
  D. Zwicker, A. Hyman and F. Jülicher, Physical Review E 92, 012317 (2015).
• ^[24]
  D. Zwicker, R. Seyboldt, C. Weber, A. Hyman and F. Jülicher, Nature Physics 13, 408 (2017).
• ^[25]
  G. Lebon, D. Jou and J. Casas-Vázquez, Understanding non-equilibrium thermodynamics (Springer, Berlin, 2008).
• ^[26]
  D. Kondepudi and I. Prigogine, Modern thermodynamics: from heat engines to dissipative structures (John Wiley & Sons, Chichester, 2014).
• ^[27]
  H. Nussenzveig, Curso de Física Básica: fluidos, oscilações e ondas, calor (Editora Blucher, São Paulo, 2018).
• ^[28]
  S. De Groot and P. Mazur, Non-equilibrium thermodynamics (Courier Corporation, Amsterdan, 2013).
• ^[29]
  I. Prigogine, Science 201, 777 (1978).
• ^[30]
  J. Vilar and J. Rubi, Proceedings Of The National Academy Of Sciences 98, 11081 (2001).
• ^[31]
  H. Callen, Thermodynamics and an Introduction to Thermostatistics (American Association of Physics Teachers, San Francisco, 1998).
• ^[32]
  D. McQuarrie and J. Simon, Physical chemistry: a molecular approach (University Science Books Sausalito, Sausalito, 1997).
• ^[33]
  L. Onsager, Physical Review 37, 405 (1931).
• ^[34]
  L. Onsager, Physical Review 38, 2265 (1931).
• ^[35]
  P. Curie, Journal De Physique Théorique Et Appliquée 3, 393 (1894).
• ^[36]
  A. Chalmers, The British Journal For The Philosophy Of Science 21, 133 (1970).
• ^[37]
  S. Salinas, Introdução a física estatística (Edusp, São Paulo, 1997), v. 9.
• ^[38]
  D. Lee, J. Huh, D. Jeong, J. Shin, A. Yun and J. Kim, Computational Materials Science 81, 216 (2014).
• ^[39]
  P. Chaikin, T. Lubensky and T. Witten, Principles of condensed matter physics (Cambridge University Press, Cambridge, 1995).
• ^[40]
  P. Flory, The Journal Of Chemical Physics 10, 51 (1942).
• ^[41]
  M. Huggins, The Journal Of Physical Chemistry 46, 151 (1942).
• ^[42]
  M. Rubinstein and R. Colby, Others Polymer physics (Oxford University Press, New York, 2003).
• ^[43]
  J. Cahn and J. Hilliard, The Journal Of Chemical Physics 28, 258 (1958).
• ^[44]
  N. Lemos, Mecânica analítica (Editora Livraria da Física, São Paulo, 2007).
• ^[45]
  H. Goldstein, C. Poole and J. Safko, Classical mechanics (American Association of Physics Teachers, San Francisco, 2002).
• ^[46]
  M. Kardar, Statistical physics of fields (Cambridge University Press, Cambridge, 2007).
• ^[47]
  J. Shin, D. Jeong and J. Kim, Journal Of Computational Physics 230, 7441 (2011).
• ^[48]
  Y. Li, D. Jeong, J. Shin and J. Kim, Computers & Mathematics With Applications 65, 102 (2013).
• ^[49]
  S. Aland, J. Lowengrub and A. Voigt, Computer Modeling In Engineering & Sciences: CMES 57, 77 (2010).
• ^[50]
  A. Diegel, C. Wang and S. Wise, IMA Journal Of Numerical Analysis. 36, 1867 (2016).
• ^[51]
  A. Diegel, X. Feng and S. Wise, SIAM Journal On Numerical Analysis 53, 127 (2015).
• ^[52]
  A. Diegel, C. Wang, X. Wang and S. Wise, Numerische Mathematik 137, 495 (2017).
• ^[53]
  D. Li and Z. Qiao, Journal Of Scientific Computing 70, 301 (2017).
• ^[54]
  D. Li, Z. Qiao and T. Tang, SIAM Journal On Numerical Analysis 54, 1653 (2016).
• ^[55]
  D. Li and Z. Qiao, Communications In Mathematical Sciences 15, 1489 (2017).
• ^[56]
  Y. Oono and S. Puri, Physical Review Letters 58, 836 (1987).
• ^[57]
  Y. Oono and S. Puri, Physical Review A 38, 434 (1988).
• ^[58]
  G. Smith, Numerical solution of partial differential equations: finite difference methods (Oxford University Press, Oxford, 1985).
• ^[59]
  M. Cates, Soft Interfaces: Lecture Notes Of The Les Houches Summer School 98, 317 (2017)
• ^[60]
  G. Dennis, J. Hope and M. Johnsson, Computer Physics Communications 184, 201 (2013).

Publication Dates

• Publication in this collection
  08 Jan 2024
• Date of issue
  2023

History

• Received
  22 Sept 2023
• Reviewed
  14 Nov 2023
• Accepted
  16 Nov 2023