A phase-field study on polymerization-induced phase separation occasioned by diffusion and capillary flow—a mechanism for the formation of porous microstructures in membranes

The performance and the application of membranes, which are usually produced from polymer solutions, are strongly determined by their porous microstructures. One important mechanism for producing the porous microstructures of membranes is polymerization-induced phase separation (PIPS). Here, we scrutinize PIPS by employing a Cahn–Hilliard–Navier–Stokes method coupling with the Flory–Huggins model. We focus on the formation of membranes via diffusion as well as capillary flow. We report several morphological evolution characteristics of PIPS: (1) an asynchronous effect, where the polymer-rich phase and the polymer-lean phase reach their equilibrium concentrations at different times, (2) a center-to-center movement and collision-induced collision of polymer-rich particles, (3) transition of network structures into polymer particles and rebuilding of network structures from polymer particles, (4) polymer ring patterns. We expect that these findings would shed light on complex microstructures of membranes and provide guidance for the fabrication of desired membranes.


Introduction
Membranes have been applied in a wide range of areas, such as microfiltration [1], ultrafiltration [2], fuel cells [3], reverse osmosis [4], and gas separation [5]. These membranes are mostly fabricated from polymer solutions, such as resorcinol/formaldehyde (RF) [6], cellulose [7], and poly-methacrylates-alcohols [8]. The functionality, reliability, and durability of the membranes strongly depend on the microstructural patterns. One underlying mechanism for the formation of the microstructures in membranes is phase separation [9], where the polymer solution decomposes into two immiscible phases: one with a highpolymer concentration and one with a low-polymer concentration. Shortly after the phase separation, the concentrated phase solidifies and transforms into membranes. This process is frequently denoted by the general term "gelation" [10].
The formation of membranes via phase separation has been extensively observed in experiments [11][12][13][14][15][16][17] (see Fig. 1). In order to provide additional insights into porous microstructures of membranes, a number of numerical models have been developed. For instance, dissipative particle dynamics simulations are carried out to study the structure and kinetics of membrane formation via thermally induced phase separation [18]. Monte Carlo simulations have been adopted to study phase separation dynamics in the presence of block polymer [19]. Besides, the Cahn-Hilliard model incorporating the Flory-Huggins theory [20][21][22] is an alternative and widely used method to investigate the phase separation forming membranes.
There are two aspects which have been seldom contemplated in the models in literature. The first one is that during the formation process of membranes, the degree of polymerization generally increases with time. At the initial stage, the free energy density of the polymer solution has a  [23] and a miscibility gap does not exist. Thus, the monomers are fully miscible with the solvent, e.g., resorcinol can be dissolved in water to a large extent and silicic acid is capable of forming a homogeneous solution with alcohol and water. As the monomers react with each other via a polycondensation reaction, oligomers are formed. Because of the increase in the degree of polymerization, two local minima appear in the free energy density [23] and these polymer chains are no longer miscible with the solvent. This process is termed as polymerization-induced phase separation (PIPS) [10]. PIPS has been widely investigated in literature and is considered to be one vital mechanism for the formation of membranes [24]. Arising from the cross-linking reaction of monomers leading to an increase in DP and the constant DP of the solvent DP = 1, the entropy contributions from the polymer and solvent are asymmetrical. The effect of this asymmetry on the kinetics of the microstructural evolution has not yet been fully revealed, though some literature have been denoted to this topic [25]. We will elucidate that this asymmetrical contribution leads to an asynchronous evolution of the polymer-rich and the polymer-lean phase.
The second aspect is the capillarity effect which is usually overlooked in previous models. The ligaments and the pores of the membranes are typically in the scale from nanometers to micrometers. At these length scales, capillary force has been proved to be more significant than other forces, i.e., gravity. In the present work, we present a Cahn-Hilliard-Navier-Stokes model which includes the time-dependent degree of polymerization as well as the capillary effect. By considering a time-dependent phase diagram and capillary flow, we report the following microstructural evolution characteristics. (1) The polymer-rich phase and the polymer-lean phase asynchronously reach the equilibrium state. (2) Capillary flow gives rise to a center-tocenter movement of the polymer-rich particles. (3) Capillary flow leads to a broader range distribution of the polymer-rich particles in comparison with the diffusion-controlled phase separation. Tiny and relatively large particles occur when capillary flow is present. (4) When pre-nucleation is considered, a polymer ring pattern appears. Our simulations are performed both in 2D and 3D by considering upper as well as lower critical temperature systems.
The rest of the paper is structured as follows: in Section 2, we present the phase diagram based on the Flory-Huggins theory. In Section 3, we depict the Cahn-Hilliard model coupling with the Navier-Stokes equation. In section 4, we show the simulation results in 2D and illustrate an asynchronous kinetics for the polymer-rich and polymer-lean phases. The size distribution of the particles with and without capillary flow is compared. In Section 5, we demonstrate three-dimensional microstructures of PIPS. We conclude the paper in Section 6.

