Unconditionally Stable Algorithm for Copolymer and Copolymer-Solvent Systems

In the time evolution simulation of a copolymer system towards its equilibrium configuration, it is common to use the Otha-Kawasaki approach for free energy and time evolution by means of a Cahn-Hilliard diffusion equation. The conventional numerical resolution is to use the cell dynamics simulation method (CDS). Although this method gives an adequate response, it is limited since it needs very small time steps to present both appropriate resolution and stability. Unconditionally stable methods have recently been used in gradient systems that provide adequate resolution and stability with a greater time step in solving Cahn-Hilliard equations. In this paper we develop and implement unconditionally stable algorithms for copolymer-solvent systems and for resolution of the time evolution of block copolymer systems under the Otha-Kawasaki functional.


I. Introduction
The self-assembling ability of copolymers on a nanometric scale has aroused great interest in the development of potential applications [1]. Bulk morphologies of diblock copolymers have been extensively studied both experimentally and theoretically [2][3][4]. Depending on the composition of the blocks that comprise it, and the χN parameter, where χ is the Flory-Huggins parameter and N is the polymerization rate, the bulk structures correspond to lamellae, gyroids, cylinder hexagonal arrangements and BCC structures of spheres [5]. Recent research has shown that confinement is a powerful tool for breaking the symmetry of a structure and obtaining new shapes, different from those ob-tained in bulk [6,7]. New studies have shown that 2D and 3D confinement induces new morphologies; for instance, copolymers immersed in cylindrical or spherical nanopores. This is produced by the combined effect of confinement and curvature [8]. The new morphologies enhance new applications, such as the development of membranes in fuel cells, membranes for virus filtration, photonic materials and drug dispensers, among others [9][10][11][12]. Taking into account, firstly, the experimental difficulties of analyzing the effects of confinement and curvature on morphology, and secondly, the success of theoretical models, it is very important to have efficient calculation tools to explore the complex configurational map of confined copolymers. Simulation of the formation processes of confined copolymer structures involves a high computational cost. Due to this, it is essential to have a highly efficient algorithm to perform numerical simulations. In copolymer systems modelled under the Ginzburg-Landau theory, time evolution is carried out through the Cahn-Hilliard equation [13]. The traditional method for numerical simulation is the cell dynamics simulation method (CDS) [14]. In this article, we will discuss the advantages and disadvantages of this method and compare it with a new numerical scheme developed by Eyre [15] for the study of gradient systems. We will also develop the implementation of said numerical scheme for solvent-copolymer and copolymer systems.

II. Cahn-Hilliard equation
The Cahn-Hilliard equation was developed in 1958 to model the phase separation process of a binary mixture [16]. This approach has been extended to many other branches of science as dissimilar as polymer systems [17], population growth [18], image processing [19], and structures formed at the bottom of a river [20], among others. The Cahn-Hilliard model is a minimalist approach that contains the elements that are essential to description of the dynamics and the equilibrium properties of a wide variety of systems. The general scheme of the Cahn-Hilliard equation is the following: Where µ represents the chemical potential and F the free energy functional that describes the system under study. Particularly for a binary system, the expression of free energy is given by: Here, φ is the order parameter representing the local density of the components. Phase separation occurs during the cooling of the system, from a disorderly, high-temperature state, to an ordered lowtemperature state consisting of two phases. The dynamics are controlled by a specific length L, which increases with a power law in time L(t) ∼ t (1/3) [20]. This Law of Growth implies an extremely slow speed at advanced times during phase separation, with a speed v ∼ ∂L ∂t ∼ t (−2/3) in the domains. There is no analytical solution of the Cahn-Hilliard equation for all time, thus numerical methods are required for resolution. A wide variety of algorithms have been developed to solve the Cahn-Hilliard equation numerically. An extensive analysis of the different methods developed can be found in the reference section [22]. The traditional method is the cell dynamics method (CDS), whose advantages and disadvantages are discussed in the next section.

III. CDS method
The standard method for numerical resolution of the Cahn-Hilliard equation is the cell dynamics simulation (CDS) method [23]. This method corresponds to a Euler numerical scheme. The convergence of the Euler algorithm strongly depends on the time discretization employed in ∆t. Above a certain value of ∆t, the system is numerically unstable, presenting "chess-board" instabilities. Consequently, the CDS algorithm is subject to the ∆t max ∼ ∆x 4 [23] restriction, where ∆x is the spatial discretization. For ∆t 1 values, the CDS method provides good convergence, although it requires excessive computational time. On the other hand, in the case of Cahn-Hilliard n = 1/3, coarsening systems present evolution of their characteristic length L(t) and shape L(t) t n . This growth characteristic has a very slow domain growth speed associated with long times. It would be ideal to have a numerical scheme that provided good convergence and allowed varyiation of the time steps, adapting it to the coarsening process. As we saw in the previous equation, the CDS method cannot fulfil this requirement since ∆t is limited.

