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 24. 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 connement is a powerful tool for breaking the symmetry of a structure and obtaining new shapes, diFerent from those obtained in bulk 6, 7. New studies have shown that 2D and 3D con nement induces new morphologies; for instance, copolymers immersed in cylindrical or spherical nanopores. This is produced by the combined e ect of con nement and curvature 8. The new morphologies enhance new applications, such as the development of membranes in fuel cells, membranes for virus ltration, photonic materials and drug dispensers, among others 9-12. Taking into account, rstly, the experimental diculties of analyzing the e ects of con nement and curvature on morphology, and secondly, the success of theoretical models, it is very important to have eficient calculation tools to explore the complex con gurational map of con ned copolymers. Simulation of the formation processes of concient 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 dis- advantages 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 low- temperature state consisting of two phases. The dynamics are controlled by a specic 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 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 di erent 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 corre- sponds to a Euler numerical scheme. The conver- gence of the Euler algorithm strongly depends on the time discretization employed in. Above a certain value of, the system is numerically unstable, presenting "chess-board" instabilities. Consequently, the CDS algorithm is subject to the tmaxx4 23 restriction, where is the spatial discretization. For 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, coars- ening systems present evolution of their characteristic length L(t) and shape L(t) ' tn. 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 ful
l this requirement since is limited.
IV. Eyre algorithm
The algorithm developed by Eyre for gradient systems has the advantage of being unconditionally stable for every 15. The algorithm guarantees stability if the free energy functional F is divided into its contractive FC and expansive FE 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 is ful lled, where Eand are the eigenvalues of the Hessian matrices of the expansive part of the energy and of the total energy, respec- tively. Implementation of the Eyre algorithm for a specic system requires evaluation of the free en- ergy 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:
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.
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:
VI. CDS vs. Eyre comparison
The CDS algorithm registers good convergence in the numerical resolution of the Cahn-Hilliard equation. However, its e ectiveness 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 m(r). Finally, the average domain size L(t) was obtained through the orientational correlation length1. The results are presented in Fig. 2 obtaining a growth law L mt0:330:01 in good agreement with the previous results. In order to compare these methods, the diference 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 df erent 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 Fs and a long-range term the short-range term is the following:
The expression above is related to the parameter that represents the copolymer. The bi 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 di erential equations for conserved order parameters.
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:
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.
this set of parameters, the structure corresponds to a lamellae structure. For the solvent, c1 = 0:3; u1 = 1;D1 = 0:5 was used. The interaction parameters were b1 = 0:07; b2 = 0:3; b3 = 0; b4 = 0:1. Fig. 4 illustrates the time evolution of the copolymersolvent system phase separation. The b1 parameter determines the preference or anity of the solvent for one of the copolymer blocks. For the value b1 = 0:07, there is a strong anity 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 b1 value to a b1 = 0:03 value is illustrated in Fig. 5. A lower b1 value results in similar anity 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 b1. Parameter b1 allows the interaction between the solvent and the copolymer blocks to be modeled. The morphological change occurs at about b1 = 0:04. This behavior coincides qualitatively 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 evo- lution 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 ux 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 nal 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.
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.
Acknowledgements - This work was supported by Universidad Nacional del Sur, and the National Research Council of Argentina (CONICET).