Phase diagrams and the Flory-Huggins model
There is a paucity of thermodynamic databases for polymer solutions forming membranes. Even though some thermodynamic databases have been produced experimentally [26], it is still a great challenge to obtain the dynamic-phase diagram which depends on the degree of polymerization. Hence, to study PIPS in the formation process of membranes, we use the model of Flory and Huggins (FH) for a polymer mixture with an average degree of polymerization DP = N and a solvent. The basic equations can be found almost everywhere in the literature on polymer solutions [27][28][29][30][31]. The starting point is the free energy density of mixing, which reads where k b is the Boltzmann constant and T denotes the temperature. The first two terms in Eq. (1) represent the entropy change of the solution due to mixing of the polymer with volume concentration φ in the solvent. The last term is the latent heat of mixing, characterized by the Flory parameter χ ¼ ZΔϵ=k B T, in which Z is the coordination number, Δϵ is the difference between the polymer-solvent interaction energy ϵ PS and the average interaction energy of the polymer-polymer and solvent-solvent, Δϵ ¼ ϵ PS À ðϵ SS þ ϵ PP Þ=2. These interaction energies are determined by the intermolecular potential. If χ < 0, the solution of the polymer in the solvent is more favorable compared with the situation, where χ > 0, which means polymer-polymer and solvent-solvent combinations are energetically favored over polymer-solvent. A polymer solution with positive χ and sufficient high DPs gives rise to a miscibility gap, leading to the decomposition of the solution into two phases, as illustrated in Fig. 2. The reciprocal dependence of χ on temperature results in a miscibility gap with an upper critical solution temperature (UCST). This can be seen by looking at the spinodal line, which is defined by the zeros of the second derivative of the free mixing enthalpy with respect to the volume concentration, yielding The plot of 1=χ or T against the volume concentration with varying degree of polymerization reveals the effect of DP on the shape of the spinodal regions. The dependence of the spinodal region as a function of degree of polymerization is shown in Fig. 2a. Besides the spinodal, the phase equilibrium concentration is defined by the binodal lines, which are derived from a common tangent construction, namely, μ p :¼ ∂ φ f j φ¼φ p ¼ ∂ φ f j φ¼φ s ¼: μ s and ðf À μ p φÞj φ¼φ p ¼ ðf À μ s φÞj φ¼φ s in which φ p and φ s are the equilibrium polymer concentrations in the polymer-rich and polymerpoor (solvent) phases and μ s , μ p are the chemical potentials. The unknowns φ p and φ s are determined by the Newton iteration method and the plot T versus φ p and φ s gives rise to the binodal line. As an example, the binodal line for N ¼ 10 is illustrated by the dashed curve in Fig. 2a.
The critical point of the miscibility gap only depends on the degree of polymerization. This can be seen by calculating ∂χ=∂φ ¼ 0 leading to and the critical temperature (inverse of χ) varies like The above outlined Flory-Huggins model with an upper critical point could in principle work for membranes. But it might be useful to also discuss another FH-model, namely, one with a lower critical point where the phase boundaries are inverted. A phase diagram with a lower critical solution temperature (LCST) is most easily described by a Flory parameter χ that depends linearly on temperature [27] χðTÞ with χ 1 being a constant. In literature, the Flory parameter is mostly written in such a form χ ¼ A þ B=T [32], with A and B being constants. Equation (5) may be interpreted as a binomial expansion of the expression χ ¼ A þ B=T with only linear terms, so that Eq. (5) is consistent with previous formulations. Such a behavior is for instance observed in systems with hydrogen bonds [27][28][29]. The relation Eq. (5) describes the polymer-solvent interaction being more strongly dependent on temperature. Performing the same analysis gives rise to the LCST spinodal regions, as shown in Fig. 2b. One should notice that the binodal and spinodal always meet at the critical point and that the binodals are always broader than the spinodals. Hence, for a given temperature and an initial low volume concentration, e.g., φ 0 < 0:1, with an increase in DP, the system firstly touches the binodal line and thereafter passes through the spinodal line. It is also worth mentioning, that polycondensation reactions can typically be described by a second order kinetics and thus the average degree of polymerization depends linearly on time as N ¼ 1 þ kt [33,34], with k being a temperature dependent reaction rate constant. Although, as mentioned above, dynamic-phase diagrams are not really available for polymer solutions, one can make the following qualitative prediction for some polymer solutions. Figure 2c, d schematically show both systems and a state point of a resorcinol-formaldehyde mixture indicated by the black circle. If one would perform a sudden temperature change in an UCST system, the state point would move into the two-phase regime, whereas in the LCST system the state point stays in the single phase regime. Experimentally, it is exactly the latter what is observed for the RF-water system. Whenever a solution, kept at a certain temperature (room temperature or higher), is transferred to a fridge, it can stay there without gelation for quite a while (weeks to months). If the RF-water system exhibits an upper critical point, this would not be possible, since a quench to a lower temperature immediately leads to a phase separation and eventually to gelation. Therefore, from an experimental point of view, the RF-water system must be one with a lower critical point, while this is not the case for silica-water-alcohol or biopolymer-solvent systems. The following discussions will consider both UCST and LCST phase diagrams.
3 Phase-field model with fluid flow In a closed isothermal isobaric system, the microstructural evolution of the polymer solution is such as to decrease the free energy functional of the system, which is written as [35]: in which V is the domain occupied by the system and κ, a gradient energy coefficient, is determined by the interfacial tension σ as κ ¼ σ= 2 R 1 À1 ðdφ=dxÞ 2 dx . The interfacial tension σ is related to the free energy density as σ ¼ 2 dφ (see the derivation in [35]), where Δf ¼ f ðφÞ À f ðφ s Þ þ μ s φ À μ s φ s and the effect of the interaction coefficient χ on the interfacial tension is incorporated in the free energy density. The concentration evolution follows the variational approach as [36]: where MðφÞ denotes the mobility and is set as Here, D is the interdiffusivity, v m is the molar volume and N A is the Avogadro constant. In general, the diffusivity D depends on φ, which can be an interpolation of the self diffusivity of polymer D p and solvent D s according to Darken's equation (for more details see [23]). For simplicity, we set a constant diffusivity in the present work. The last term ξ in Eq. (7) represents noise fluctuations. Aiming at a more realistic scrutiny on the microstructural evolution, in which liquid phases are involved, the contributions to the reduction of the free energy are attributed not only to diffusion, but also to the fluid flow.
The mass flux generally is comprised by two synergistic contributions: the diffusion arising from the chemical potential gradient and the convection φu describing the mass transport via the mean fluid velocity u. With the generalized mass flux, J ¼ ÀMðφÞ∇ðδF =δφÞ þ φu, the mass conservation is rewritten as [37]: In the classic Cahn-Hilliard model, the kinetic energy arising from the capillary flow has not been taken into account. However, many in-depth investigations [38][39][40][41][42] have shown the importance and the necessity for the consideration of convection. In this way, the decrease rate of the free energy caused by convection is ðdF =dtÞ convection ¼ R V ðδF =δφÞð∂ t φÞ convection dv, which is equivalent to À R V ðδF =δφÞðu Á ∇φÞdv with the aid of the incompressible condition ∇ Á u ¼ 0. Thus, the soaring rate of the kinetic energy R V f s Á udv equates to the decline rate of the Gibbs free energy and the capillary force f s is formulated as f s ¼ Àφ∇ðδF =δφÞ [43]. This derivation does conform to Noether's theorem where a stress tensor is yielded [40,[43][44][45][46] in which I is the identity matrix. Furthermore, a generalized momentum balance equation [43,47] can also be obtained by including the capillary force in the Navier-Stokes equation where ρ, η and p are, respectively, the density, the viscosity, and the pressure. The density and viscosity of the solutions are interpolated over the individual species as: ρðφÞ ¼ ρ p hðφÞ þ ρ s ½1 À hðφÞ; ð11Þ where ρ p (η p ) and ρ s (η s ) are the densities (viscosities) of the polymer and solvent, respectively. hðφÞ is an interpolation function such that hðφÞj φ¼φ p ¼ 1 and hðφÞj φ¼φ s ¼ 0. Equation (8) coupling with Eq. (10) is named as the phase-field model. A crucial feature of this model is that it does obey the second law of thermodynamics as the sum of the free energy and the kinetic energy (E v ðφ; uÞ) declines with time, e.g., dE v ðφ; uÞ=dt 0 [37,48] Equations (8) and (10) are non-dimensionalized by choosing d 0 and d 2 0 =D as the dimensionless factors for the space and time, respectively, where d 0 is the characteristic length. The two evolution equations are solved by the finite difference and explicit Euler methods on a staggered mesh. Periodic boundary conditions are applied in all dimensions. The parameters for the discretization are