IV. Eyre algorithm
The algorithm developed by Eyre for gradient systems has the advantage of being unconditionally stable for every ∆t [15]. The algorithm guarantees stability if the free energy functional F is divided into its contractive F C and expansive F E parts. The contractive term is evaluated implicitly and the expansive term, explicitly. The numerical scheme for resolution of the Cahn-Hilliard equation is the following: Here, the C and E superscripts indicate division of the free energy into contractive and expansive. The stability condition is guaranteed if max({λ E }) ≤ 1 2 min({λ}) [24] is fulfilled, where λ E and λ are the eigenvalues of the Hessian matrices of the expansive part of the energy and of the total energy, respectively. Implementation of the Eyre algorithm for a specific system requires evaluation of the free energy that characterizes the system.

V. Eyre algorithm applications to block-copolymer systems
In this section, we will develop the Eyre numerical scheme for the Ginzburg-Landau free energy functional of a diblock copolymer system. In a diblock copolymer system the expression of free energy is the following [25,26]: Where φ represents average composition. Green's function G(r, r ) verifies the condition −∇ 2 G = δ(r − r ) and W (φ) represents the double-well potential.
The mesoscopic order is the result of competition between a short-range attractive interaction, corresponding to the gradient term of the equation, and a long-range repulsive interaction, corresponding to the long-range term that is imposed to avoid phase macroseparation, and to set the copolymer periodicity. The characteristic wavelength is λ = ( 1 α ) 1/4 . After dividing the energy into its contractive and expansive terms, in accordance with the division proposed in the previous section, and evaluating the expansive part explicitly and the contractive part, implicitly, we obtain: The stability condition for this system results in:

VI. CDS vs. Eyre comparison
The CDS algorithm registers good convergence in the numerical resolution of the Cahn-Hilliard equation. However, its effectiveness is limited by requiring an extremely small time step ∆t. In this section, the CDS and Eyre algorithms for numerical resolution of the Cahn-Hilliard equation for modeled diblock copolymer systems are compared with the Ginzburg-Landau free energy functional system, according to the result obtained in the previous section. The time evolution of a 2D system with a 1024x1024 grid with spatial discretization of ∆x = 0.5 periodic boundary conditions was simulated. A pseudo-spectral method was used for spatial derivatives [27]. For the copolymer, g = 1, α = 1, φ = 0.1 and τ = 2.1 were used. The bulk equilibrium structure corresponds to a BCC formation of spheres, while in 2D systems it corresponds to an hexagonal arrangement, as can be seen in Fig. 1, where the time evolution of the system is illustrated, starting from a disordered state. The value ∆t = 0.03 was assigned to the CDS method. The evolution of the system and the growth of the size of the domains was measured using standard image processing techniques to identify the center of each disk in the hexagonal pattern. A Delaunay triangulation [28] was then used to determine the inter-bond orientation between near-neighbor disks θ(r). Finally, the average domain size L(t) was obtained through the orientational correlation length 1 . The results are presented in Fig. 2 obtaining a growth law L ∼ t 0.33±0.01 in good agreement with the previous results. In order to compare these methods, the difference between the numerical solutions obtained by each algorithm for the above-mentioned system was calculated. The solution of the CDS method was taken as a reference. The energy of the system was calculated at each time step. Remember that in gradient systems the energy decreases as time progresses. The solutions were compared to equal values of free energy, and the E error associated with the Eyre method was calculated by expressing [30]: by the Eyre method. The notation <> indicates an average over all points of the spatial grid. The results for different ∆t values of the Eyre algorithm are presented in Fig. 3, where the error value for a given ∆t corresponds to an average over 100 comparisons of the Eyre solution with regard to the solution taken as a reference to equal energy values. For a three percent bounded error, the Eyre algorithm was 120 times faster than the CDS method.

VII. Free energy of a copolymersolvent system
We can express the free energy of a copolymersolvent ternary system based on the order parameters η and φ, where φ represents the composition of the copolymer mixture, and η, the solvent [31,32]. The energy is composed of a shortrange term F s [φ, η] and a long-range term F l [φ] The short-range term is the following: Where D 1 and D 2 are phenomenological parameters related to the interface energy. The expressions f η (η), f φ (φ) refer to double-well potentials, explicitly: Where c i and u i are phenomenological constants. The expression W (φ, η) indicates the interaction potential between the two order parameters. It is given by: The long-range term is the following: The expression above is related to the φ parameter that represents the copolymer. The b i parameters of the interaction potential characterize the interactions between block copolymers and the solvent. The evolution of the order parameters φ and η represents the evolution of the copolymer and the solvent, respectively. Time evolution results in a set of coupled differential equations for conserved order parameters.
Where F = F s + F l is the full energy.

VIII. Eyre algorithm for the Cahn-Hilliard equation
The Eyre algorithm previously developed for a copolymer system can be adapted to solve the Cahn-Hilliard equation set that models the time evolution of the copolymer and solvent equations 15 and 16. The expression of the algorithm is the following: With the following stability conditions:

IX. Numerical simulations
In order to evaluate the numerical resolution obtained by the Eyre method, the time evolution of a 2D copolymer-solvent system was simulated. The system was simulated with a 512x512 grid with ∆x = 0.5, ∆t = 1.5 and periodic boundary conditions for both order parameters. A pseudo-spectral method was used for spatial derivatives [27]. The initial conditions correspond to random fluctuations in both order parameters. For the copolymer, c 2 = 0.1, u 2 = 1, D 2 = 1, α = 1, φ = 0 was used. For   this set of parameters, the structure corresponds to a lamellae structure. For the solvent, c 1 = 0.3, u 1 = 1, D 1 = 0.5 was used. The interaction parameters were b 1 = 0.07, b 2 = 0.3, b 3 = 0, b 4 = 0.1. Fig. 4 illustrates the time evolution of the copolymersolvent system phase separation. The b 1 parameter determines the preference or affinity of the solvent for one of the copolymer blocks. For the value b 1 = 0.07, there is a strong affinity with one of the blocks, which results in the formation of copolymer micelles. It should be noted that in most parts of macrodomains, the lamellar domains tend to be parallel to the macrophase interfaces and the domains are surrounded by a thin layer of one block. The decrease in the b 1 value to a b 1 = 0.03 value is illustrated in Fig. 5. A lower b 1 value results in similar affinity between the copolymer blocks and the solvent. As seen in Fig. 5, the interface of the copolymer domains is uniformly composed by both blocks. For this value the stripe domains emerge perpendicularly to the macrophase interfaces. Thus a clear morphological transition is brought about by changing the interaction parameter b 1 . Parameter b 1 allows the interaction between the solvent and the copolymer blocks to be modeled. The morphological change occurs at about b 1 = 0.04. This behavior coincides qualitatively Nano-droplet of copolymer in lamellae phase obtained by the dewetting process. The image shows the stepwise shape of the droplet induced by competition between the surface tension on the liquid surface and the lamellae structure, forming the copolymer droplet. The image was extracted from the work of Croll et al. [33].
with the experimental results [32]. To illustrate this point we simulate the formation of a nano-drop on a rigid substrate. To simulate the temporal evolution of a copolymer droplet deposited on a rigid substrate, the air copolymer system was modeled as a ternary system (solvent copolymer), where the air was treated as a bad solvent. The system was simulated with a 512x512x256 (x, y, z) grid with ∆x = ∆y = ∆z = 0.3 and ∆t = 1.5. Periodic boundary conditions in the x and y directions and no flux in the z direction were used. The Eyre algorithm developed in the previous section was used for time evolution, and spatial derivatives were resolved using a pseudo-spectral method. Starting from a disordered state, the time evolution is illustrated in Fig. 6. The final equilibrium structure is illustrated in Fig. 7. Experimental work in the literature shows a similar formation of structures within the nano-droplets. Fig. 7, insert b, shows the image of a nano-droplet obtained by the dewetting process of a diblock polystyrene-polymethyl methacrylate (PS-PMMA) copolymer in the lamellae phase. Details can be found in the work of Croll et al. [33]. The images show good qualitative agreement with the nano-drop obtained by numerical simulation.

X. Conclusion
In summary, the numerical scheme developed by Eyre for unconditionally stable gradient systems was developed and implemented in the simulation of block copolymer systems. Numerical evaluations allow estimation of a higher resolution speed of up to 100 times the traditional method used, called CDS. Subsequently, the development for the simulation of a solvent-copolymer system was extended. In this case, examples of application and potential use were presented in the dynamics resolution of micelle formation and free interfaces, among others. These will be evaluated in later work.