Simulated 2D structures
In this section, we perform simulations on a 2D domain with a size of 300 Â 300. We consider a liquid mixture initially with a polymer concentration of φ 0 ¼ 0:06 at a dimensionless temperature of k b T=ZΔϵ ¼ 1. This setup is denoted by the blue circle in Fig. 3. A Gaussian distributed noise is introduced to perturb the homogeneous liquid phase, as shown in Fig. 4a, e. For N < 10, this concentration is outside the spinodal and the disturbed liquid phase gets homogenized with time. Presently, DP is set to be 10 at the beginning such that the perturbation can induce spinodal decomposition immediately. Arising from polymerization, DP increases with time following the expression N ¼ 10 þ kt, where k ¼ 2 Â 10 8 s À1 . Because of the rising in DP, the spinodal region shifts downwards and leftwards (Fig. 2b). With time, the initial concentration relatively moves into the interior of the spinodal region with increasing driving force, which leads to a typical self-similar worm-structure (Fig. 4b). This selfsimilar structure has a rather narrow distributed polymersegment, which is a characteristic of the early stage of spinodal decomposition. This feature is inherited to the subsequent microstructures (Fig. 4c, d), which possess quasi-homogeneous polymer droplets.
Per contra, with the presence of capillary flow, the spinodal decomposition rate is significantly enhanced, as shown in Fig. 4f. It is noted that the worm-like structure also occurs in this case (not shown), but it has transformed into droplets much earlier than in the simulations without capillary flow. Moreover, the capillary flow provokes directional movements of percolated polymeric droplets (see Section 5 and [23] for more details), in contrast to the stochastic Brownian motion. This directional motion, which is not evident in Fig. 4a-d, prompts center-tocenter impingements for the polymeric droplets. As a result, the size distribution of the droplets is more inhomogeneous than the diffusion-controlled case, as can be seen in Fig. 4g, h. The fluid velocity field illustrating the motion and collision of droplets during phase separation is shown in Fig. 4i-l. It is noted that high-fluid velocities occur at the positions where collisions take place. This phenomenon is especially pronounced in Fig. 4k, where the fluid velocity in the highlighted region is much higher than the ones in other regions, as demonstrated by the color bar. Because of this large contrast of the velocities, fluid flow outside the highlighted region is not visible. Strong flows are observed in regions of collision events due to the negative curvature established at the neck when two droplets tend to impinge. The strong flow in Fig. 4k corresponds to a Péclet number around 10, which is at least one magnitude higher than the ones in Fig. 4i, j, l. Being ascribable to the large difference in the curvature between the neck and other convex surface of the droplets, a natural flow is induced to accelerate the collision process. This phenomenon is the so-called collisioninduced collision. In comparison with the Brownian motion flow estimated by the Stokes-Einstein equation D s ¼ k b T=ð6πηrÞ % 10 À10 m 2 /s, this surface tension directed flow is much stronger. Figure 5a illustrates the maximum (cyan) and minimum (orange) volume concentration as a function of time, corresponding to Fig. 4a-d. The cyan and orange dashed lines describe the time-dependent equilibrium concentration in the polymer-rich and polymer-lean phases, respectively. The dot-dashed lines show the time-elapsed spinodal concentrations. As can be seen from the solid lines, the polymer and the solvent phases start to separate from each other evidently after t ¼ 37:5 ns, resulting in the formation of polymeric droplets (Fig. 4c). A noteworthy feature is that when the concentration in the solvent phase (orange line) reaches the equilibrium value, the polymeric phase is still far from equilibrium. In order to interpret this asynchronous effect, we calculate the evolution rate of the polymer and solvent phases as a function of time based on the data in Fig. 5a and the results are shown in Fig. 5b. As we can see, the absolute value of the evolution rate of the solvent phase (orange line) increases and reaches a peak at the time around t ¼ 37:5 ns. This time is the moment when the solvent phase moves outside the spinodal region, as indicated by the crossing point of the orange solid and dot-dashed lines in Fig. 5a. Shortly after the time t ¼ 37:5 ns, the evolution rate of the solvent phase rapidly decreases and tends to zero. The evolution rate of the polymer phase as a function of time (cyan line) follows a similar behavior as the solvent phase, but is much greater than the latter one.
In order to understand these two different evolution rates, we estimate the second derivative of the free energy density d 2 f =dφ 2 ¼ 1=½NðtÞφðtÞ þ 1=½1 À φðtÞ À 2χ by using the concentration φðtÞ in Fig. 5a and the time evolution of the DP: NðtÞ ¼ 10 þ kt. The motivation for this estimation is that d 2 f dφ 2 is the thermodynamic driving force for the phase transition. Supposing that φ 0 is the initial concentration and δφ is the perturbation, the difference of the free energy density between the perturbed and initial states is The value of d 2 f dφ 2 with time is displayed in Fig. 5c. At the beginning stage, the values of d 2 f dφ 2 for the polymer (cyan) and the solvent (orange) phases both are negative, so that phase separation takes place in order to reduce the free energy. For the polymer-rich phase, its concentration gradually increases with time undergoing an "uphill" diffusion process. The summit corresponds to the peak of the free energy curve with the most negative of d 2 f dφ 2 (cyan arrows in Fig. 5c, d) and the maximum growth rate of the polymer phase in Fig. 5b.
Succeeding the cyan arrow near the cyan solid line in Fig. 5b, the subsequent three black arrows denote several abrupt increases in the maximal polymer concentration. The microstructures just before and right after these three discontinuous points are displayed in Fig. 6a-f, respectively. A heedful scrutiny of these microstructures shows that these discontinuous points are due to the fact that the maximal concentration shifts from one droplet to a different one, as highlighted by the squares. For instance, the droplets with  Fig. 6a, b are distinct, which leads to the first discontinuous point in Fig. 5b. The change of the maximal concentration from one droplet to another one is as a result of distinct Ostwald ripening rates, which depend on the size and quantity of the neighboring droplets. From 39.925 to 45.8 ns, the droplets owing the maximal concentration are identical (Fig. 6b, c), so that at this stage, the curve in Fig. 5b is smooth.
For the solvent-rich phase, its concentration evolves towards the left spinodal point (orange arrow in Fig. 5d) through "downhill" diffusion, so that in Fig. 5c, the solvent phase (orange line) does not pass across a maximum separation rate with a lowest negative of d 2 f dφ 2 before achieving the left spinodal point. After going through the spinodal point, d 2 f dφ 2 becomes positive and significantly increases with a slight decrease in the concentration, arising from the fact that the free energy curve is extremely compact at the solvent side. Such a dramatic increase in d 2 f dφ 2 gives rise to a maximum evolution rate of the solvent phase, as represented by the peak of the orange curve in Fig. 5b. It is emphasized that this peak for the evolution rate of the solvent phase differs from that of the polymer phase which is engendered from a maximum negative of d 2 f dφ 2 . As can be seen in Fig. 5c, after the spinodal, d 2 f dφ 2 of the solvent phase is remarkably larger than the one of the polymer phase. Due to such a significant difference in the thermodynamic driving The microstructures before and after the three discontinuous points in Fig. 5b force, the solvent phase attains the equilibrium value much earlier than the polymer phase. This asynchronous effect is clearly visible in Fig. 5b that when dφ=dt for the solvent phase is completely zero, dφ=dt for the polymer phase is still above the horizontal dashed line. Because the difference in the driving force d 2 f dφ 2 for the polymer and the solvent phases increases with DP, it is expected that the asynchronous evolution is more pronounced with larger DPs. The analysis for the evolution rates in Fig. 4e-h is skipped, since the behaviors are quite similar except that both phases reach the equilibrium concentration earlier than the one in Fig. 4a-d. Figure 7a shows the radius (R p ) distribution of the polymer particles for the microstructures in Fig. 4d, h both at the time t ¼ 47:5 ns. In the case without convection, the radius of the particles locates in a relatively narrow range around R p ¼ 3 nm. In contrast, with convection, the distribution is much broader and there are two peaks in the histogram, one at R p ¼ 3 nm and the other one at R p ¼ 5 nm. In addition, as a result of the motion and collision mechanism, the number of the total polymer particles is significantly reduced. In order to reveal more details about the effect of the convection on the particle size distribution, we plot the particle radius distribution with and without convection with the same maximum volume concentration φ m . As illustrated Fig. 7b, c, φ m ¼ 0:5, 0.6 and 0.69 are depicted by the blue, red, and green histograms, respectively. In both cases with and without convection, the range of the distribution is nearly the same. Without convection, from φ m ¼ 0:5 to φ m ¼ 0:6, spinodal decomposition is the dominating evolution mechanism, so that an evident peak occurs at R p % 2:6 nm in the red histogram. Thereafter, Ostwald ripening plays a more important role than the spinodal mechanism. As a consequence, particles with R p > 2:6 nm are gradually formed at the expense of small particles with R p < 1:5 nm. However, in the case with convection, in the whole evolution range from φ m ¼ 0:5 to φ m ¼ 0:68, spinodal is always the decisive mechanism for the microstructural evolution. The peak at R p % 2:6 nm consecutively rises up, which is a typical characteristic of the spinodal decomposition. Although large particles with R p > 3 nm are observed, the underlying formation mechanism is caused by convection-induced collision, rather than Oswald ripening in Fig. 7b. This is evidenced as follows: (1) Small particles with R p 2 ð1:0; 2:2Þ nm are not obviously consumed from the red to the green histograms in Fig. 7c. (2) The big particles with R p > 2:6 nm exhibit a smooth distribution in the green histogram of Fig. 7b, since these large particles are transformed continuously from the small ones by diffusion. However, in Fig. 7c, the distribution of large particles with R p > 2:6 nm is relatively dispersed because of occasional collisions.
In order to further understand the microstructural evolution with convection, we analyze the quantity N of the polymer particles as a function of time in the time interval t 2 ½30; 60 ns and the results are shown in Fig. 8. As illustrated in Fig. 8j, the evolution is divided into three stages. From 30 to 33 ns (violet region), the number of the polymer particles rapidly increases with time. At this stage, the concentration in the particles constantly increases with time because of spinodal decomposition, leading to the growth of the particles. The size of these particles is relatively small (see Fig. 7c). A few large particles with radius > 3 nm are obtained due to occasional collisions, as highlighted by the white circles in Fig. 8b (before collision) and Fig. 8c (after collision). In the second stage (cyan region), the concentration in the polymer particles does not increase significantly and the quantity of the particles is reduced dramatically engendering from frequent collisions. The particles before and after collisions are highlighted by different colored circles in Fig. 8d, e, respectively. In the third stage, collision still takes place, which results in a further decrease in the total number of the particles, but the impingement and clash events are not as often as that in the second stage due to limited quantity. Hence, the large particles in the third stage mostly stem from the consecutive collisions in the second stage. For instance, the large particles with a radius of around 5 nm in Fig. 7a are produced as a consequence of the collision of two small particles with a radius of about 3 nm shown in Fig. 7c. After the collision, the volume of the large particles is slightly greater than the summation of two small particles, since due to the polymerization, the large particle is still consuming the polymer species in the matrix to increase its volume. After the third stage, collision continues until gelation. It is noted that in contrast to the Ostwald ripening mechanism for the formation large particles in the case only with diffusion, collision is the main reason for the generation of large particles when convection is present. For the Ostwald ripening mechanism, the mass flux is perceptually from the small particle to the large one via diffusion. When convection occurs, an additional flux uφ is added to the diffusion flux. This convection flux sometimes works against the diffusion, leading to the drift of the particles. More details are discussed in Fig. 13.

Three-dimensional studies
In this Section, we focus on phase separation in threedimensions with different initial concentrations for an UCSD system. A dimensionless parameter C :¼ R g Td 2 0 =ðρD 2 v m Þ is introduced to characterize the intensity of the capillary flow. The simulations are performed on a cubic domain with a size of 400 Â 400 Â 400. Figure 9 shows the spinodal lines for a series of DPs for an UCSD system with a Flory parameter χ ¼ 1:5. The horizontal dashed line represents the storage temperature of the sample T g ¼ 1:6 and the cyan circles along the dashed line denote four different initial concentrations φ 0 ¼ 0:05, 0.1, 0.2, and 0.3 for the simulations. For N ≲ 7 (red line), the storage temperature is greater than the critical temperature and the polymer solution remains a homogeneous mixture. With an increase in DP, the polymer solution with φ 0 ¼ 0:3 firstly goes inside the spinodal region. The other The horizontal dashed line denotes the storage temperature T g ¼ 1:6 and the circles represent the initial polymer concentrations for the three-dimensional studies three polymer-solutions will also be inside the spinodal region when DP is sufficient large. The cyan dashed line depicts the binodal line for N ¼ 20. Figure 10 illustrates the time evolution of the microstructures for three different setups: (1) only with diffusion, (2) with diffusion and weak capillary flow (C ¼ 10), and (3) with diffusion and strong capillary flow (C ¼ 1000). At the beginning, the concentration is uniformly set as φ 0 ¼ 0:3 at a dimensionless temperature T g ¼ 1:6 and with a DP of N ¼ 1. A Gaussian distributed noise is imposed to the homogeneous solution. The initial states of the three cases (1)-(3) are shown in Fig. 10a, d, g. If DP is kept at N ¼ 1, the noise will be smoothened out by diffusion and a homogeneous polymer solution is achieved. However, due to the polymerization, i.e., the increase of DP, the polymer solution goes into the spinodal region, as depicted in Fig. 9. Thereafter, the system decomposes into two separate phases. The interface of the polymer-rich and the polymerpoor phases is represented by the cyan color in Fig. 10; the yellow surface is an isosurface standing inside the polymerrich phase. Here, the interface is defined by the condition are, respectively, the minimum and maximum concentrations in the whole domain. The microstructure further evolves with time and gets coarser (see Fig. 10c, f, i), caused by surface area minimization. In the case (1) only with diffusion, a bicontinuous cluster structure is formed. When a weak capillary flow is considered (case (2)), the ligaments of the bi-continuous cluster structure are thicker than the previous one. In the case (3) with a stronger capillary flow, tiny polymer-particles (cyan) as well as small pores (yellow) both are observed besides polymer-ligaments. Figure 10j illustrates the volume fraction of the polymerrich phase v p (solid line) and the permeability of the  (1)), d-f with diffusion and weak capillary flow (C ¼ 10, case (2)), and g-i with diffusion and strong capillary flow (C ¼ 1000, case (3)). The cyan surface denotes the interface of the polymer-rich and the polymer-poor phases. The yellow one illustrates an isosurface inside the polymerrich phase. j The volume fraction of the polymer-rich phase v p and the permeability K as a function of time for the case (3). The diagram is divided into two stages: I and II. k Sketch of the dynamic-phase diagram in the stages I and II (Color figure online) microstructure K (dashed line) as a function of time for the case (3), corresponding to the microstructures in Fig. 10g-i. The volume fraction increases (stage I), followed by a decrease (stage II) with time. This observation is in contrast to the conventional spinodal decomposition, where the volume fraction monotonically converges to the equilibrium value. The nonmonotonic temporal behavior of the volume fraction can be explained by using the dynamic T-φ phase digram, as sketched in Fig. 10k. The red circle depicts the initial concentration φ ¼ φ 0 , while the cyan and orange circles represent the equilibrium concentrations in the stages I and II, respectively. The horizontal dashed line denotes the storage temperature T g . In the stage I, the equilibrium volume fraction of the polymer-rich phase is given by > 0:5 according to the lever rule. The initial noise gives rise to a volume fraction v p ¼ 0:5, which is fair for both phases, as can be seen in Fig. 10j. In order to reach the temporal equilibrium v I p , v p has to increase in the stage I. Due to the change of DP, the phase diagram shifts upward and leftward, as drawn by the orange line in Fig. 10k. As a consequence, the volume fraction tends to fulfill the following expression, v II p ¼ in the stage II. Since the shift of the equilibrium concentration in the polymer side is greater than the one at the solvent side, i.e., φ II p À φ I p > φ I s À φ II s , the volume fraction v II p is less than v I p . Therefore, v p has to decrease in the stage II. Similar nonmonotonic evolutions of the volume fraction (not shown) are also observed in the cases (1) and (2). The permeability of the microstructures is measured at five different time steps, as shown by the square symbols in Fig. 10j. For the permeability measurement, fluid flow simulations are conducted by solving the Stokes equation with the assumption of low Reynolds number, i.e., Re ( 1. By imposing a pressure gradient at the boundaries, the resulting velocity field in the pore space is used for the calculation of the permeability by applying Darcy's law (see [49] for more details). In the stage I, more closed pores are developed in the microstructure with time due to capillary flow (see Fig. 10h). These closed pores reduce the remaining flow paths inherently and therefore the permeability decreases. In the stage II, the closed pores move in space and interconnect with each other. As a result, the amount of open-pores increases and the classical Kozeny-Carman (K-C) equation [50] can be applied, , where a is the K-C coefficient and S s denotes the specific surface area. Since v p and S s both decrease with time in the stage II, the permeability rises. As observed, interrupting the process of phase separation at different points in time results in varying patterns and, accordingly, permeability is changing. Furthermore, two different types  Fig. 11. For φ 0 ¼ 0:05 (Fig. 11a), PISD leads to the formation of a number of polymer droplets in the case (1) and some droplets are connected establishing a small cluster. With an increase in C (Fig. 11b, c), a number of relatively large polymer droplets appear and the amount of the connected droplets are less than the ones in the case (1), as can be seen in the insets of Fig. 11a-c. For an initial concentration φ 0 ¼ 0:1, PISD gives rise to the formation of a cluster structure in the case (1), as shown in Fig. 11d. The cluster structure transforms into a structure mixed with clusters and droplets in the case (2) (Fig. 11e) and completely into a droplet-structure in the case (3) (Fig. 11f). For an initial concentration φ 0 ¼ 0:2, cluster structures are also observed in the cases (1) (Fig. 11g) and (2) (Fig. 11h). The microstructure in the latter case is coarser than the one in the former case. Apart from clusters, a number of very tiny droplets appear in the case of a strong capillary flow, as shown in Fig. 11i. It is noteworthy that tiny droplets are also observed in Fig. 11f. We stress that the formation of the tiny droplets can only be observed when the polymerization is coupled with the capillary flow. The detailed reason is reported in the following.
In order to figure out the reasons for the distinct microstructures in the cases (1), (2), and (3), a heedful look has been taken to scrutinize the evolution process of PISD. Figure 12a-d illustrates the microstructures at four different time steps with an initial concentration φ 0 ¼ 0:1 and C ¼ 1000. The blue-yellow curves are the stream lines of the capillary flow. Figure 12e-h corresponds to a subregion in Fig. 12a-d. Figure 12i-l is sectional view of the microstructures in Fig. 12e-h. The red-blue arrows depict the strength and the direction of the capillary flow. The value of the flow velocity is given by the color legends at the right side. By tracing these two droplets in Fig. 12i-l, it is observed that the two droplets moves in a particular direction towards each other. In contrast to the random Brownian motion, this kind of directional approach has also been observed in other simulations with different initial concentrations in the presence of capillary flow. This directional movement is driven by the nonuniform concentration around the surface of the polymer particles, which emanates from the interactions between the particles. For instance, when two particles are different in size in proximity, the polymer solute is transfered from the small particle to the big one passing through the solution due to the difference in curvature (Ostwald ripening). The diffusion flux around the particles is inhomogeneous since the length of the diffusion path is distinct, i.e., the diffusion path through the gap between the particles is shorter than the one outside the gap. The inhomogeneity of diffusion flux engenders a concentration gradient around the particles. This concentration gradient would be dispersed by diffusion, if there were no convection. However, when convection is more pronounced than diffusion, the established concentration gradient gives rise to a thermodynamical driving force for the motion. In a two-particle system, the concentration gradient as well as the induced motion can be analytically solved in the bipolar coordinate, as initially done by Golovin et al. [51] and extended in our previous works [38,47]. In a multiparticle system, the concentration gradient can only be numerically simulated rather than analytically computed. It is noted that this directional motion is only significant when the particles are in proximity where the interactions between them are of relevance. In a dilute system where the particles are far from each other, Brownian motion is a major effect in comparison with the capillary-driven motion. Figure 13 presents some details for a better understanding on the PISD microstructures. Figure 13a-c, d-f illustrate two transition phenomena: "percolation-to-cluster" (PTC) and "cluster-to-percolation" (CTP), respectively. In the former case, clusters transform into discrete polymer droplets. This transformation typically takes place at the early stage of PISD for φ 0 ¼ 0:05 and φ 0 ¼ 0:1 either with or without capillary flow. The initial concentration noises lead to the formation of a cluster structure with a volume fraction of 50%. Subsequently, the volume fraction of the cluster structure has to decrease with time because of the lever rule. The decrease of the volume fraction initiates the PTC phenomenon. The droplet-structures in Fig. 11a, b, e and at the time before Fig. 11d (not shown) are originated from the PTC mechanism. After PTC, the resulting polymer particles grow with time because of the two following reasons. The first one is the Ostwald ripening mechanism where the big particles increase their sizes at the expense of the small particles. The second reason is that due to the highly asymmetrical phase diagram (see Fig. 9), the polymer-rich and the solvent-rich phases asynchronously arrive at their equilibrium values. The solvent-rich phase reaches the equilibrium earlier than the polymer-rich phase and a chemical potential gradient between the two phases is established. The resulting chemical potential is an additional driving force for the growth of the polymer-rich particles. When the growing particles are sufficient large, they get in touch with each other, giving rise to CTP, as shown in Fig. 13d-f. The CTP transition is only observed in the case (1) (only diffusion) or case (2) (weak capillarity) for φ 0 ¼ 0:05 and φ 0 ¼ 0:1, and has been directly captured in experiments [52]. In the case (3) with a strong capillary flow, the polymer particles arising from PTC move in space, as discussed in Fig. 12, and coalesce into relatively large particles. Thus, CTP is unlikely to occur in the case (3), where convection is dominated.
Instead of PTC and CTP in the cases (1) and (2), it is particularly noted that a strong capillary flow facilitates the appearance of little bitty polymer droplets (Fig. 11f, i) or pores (Fig. 10i) in the microstructures, which is manipulated by two different mechanisms. The first one is based on the Plateau-Rayleigh (PR) theory, as shown in Fig. 13g-i, where a fluid jet subjected to perturbations can spontaneously decompose into several droplets in order to minimize the surface area. According to the PR theory, a fluid jet evolves into a droplet only for a long wavelength perturbation, i.e., the perturbation wavelength is greater than the perimeter of the jet. A long wavelength perturbation is more prone to appear in a large scale simulation rather than in small domain simulations. Thus, the revelation of the tiny pore shown in Fig. 13g-i benefits from the present large scale simulations [30,31].
A second mechanism for the formation of tiny droplets is pictured in Fig. 13j-l, where a small droplet spontaneously appears in between the big droplets and grows with time. The explanation for this phenomenon is as following: Without capillary flow, diffusion is responsible for the mass transfer from small droplets to large droplets according to the mechanism of Ostwald ripening. The diffusion fluxes are sketched by the black arrows in Fig. 13m, where the mass flows from the droplet 1 to droplets 2-5. The diffusion flux between droplets 1 and 2 is greater than the one between droplets 1 and 3, since the larger curvature difference brings a larger diffusive flux. The difference in the diffusion fluxes further leads to a concentration gradient in Fig. 13 Some characteristic PISD microstructures: a percolation-tocluster (PCT) (a-c) and cluster-to-percolation (CTP) (d-f) transition phenomena. These two transitions are observed in the simulation with φ 0 ¼ 0:05 and C ¼ 0. g-i The formation of a tiny pore based on the Plateau-Rayleigh theory. This observation is taken from the simulation with φ 0 ¼ 0:3 and C ¼ 1000. j-l The appearance and the growth of a tiny droplet due to the presence of capillary flow. This phenomenon is taken from the simulation with φ 0 ¼ 0:1 and C ¼ 1000 and is in contrast to the Ostwald ripening mechanism, where small droplets are always consumed by big droplets. m, n Sketches for the reason why a small droplet can survive and grow when capillary flow is present the matrix. Even though the radii of droplets 2 and 3 are equal, a concentration gradient in the matrix can be established if the distances between 1 and 2 and between 1 and 3 are different. The resulting concentration gradient subsequently induces a convection flux, as illustrated by the red arrows in Fig. 13n. When the net flux at one spatial point in the matrix is positive, a polymer nucleus appears. Because of the increase of DP and the resulting asymmetrical phase diagram, when the matrix has already reached the equilibrium concentration, the newly formed polymer nucleus is still far from the equilibrium, as aforesaid. As a consequence, the new polymer nucleus grows with time caused by the chemical potential difference between the polymer nucleus and the matrix. This abnormal growth of tiny droplets is in contrast to the quintessential Ostwald ripening mechanism, where big droplets always consume little droplets. The mechanism for the formation of the small pores in Fig. 10i is the same as that of the small droplets. The difference is that in the former case, the pore is the minority phase, whereas in the latter case, the polymer particle is the minority phase.
As aforementioned, polymerization shifts upward and leftward the spinodal region, which results in a decomposition for the polymer solution with a concentration initially outside the spinodal region. However, prior to the spinodal decomposition, the system may pass through the binodal region wherein the classic nucleation could take place. In order to shed light on the effect of the classic nucleation on the microstructural evolution, we have carried out two simulations, one without pre-existing nuclei and the other one with pre-existing nuclei. The time evolution of these two simulations are shown in Fig. 14. In the former case, polymerization gives rise to the phase separation of the polymer solution. The positions for the formation of the polymer-rich phase depend on the random concentration noise imposed initially. In the latter case, a number of polymer-rich nuclei (Fig. 14k, p) pervade the domain at the initial stage. Arising from the polymerization, these nuclei grow with time caused by the thermodynamical driving force of binodal decomposition, as can be seen in Fig. 14l, q. With a further increase in time, the system goes inside the spinodal region and the polymer solution (dark blue phase) is energetically unstable. It is observed that in this case with pre-existing nuclei, the phase separation takes place symmetrically around the growing polymer nuclei, which yields a ring pattern, as shown in Fig. 14m, r. When the polymer nuclei grow, these rings expand synchronously until they get in touch with each other and thereafter transform into droplet-structures by coalescence. The microstructures right after the coalescence of the ring patterns and after a certain long time are depicted in Fig. 14n, o, respectively. As highlighted by the yellow square in Fig. 14s, the PISD microstructure with pre-existing nuclei is more regular than the one without pre-existing nuclei. Simulations with capillary flow are also performed for the second case and similar ring patterns are achieved. These polymer ring patterns are attributed to PISD as well as the binodal decomposition.
A typical distribution of the concentration around a preexisting nucleus at different time steps is shown in Fig. 14u. The space position for the concentration is sketched by the orange line in the inset. As shown by the violet line, the growing nucleus perturbs the energetically unstable liquid matrix, causing the occurrence of a highconcentration region (X ≳ 15) in front of the nucleus. Thereafter, there are two competing effects for the concentration evolution. One is a further increase in the concentration in the high-concentration region, as depicted by the green, blue, and orange lines. The driving force for this evolution is spinodal decomposition. The other one is the mass transport from the high-concentration region to the low-concentration region between the ring and the nucleus. The driving force for this evolution is binodal decomposition. A consequent result of this evolution is that the ring is pushed outward and the area of the low concentration increases, which can be seen from a comparison for the concentration curves at times t ¼ 2:7 ns and t ¼ 6:4 ns. With time, the ring impinges with other neighboring rings and breaks into pieces of irregular particles (see, e.g., Fig. 14s). The spheroidization of these irregular particles leads to an occupancy of the low-concentration region by the relatively high-concentration irregular particles, so that the low-concentration region shrinks with time, as described by the blue and orange lines. Meanwhile, the concentration in these irregular particles increases with time arising from spinodal decomposition. This process is until that the nucleus connects with these particles, forming an irregular cluster structure (Fig. 14t), which is in contrast to the structure with spherical particles in Fig. 14j.
The porosity as a function of time for these two cases with and without pre-existing nuclei is illustrated in Fig. 14v. At the same time, the porosity in the former case (green line) is less than the one in the latter case (violet line). There are two main reasons for this observation. (1) Because of the supersaturated matrix, the pre-existing nuclei grow with time. This gives rise to a higher volume fraction of the polymer-rich phase and hence a lower porosity. It is noteworthy that the initial setups with a small amount of pre-existing nuclei do not have an evident impact on the porosity at the time t ¼ 0, as can be seen in Fig. 14v. (2) The second reason is that phase separation in the case with pre-existing nuclei occurs earlier than the scenario with random noise. This is explained as follows. The pre-existing nuclei grow with time, perturbing the liquid matrix. The wavelength of this perturbation is greater than the critical wavelength of spinodal decomposition. In this case, around all the growing nuclei, phase separation occurs, as can be seen in Fig. 14r. In contrast, the random noise contains perturbations with wavelength less and greater than the critical wavelength of spinodal decomposition. In this scenario, different wavelengths compete with each other and eventually, only those wavelengths greater than the critical one lead to phase separation. Such a competition delays the spinodal decomposition, as can be seen from the comparison between Fig. 14c, m. Due to the delay in the phase separation, at the same time, the volume fraction f-j and p-t are sectional views of a-e and k-o, respectively. In the former case, the formation sites of the polymer-rich phase are determined by the random concentration noise imposed initially. In the latter case, the pre-existing polymer-nuclei grow with time arising from the thermodynamical driving force of binodal decomposition.
Thereafter, PISD takes place around the pre-existing polymer-nuclei and a ring pattern is formed. With time, the polymer rings get in touch with each other and transfer into a droplet-structure by coalescence. u A typical concentration distribution at different time steps around a pre-existing nucleus. The inset sketches the space position of the concentration. v The time evolution of the porosity for the microstructures with and without pre-existing nuclei of the polymer-rich phase in the latter case is less than the one in the former case.

Conclusion and outlook
In conclusion, we have investigated PIPS that is one important mechanism for the formation of porous microstructures in membranes. Our work is based on the Flory-Huggins theory, which takes the degree of polymerization into account for the free energy mixing. Since the increase in the degree of polymerization with time, the free energy and the phase digram become highly asymmetric. When the solvent phase has already moved outside the spinodal region and attained the equilibrium concentration, the polymer phase is still inside the spinodal region. This observation reveals different diffusion mechanisms for the solvent and the polymer phases. That is, when the polymer phase is governed by an abnormal diffusion mechanism with negative d 2 f dφ 2 , the solvent phase is controlled by a normal diffusion with positive d 2 f dφ 2 . Since the driving force for the concentration evolution in the solvent phase after passing through the spinodal point is much greater than the one in the polymer phase, these two phases reach the equilibrium concentration asynchronously, which leads to an asynchronous evolution of the polymer and solvent phases. This asynchronous effect gives rise to the morphological transition from continuous networks into dispersed droplets.
By analyzing the size distribution of the polymer particles without and with convection, we found different mechanisms for the formation of large particles. In the former case, the microstructure is coarsened by Ostwald ripening. In the latter scenario, convection-induced collision is responsible for the appearance of large particles.
We presently consider the polymerization-induced demixing and the system in the liquid state. These liquid structure consequently solidifies into gels. A future work is to simulate the gelation process for polymer droplets. As mentioned above, the polymer-rich droplets emanating from binodal or spinodal decomposition swim in the liquid solution and move due to the coaction of Brownian and Marangoni mechanisms. As a result of the motion, collision and coalescence, one would obtain relatively large polymer particles with a lower surface area. However, the evolution picture with the presence of gelation is totally different. As demonstrated in Fig. 15a, initially a polymer droplet arising from phase separation contains oligomers of different degree of polymerization (here we took as an example RF-oligomers indicated by the hexagons for the resorcinol ring connected by methylene bridges). Looking at the phase diagram, one typically expects a polymer fraction inside the droplets of 50% or larger. Once the percolation threshold (equilibrium concentration in the droplets) is reached, the polymer droplets undergo a gelation transition, as sketched by the red chained polymer in Fig. 15b. These gelatinizing droplets own solid-like properties, such as elasticity, which prevents from complete coalescence between neighboring polymer droplets. They can build with other solid-like particles a network on coalescence. During the gelation process, which requires a diffusional rearrangement of oligomers inside the droplets, also new polymer-rich droplets are continuously produced from the matrix by spinodal decomposition. These droplets drift in the matrix and impinge with the solid-like particles simultaneously with gelation when reaching the percolation threshold. This process of droplet formations followed by movement and gelation iterates and by this means the previously established solid-like networks are reinforced. These networks store elastic energies, which shall be considered in future works as well as the kinetics of the gelation transition inside the droplets.
As previous works [20], we assume a constant viscosity in our model for the hydrodynamic simulations. The change of the viscosity with time is important for the phase separation in polymer solutions, especially in the gelation process, and will be addressed in the future. Fig. 15 Sketch of the gelation process inside a droplet emanating from the phase separation. The gray hexagons inside a droplet schematically depict the resorcinol rings connected at the 2-position via a methylene bridges to oligomers of different degree of polymerization. After formation of liquid like droplets, the oligomers rearrange and form clusters inside the droplet since the percolation threshold is reached. The red chain in the left picture schematically shows such a cluster in a gelled droplet having then solid-like properties (Color figure online) Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.