Density and concentration field description of nonperiodic structures
DDensity and concentration field description of nonperiodic structures
Andreas M. Menzel ∗ Department of Physics, University of Illinois at Urbana-Champaign,Loomis Laboratory of Physics, 1110 West Green Street, Urbana, IL 61801, USA (Dated: April 15, 2019)We propose a simple nonlocal energy functional that is suitable for the continuum field charac-terization of nonperiodic and localized textures. The phenomenological functional is based on thepairwise direction-dependent interaction of field gradients that are separated by a fixed distance. Inan appendix, we describe the numerical minimization of our functional. On that basis, we investigatethe kinetic evolution of thread-like stripe patterns that are created by the functional when we startfrom an initially disordered state. At later stages, we find a coarse-graining that shows the samescaling behavior as was obtained for the Cahn-Hilliard equation. In fact, the Cahn-Hilliard modelis contained in our characterization as a limiting case. A slight modification of our model omits thiscoarse-graining and leads to nonperiodic stripe phases. For the latter case, we investigate the tem-poral evolution of the defects (end points) of the thread-like stripes. In view of possible applicationsof this functional, we discuss a possible characterization of polymeric systems and vesicles. Thestatistics of the growth of the thread-like structures is compared to the case of step-growth poly-merization reactions. Furthermore, we demonstrate that the functional may be applied for the studyof vesicles in a continuum field description. Basic features, such as the tendency of tank-treading insimple shear flows and parachute folding in pipe flows, are reproduced.
PACS numbers: 81.16.Rf,83.80.Qr,87.16.D-,81.16.Dn
I. INTRODUCTION
For several decades, concentration and density field de-scriptions have been used to model the behavior of pe-riodic micro- and mesoscopic systems. A prominent ex-ample is given by the Ohta-Kawasaki characterization ofmesophases as they are observed in block copolymers [1].This description was derived by averaging over the localpositions of the monomers that form the different blocksof the polymers. The procedure led to an effective contin-uum characterization of the observed periodic structuresby the famous Ohta-Kawasaki energy functional. In thecase of incompressible diblock copolymers, the order pa-rameter is given by the concentration of one of the twomonomer types that form the different blocks. This iswhy we call this sort of characterization a concentrationfield description in the following.Since then, a lot of studies have demonstrated the ef-fectiveness of the approach. For instance, the kineticevolution of the periodic mesophase textures from thedisordered (spatially homogeneous) state has been inves-tigated [2]. Also the kinetic evolution from one textureto the other after temperature quenches has been ana-lyzed [3, 4]. Elastic moduli characterizing deformationsof the textures were determined [5, 6], and rheologicalproperties of such structures were studied [7–11].In a recent model of similar spirit, an energy functionalhas been postulated to characterize the behavior of crys-talline structures through a density field description. We ∗ Current address: Max Planck Institute for Polymer Research,P.O. Box 3148, 55021 Mainz, Germany; email: [email protected] refer to the phase field crystal model suggested by El-der and Grant [12, 13]. Periodicity of the density fieldis put into the model by energetically stabilizing at leastone nonzero wave vector in the ground state. The mag-nitude of the subsequent periodic density field is inter-preted as the averaged probability to find a constituentof the crystal. Short timescales of motion are regardedas averaged over. Because of this, the approach is veryeffective in modeling processes in crystals on timescalesthat are long from a molecular dynamics point of view.We refer to such a model as a density field descriptionbecause the order parameter is related to a mass densityrather than a concentration field.Kinetic evolution of periodic patterns from a disor-dered state and corresponding phase diagrams have beenanalyzed in the framework of the phase field crystalmodel [13]. Furthermore, the approach has been usedto describe the motion of defects in crystalline structures[14, 15], including avalanche events [16]. The phase fieldcrystal model has been shown to be a starting point formultiscale simulations of multidomain structures [17, 18].In addition, the suitability for the characterization of epi-taxial growth and crack propagation has been pointedout [12, 13]. It has been demonstrated that averagingthe trajectories obtained from molecular dynamics simu-lations indeed leads to density fields that can be capturedby the phase field crystal model [19]. The connection toclassical density functional theory is of current interest[20, 21].Both of the two approaches outlined above can beinterpreted as describing a sort of phase separation onmicro- or mesoscopic length scales. This means thatmacroscopic coarse-graining of the observed structures,which sets in, for example, in the case of the Cahn-Hilliard model [22], is suppressed. The structures that a r X i v : . [ c ond - m a t . s o f t ] A ug are obtained are periodic in space. Or, at least, if wethink of defect-rich and/or multi-domain structures, theenergetically favored ground state is periodic in space. Itis our goal in this paper to describe a sort of microphaseseparation that does not feature this tendency and isnonperiodic in space. Whereas the phase field crystalmodel has been successful to characterize processes incrystalline materials, our approach may be interesting tomodel amorphous materials such as polymers.In the next section we phenomenologically introducesuch a nonlocal energy functional and motivate the con-cept behind it. We explain how we numerically solve thecorresponding kinetic equation in the appendix. Afterthat, in section III, we give an overview on the texturesthat are obtained from minimizing the energy functional.In particular, we focus on the time evolution of the struc-tures as they develop from the disordered state. At thisstage, coarse-graining to macrophase separation still setsin. We therefore present a modification to our model insection IV, which leads to stabilized stripe phases that donot phase separate macroscopically. Some possible appli-cations, predominantly in the field of soft matter systems,are pointed out in section V, before we conclude. II. SIMPLE NONLOCAL ENERGYFUNCTIONAL
As pointed out above, we now introduce a simpleenergy functional that produces localized nonperiodicstructures. The functional is introduced in a phenomeno-logical way. However, we will shortly explain the conceptbehind the chosen functional form using the illustrationin Fig. 1. After that, we will point out the sort of tex-tures obtained. We will restrict ourselves to two spatialdimensions in this paper.Our goal is to use a single scalar order parameter field ψ ( r ) to describe localized nonperiodic objects, as for ex-ample the ones shown in Fig. 1(a). ψ ( r ) denotes a con-centration or density field, similar to the continuum fielddescription on block copolymers [1] and the phase fieldcrystal model [13]. The spheres and stripes shown inFig. 1(a) feature a characteristic size d that we externallyimpose. It corresponds to the diameter and thickness ofthe spheres and stripes, respectively. In this example,the distances between the objects do not follow a peri-odic rule. More details are given below.Without loss of generality, we identify the black andwhite regions in Fig. 1(a) with values ψ = − ψ = +1, respectively. Fig. 1(c) shows the density orconcentration profile along a horizontal cut through oneof the stripes.We now consider a single isolated up-gradient, asshown, for example, in the left half of Fig. 1(c). It wouldbe energetically unstable in periodic approaches like theones introduced above [1, 13] because it is not embeddedin an overall periodic structure. Our idea is to balanceeach such local up-gradient by a down-gradient that is i j (a) i -1012 ψ i , j = , ψ i , ¯ j (b) profile at j=15: ψ i,j =15 vertically averaged profile: ψ i, ¯ j
70 75 80 85 90 95 100 i -1012 ψ i , j = , ψ i , ¯ j (c) ∆ =14 profile at j=15: ψ i,j =15 vertically averaged profile: ψ i, ¯ j FIG. 1: (Color online). (a) Textures observed from a numer-ical iteration of Eqs. (1) and (2) forward in time, startingfrom a homogeneous disordered state. Brightness indicatesthe value of ψ , where black and white correspond to valuesaround − ×
20. Only a domain of size 150 × ψ i,j =15 along a horizontal line at height j = 15 inFig. (a) (+ symbols) and ψ i, ¯ j averaged over the j -directionof Fig. (a) ( × symbols). (c) Zoom into the area in (b) sur-rounded by the dashed rectangle. Gradients ∇ ψ at the latticesites pairwisely sum up to zero, as indicated for one pair. Thedistance between the points of each such pair is ∆ = 14 latticedistances. (Technical details: ϑ = 5, α = 3, d = 15, L = 1,¯ ψ = − .
4, lattice distance dx = 1, 800000 time steps of stepsize dt = 0 . t = 640). located a distance d away in the direction the gradientis pointing to. This direction is given by the unit vector ∇ ψ/ (cid:107)∇ ψ (cid:107) . In this way, we form localized pairs of gra-dients that energetically balance each other, having thesame magnitude but opposite sign and being separatedby an imposed distance d ( ≈ ∆ in Fig. 1(c), see below).Translating this into a formula, the energy functional F that we will investigate is based on the energy density F ( r ) = ϑ (cid:0) ψ | r − (cid:1) + α (cid:104) ∇ ψ | r + ∇ ψ | r + d ∇ ψ (cid:107)∇ ψ (cid:107) (cid:105) . (1)Here, ϑ is a measure for the inverse temperature. Thepart in square brackets is weighted by a second constantparameter α and contains the impact of the gradient field ∇ ψ . Here, the first term refers to ∇ ψ at position r andis local. The second term, however, refers to the valueof ∇ ψ at a position of distance d away from r in the di-rection that ∇ ψ points to. This latter term is thereforea nonlocal contribution. Expression (1) satisfies the nec-essary symmetry requirements. To guarantee thermody-namic stability, we must set ϑ > α >
0. We obtainthe energy functional F{ ψ ( r ) } via F = (cid:82) F ( r ) d r .The kinetic equation for the time evolution of the con-served field ψ ( r ) follows from minimizing the energy F , ∂ψ∂t = L ∆ δ F δψ . (2)Mostly in this paper, the parameter L will be set equalto one, L ≡
1, which can always be achieved by rescalingthe time. On the right-hand side of this equation, theLaplacian reflects the fact that ψ is a conserved quan-tity. The nonlocal contribution of Eq. (1) is implicitlycontained in the functional derivative δ F /δψ .We must remark at this point that, through the form ofEq. (2), we limit our investigations to the regime of relax-ation dynamics. Processes on time scales that are shorterthan the diffusive one have to be considered as being av-eraged over. Such processes cannot be resolved explicitlyby this approach. As a benefit from this coarse-grainedprocedure, however, longer time scales are accessible bynumerical calculations. The same approach was used inthe introduced studies on block copolymer systems [2–4] and the phase field crystal model [12–14, 17, 19–21],or, in a similar spirit, also in classical dynamical densityfunctional theory [21, 23–25].Due to the nonlocal contribution, we could not per-form the functional derivative on the right-hand side ofEq. (2) analytically. We present our way of numericalevaluation of the functional derivative in the appendixand the supplemental material [26]. The central point inour discretized numerical evaluation is to track the con-nections between the points r and r + d ∇ ψ (cid:107)∇ ψ (cid:107) . Due to thediscrete nature of the numerical results, we will replacethe continuous position vector r by indices i and j .In the current form, Eq. (1) contains three parameters ϑ , α , and the distance d . It follows directly that α canbe scaled out. For practical reasons, we set α ≡ d →
0, our expression for F leads to the famousCahn-Hilliard equation [22].A first illustrative example created from our nonlocalenergy functional is actually given by Fig. 1. It showsthe result of an iteration of Eqs. (1) and (2) forward intime on a rectangular domain of large aspect ratio andperiodic boundary conditions. The axis labels i and j index the grid sites of the square calculation grid. Weinitialized the field ψ as ψ ij = ¯ ψ + η ij , with ¯ ψ = − . η ij a random field of zero mean, Gaussian distribu-tion, and small amplitude. Sphere and stripe domainsare obtained. Their diameter follows from the value ofthe parameter d , which we set to d = 15. We display in Fig. 1(b) the density profile ψ i,j =15 ona horizontal line at j = 15, and a vertically averagedprofile ψ i, ¯ j = (cid:80) j =0 ψ ij /
21. The values of the two pro-file functions are equal for the stripes and differ for thespheres.Fig. 1(c) is a zoom into the area of Fig. 1(b) that ismarked by the dashed rectangle and depicts the conceptbehind Eq. (1). The first part with the coefficient ϑ en-forces concentrations of ψ = ± ϑ >
0. All other val-ues increase this contribution to the energy of the system.The second part with the coefficient α leads to an ener-getic penalty for gradients ∇ ψ , unless there is a gradientof opposite sign, but same magnitude, located at a dis-tance d away in the direction that ∇ ψ points to. Indeed,this pairwise annihilation of gradients occurs as indicatedfor one pair in Fig. 1(c). (The distance is ∆ = 14 andnot d = 15 because of an integer cast of the term d ∇ ψ (cid:107)∇ ψ (cid:107) in our numerical calculation to index the discrete latticesites). In this way, no energetic penalty arises from thegradient terms in Eq. (1).In summary, we do not confine ourselves to local ex-pressions in the gradient fields. Such local expressionshave the tendency to enforce periodicity throughout thewhole space. This becomes clear from the Fourier trans-forms of the corresponding energy densities. Usually, acertain band of wave numbers is preferred homogeneouslyover the whole space. Our goal is to generate objects thatare localized in space and not embedded in an overalltendency of periodicity. This is different from multigrainstructures or defect-rich textures since in the latter casesa periodic order is still energetically favored. In otherwords, we wish to enforce phase separation on a meso-scopic length scale that leads to amorphous textures asobserved in non-crystallized polymeric systems, for ex-ample, or localized objects such as vesicles floating in afluid environment, using one scalar order parameter field. III. TIME EVOLUTION OF THE TEXTURES
In the following, we describe the kinetic evolution ofthe patterns that are generated by Eqs. (1) and (2). Sincewe compare to results obtained for the Cahn-Hilliardequation, we refer to the order parameter ψ as a concen-tration field. We will start our numerical considerationsin all cases from a spatially homogeneous concentrationthat is modulated only by a random field of zero mean,Gaussian distribution, and small amplitude. The numer-ical calculations were performed in the way pointed outin the appendix and the supplemental material [26].We show in Fig. 2 a time series of phase separationfor a mean concentration ¯ ψ = 0. Furthermore, we chose ϑ = 5, α = 3, d = 5, and a grid size of 512 ×
512 lat-tice sites. The brighter the color in the pictures, thehigher the concentration (for each picture the concentra-tion values were normalized by their maximum value).Fig. 3 represents the analogous time series for a meanconcentration ¯ ψ = − .
20 40 60 80120 200 400 8001600 3200 4800 8000
FIG. 2: Time series of the kinetic evolution described by Eqs. (1) and (2), with a mean concentration ¯ ψ = 0. The gray scalerepresents the concentration values, where white and black correspond to the maximum and minimum concentration at eachtime, respectively. Parameter values were set to ϑ = 5, α = 3, d = 5, the grid size to 512 ×
512 sites of distance dx = 1. Timesteps were of size dt = 0 . First, thread-like stripe structures form. Shorterthreads connect to longer threads. The reason for thesesteps of forming longer objects is the higher energy den-sity at the ends of the threads: there is no concentrationgradient within the threads that could balance the lon-gitudinal concentration gradient at the end points. Thisincreases the local energy density (compare Eq. (1) andFig. 1). The upper part of Fig. 4 highlights the spatialdistribution of the energy density at an early stage of pat-tern evolution for the case ¯ ψ = 0. Ends of the threadsclearly stand out. We further investigate this initial pro-cess of forming and connecting to longer threads in thenext section.After that, another process of coarse-graining sets in,which leads to macroscopic phase separation. Blobs startto form at defects in the thread structure and then pref-erentially grow and coarse-grain along the thread con-nections. The smaller the mean concentration ¯ ψ , the lessconnected are the blobs by threads during the processof coarse-graining. As we will show below, this coarse- graining process is similar to the one observed for theCahn-Hilliard model.Energetically, we can explain this coarse-graining to amacroscopic phase separation in the following way. Thestripe patterns as they evolve from Eqs. (1) and (2) corre-spond to a metastable state. This is true even if the high-energetic end points of the threads are excluded from ourconsiderations. Due to pairwise annihilation of the gradi-ents at the boundaries of the stripes, the term with coeffi-cient α can be minimized in Eq. (1). However, the valuesat the boundaries of the stripes deviate from ψ = ± ϑ -contribution.For the blobs, the situation is just the other wayaround. The gradients at the surfaces cannot be anni-hilated, so the α -term in Eq. (1) increases in magnitude.This leads to the high energy density at the surfaces de-picted in the lower part of Fig. 4. Yet, this scenario isenergetically beneficial. In most of the interior of theblobs we find ψ ≈
1, which implies F ≈
0. Altogether,this decreases the total energy of the system due to the
20 40 60 80120 200 400 8001600 3200 4800 8000
FIG. 3: Time series of the kinetic evolution described by Eqs. (1) and (2), with a mean concentration ¯ ψ = − .
4. The othersettings were the same as described in the caption of Fig. 2. high number ratio of interior to surface points.To characterize the coarse-graining process more quan-titatively, we determined the structure function S ( k , t ) = 1 N (cid:42)(cid:88) r (cid:88) r (cid:48) e − i k · r (cid:2) ψ ( r + r (cid:48) , t ) ψ ( r (cid:48) , t ) − ¯ ψ (cid:3)(cid:43) . (3)Here, N is the total number of lattice sites. The sums runover all lattice sites ( r and r (cid:48) correspond to index pairs( i, j ) and ( l, m ), respectively, that index the lattice sites).Likewise, the vector k refers to discrete lattice sites of thecorresponding discrete reciprocal lattice. (cid:104) . . . (cid:105) denotesan ensemble average. In our case, we average over threeindependent numerical runs that only differ in the initialconditions. More precisely, the three numerical calcula-tions start from different spatial realizations of the nar-row Gaussian distribution of concentrations around thespatially homogeneous state.From S ( k , t ), we calculate the circularly averaged structure function S ( k, t ) = 1 k (cid:48) < k + dk (cid:88) k (cid:48) ≥ k S ( k (cid:48) , t ) . (4) k (cid:48) that satisfythe requirement k ≤ k (cid:48) < k + dk , where we set dk = 2.Despite the anisotropy in our numerical results, the circu-larly averaged S ( k, t ) is suitable for what we demonstratebelow.In Fig. 5, we plot the structure function S ( k, t ) thatcorresponds to the case of ¯ ψ = 0 and d = 5 (compare thetime series in Fig. 2). Fig. 5(a) shows S ( k, t ) at the earlytimes, when the stripe texture leads to a peak around k ≈
65 (we do not display the higher order peaks). Thispeak indicates a periodic length of about 512 / ≈ dx of the thread-like stripe structures.Comparing with the first row of pictures in Fig. 2, weconfirm this length scale: the threads have a thicknessof d = 5, and, where they are compactly packed, areseparated by roughly a distance d/
2, which gives d + d/ ≈
8. As time proceeds, this peak becomes smaller
FIG. 4: (Color online). Energy density F calculated accordingto Eq. (1) for two of the snapshots presented in Fig. 2: theones taken at time t = 20 and at time t = 8000. The reddishcolor marks regions of high energy density. Clearly, the areasaround defects and at the surface of the blobs stand out. Thepictures show averages over 100 iterations after the given timeto smooth the fluctuations that arise from the discrete natureof the numerical implementation. Only a domain of 170 × until it has vanished at t ≈ k ≈
65 decreases, another peak atsmaller k -values evolves. This peak corresponds to theonset of the second stage of coarse-graining. Fig. 5(b) il-lustrates the development of S ( k, t ) at times longer than t = 800. The peak shifts to smaller k -values, indicat-ing a growth of the characteristic size of the observedstructures. In particular, this is reflected by the last rowof pictures in Fig. 2 and indicates a macroscopic phaseseparation.To discuss the latter statement in more detail, wedemonstrate that the same scaling behavior is found asfor the later stages of coarse-graining in the Cahn-Hilliardmodel [22, 27]. We follow Ref. [27] and define a normal- k S ( k , t ) (a) k S ( k , t ) (b) k/k ( t ) k ( t ) s ( k , t ) (c) FIG. 5: (Color online). Circularly averaged structure function S ( k, t ) for a texture evolved from Eqs. (1) and (2) for ¯ ψ = 0.The curves correspond to times t as indicated on the upperright of each plot. See Fig. 2 for a real space representationand the parameter values. (a) Temporal evolution of S ( k, t )during the early stage. The decreasing first order peak at k ≈
65 represents the thread-like stripe pattern, the increasingpeak at lower wave numbers develops due to the macroscopiccoarse-graining. Only data S ( k, t ) <
100 are shown. (b)Evolution of S ( k, t ) at the later stage of coarse-graining. Asthe blobs grow in size, the peak shifts to lower k -values. Thedashed lines are guides to the eye only. (c) Data collapse forthe later stage of coarse-graining (see text and Ref. [27]).
40 240640 1600
FIG. 6: Time series of the kinetic evolution described byEqs. (1) and (2), with a mean concentration ¯ ψ = − .
7. Theother settings were again the same as described in the captionof Fig. 2. However, a nucleation core of 5 × ψ = 0 . ized structure function s ( k, t ) = S ( k, t ) N [ (cid:104) ψ (cid:105) − (cid:104) ψ (cid:105) ] , (5)as well as a characteristic wave number k ( t ) = (cid:80) k k s ( k, t ) (cid:80) k s ( k, t ) . (6)When we plot k ( t ) s ( k, t ) as a function of k/k ( t ), thedata for the later times collapse as displayed by Fig. 5(c).After the analysis of this scaling behavior, we furtherincreased the magnitude of the mean concentration | ¯ ψ | .For values | ¯ ψ | ≥ .
6, we found that the patterns didnot form spontaneously. This can also be understoodfrom a comparison with the Cahn-Hilliard equation. Inthe Cahn-Hilliard model, an initially homogeneous stateonly becomes linearly unstable for mean concentrations | ¯ ψ | ≤ / √ ≈ . | ¯ ψ | ≥ .
6. Fig. 6shows an example, where we set ¯ ψ = − .
7. The otherparameters were chosen the same as in Figs. 2 and 3.However, at t = 0, a spot of 5 × ψ = 0 .
5. Indeed, we thenobserve a process of growth after providing the initiat-ing nucleation core. Coarse-graining starts at the mostunstable points in the threads of Fig. 6.Finally, we display in more detail an example of howdefects can form, which then initiates the coarse-graining. (a) (b) (c)(d) (e) (f)
FIG. 7: Snapshots of 57 ×
57 lattice sites from a numericalintegration of Eqs. (1) and (2), with a mean concentration¯ ψ = − .
4. We set the parameter values to ϑ = 5, α = 3, d = 15, the total grid size was 512 ×
512 sites of distance dx = 1, and time steps were of size dt = 0 . t that the different pictures correspond to are (a) 5760, (b)6400, (c) 6480, (d) 6760, (e) 7200, and (f) 8000. The pictures in Fig. 7 show snapshots from a numeri-cal calculation of a mean concentration ¯ ψ = − . d = 15. The series starts from twothreads (a), the energy density at the ends of which isincreased as noted above. First, the “reactive” end of theshorter thread couples into the side of the other thread(b). The resulting structure (c) is frustrated since the re-quirement of pairwise annihilation of the gradients can-not be satisfied around the connection point. This higherenergy density around the crosslinking point initiates thecoarse-graining (d). A growing blob incorporates the re-maining thread texture (e). Finally, we end up with aspherical object (f) that looks reminiscent of a small vesi-cle.We shortly remark at the end of this section, that thepreferentially diagonal, horizontal, and vertical orienta-tion of the threads is due to our discretization schemeof the gradient field as described in the appendix. It re-flects the underlying square lattice that acts similar toan external orienting field. Varying the value of the pa-rameter d , we checked that the coarse-graining is not dueto the interplay with the discrete numerical grid. Still,the amorphous and localized character of the structuresis evident, in particular for lower mean concentrations ¯ ψ (Figs. 3 and 6). We checked the numerical scheme byvarying the size of the time steps dt , the distance of thelattice sites dx , the characteristic length d (also to non-integer values), as well as the discretization scheme of theLaplacian in Eq. (2). IV. STABILIZED STRIPE PHASES
To make our model more attractive for the descrip-tion of amorphous and localized structures, we needto stabilize the stripe state against macroscopic coarse-graining. This becomes clear when we consider surfac-tant molecules as an example. The white regions inFigs. 2 and 3 are identified with the hydrophobic compo-nents, the black regions with the hydrophilic components.Then macroscopic phase separation is not allowed. Thesurfactant molecules cannot simply chemically split theirhydrophilic from their hydrophobic parts. At the sametime, each gradient in ψ is associated with an accumu-lation of the hydrophilic component on the up-side andthe hydrophobic component on its down-side.In the following minimal approach, we only considerthe accumulation on the down-side of each gradient. Wedo this by adding an additional contribution F β to theenergy density in Eq. (1), F β ( r ) = β (cid:104) ψ | r − b ∇ ψ (cid:107)∇ ψ (cid:107) + 1 (cid:105) . (7)Due to this contribution, each gradient ∇ ψ favors a valueof ψ ≈ − r , ifthe value of ψ deviates from ψ = − b away in the direction − ∇ ψ (cid:107)∇ ψ (cid:107) | r . The interplay with the α -contribution in Eq. (1) energetically disfavors macro-scopic coarse-graining.We can understand the latter statement by looking atthe fine structure of the blob boundaries in Figs. 2, 3, 6,and 7. The α -term in Eq. (1) enforces the formation ofconcentration shells around the blob, which correspondsto the scheme depicted by Fig. 1. It leads to concentra-tion gradients at the inner shell boundary pointing intothe inward direction. However, the inside concentrationof the blob is ψ ≈ +1, which leads to a high energeticpenalty due to the β -contribution in Eq. (7). There-fore, macroscopic coarse-graining is energetically inhib-ited. Details of the numerical implementation of the β -term are given at the end of the appendix and in thesupplemental material [26].Figs. 8 and 9 illustrate the stability of the stripe pat-terns when Eq. (7) is included. The mean concentrationswere given by ¯ ψ = 0 and ¯ ψ = − .
4, respectively. Weset the other parameter values to ϑ = 5, α = 3, d = 5, β = 1, and b = 3. Again, the grid size was 512 ×
20 8000
FIG. 8: Stability of the stripe state during a kinetic evolutiondescribed by Eqs. (1) and (2), supplemented by Eq. (7), for amean concentration ¯ ψ = 0. In contrast to Fig. 2, no macro-scopic coarse-graining occurs. Parameter values were set to ϑ = 5, α = 3, d = 5, β = 1, and b = 3. The grid size was512 ×
512 lattice points of distance dx = 1, and time stepswere of length dt = 0 .
20 8000
FIG. 9: Stability of the stripe state during a kinetic evolutiondescribed by Eqs. (1) and (2), supplemented by Eq. (7), fora mean concentration ¯ ψ = − .
4. No macroscopic coarse-graining occurs, in contrast to Fig. 3. Parameter values asgiven in the caption of Fig. 8. points are generated. The latter hardly occurs at lowerconcentrations (compare Figs. 8 and 9 at time t = 8000).We have to note that the influence of the underlyingsquare grid becomes rather obvious at the later numericaltimes.In order to describe the process of coarse-graining tothe longer threads more quantitatively, we analyzed thekinetics of the end points. We refer to these points asdefects in the following. There are several qualitativelydifferent ways defects can evolve. Defects can appeartogether with the threads when the latter form from thehomogeneous disordered state. This process mainly takesplace at the very early stage of our numerical calculationsand will not be considered. Apart from that, defects canappear when crosslinking points disconnect (inversion ofthe process shown in Fig. 7(a)-(c)). Two defects can uniteto one defect when threads become very short and theend points of the same thread meet. We have checkedthat this process does not play a major role in whatfollows. Defects further vanish by forming crosslinkingpoints (in analogy to Fig. 7(a)-(c)), or by connection of t l n ( N d e f ) numerical dataline fit to intermediate regime t › ∆ r fi FIG. 10: (Color online). Number of defects (end points of thethreads) N def in the stripe patterns as a function of numericalcalculation time t in a semilogarithmic plot. As marked bythe straight line, the defect number decays exponentially ina broad regime from t = 800 to t = 6400. The inset showsthe absolute displacement distance of the defects ∆ r per timeinterval ∆ t = 20 and averaged over the whole numerical sam-ple. Dashed vertical lines mark the interval used for the linefit. Parameter values as given in the caption of Fig. 8 two threads to one at the end points. The latter processdominates particularly at lower mean concentrations ¯ ψ .Fig. 10 shows the number of the defects N def on acalculation grid of 512 ×
512 lattice sites as a functionof time. The mean concentration was ¯ ψ = 0, and theother parameter values were chosen as mentioned in thecaption of Fig. 8. Directly after the initial quench fromthe homogeneous state to a temperature given by ϑ = 5,the number of defects decreases rapidly. As given by thestraight line in the semilogarithmic plot, we then founda broad regime from about t = 800 to t = 6400, wherethe defect number decreases exponentially. At later timesthe annihilation of defects seems to saturate. Apart fromthat, we tracked the motion of the defects. In the regimeof exponential decay, the average distance of motion ofthe defects (cid:104) ∆ r (cid:105) was relatively constant with about 0 . t = 20. The correspond-ing standard deviation was less than 0 .
04 lattice sites.We have plotted the time dependence of (cid:104) ∆ r (cid:105) in the in-set of Fig. 10. There, the dashed vertical lines mark theinterval that we used for the line fit in the main part ofthe plot.We performed three independent numerical runs withdifferent initial spatial Gaussian fluctuations around thehomogeneous state. Out of these three runs, we chosethe one where the defect number saturated the latest.The other two samples with earlier defect number sat-uration did not reveal a well-defined intermediate expo-nential regime.At the end of this section we point out that variationsof the parameters ϑ , d , β , and b allow for further, quali-
20 60 1320 8000
FIG. 11: Time series of the kinetic evolution of Eqs. (1),(2), and (7), with a mean concentration ¯ ψ = 0. Parametervalues were set to ϑ = 8, d = 15, β = 0 .
8, and b = 10,the grid size was 200 ×
200 lattice sites of distance dx = 1,and time steps were of size dt = 0 . d . Then partsof the threads grow thicker than d , and finally a bubble-likestructure of homogeneous bubble size and wall thickness d emerges. tatively different structures. We give one example in thetimeseries of Fig. 11. There, the mean concentration wasset to ¯ ψ = 0, and the other parameters were chosen as ϑ = 8, d = 15, β = 0 .
8, and b = 10. For higher thick-nesses d , we usually observe that an interlaced structureof thinner stripes forms from the spatially homogeneousstate (first picture of Fig. 11). These stripes then becomeunstable with respect to the formation of the threads ofthickness d . In the case of the current set of parameters,the threads then partially thicken to threads of a crosssection higher than d (third picture of Fig. 11). Finally,the latter transform into a foam-like texture with wallsof thickness d and uniform bubble size. The macroscopicphase separation and coarse-graining into blobs as ob-served in Figs. 2, 3, and 6 does not occur for this set ofparameters. V. POSSIBLE APPLICATIONS
In this section we discuss two possible applications ofour energy density (1) and (7) to model systems in thefield of soft condensed matter. In the first case, ψ ( r ) isinterpreted as a density field, which is therefore of thesame spirit as is the phase field crystal model [13] forperiodic systems. In the second part, ψ ( r ) should ratherbe interpreted as a concentration field. This can be seenin analogy to the concentration field description of blockcopolymer systems [1] in the periodic case.From the density field point of view, as already men-tioned in the introduction, it appears natural to useEqs. (1) and (7) for a density field characterization ofpolymeric systems. ψ ( r ) is then connected to the timeaveraged probability of finding a segment of a polymerchain at position r , in analogy to the phase field crystalapproach [19]. Our model then describes polymeric sys-tems on the length scale of a segment and on diffusivetime scales.Two steps on this way to a realistic microscopic con-tinuum characterization of polymer systems seem to beimportant. First, we need to work out a more isotropic0discretization scheme than the one presented in the ap-pendix. This is necessary to decouple the orientation ofthe threads from the orientation of the calculation grid.Second, several topological properties of polymer systemscan only arise in three dimensional space. This concernstopological entanglements that are not possible in twospatial dimensions. Special, yet illustrative, examples in-clude mixtures of ring and linear polymers, where therings can be threaded by the linear chains [28, 29], andslide-ring gels [30–32]. It will therefore be interesting togeneralize our approach to three dimensional space.These points are important, e.g., for rheological consid-erations. Here, we leave them for future work. Instead,we now concentrate on the kinetics of the reaction pro-cess illustrated for example in Fig. 9. We note that theprocess of phase separation and formation of the stabi-lized stripes reflects several characteristics of step-growthpolymerization reactions [33, 34]. In these reactions, allreactants/monomers can participate in the process fromthe beginning and are quickly incorporated into shortchains. The chains can grow at both ends, and shortchains can connect to longer ones. This means that thechain length continuously grows with increasing reactiontime. Long reaction times are necessary to obtain highdegrees of polymerization (long chains). All of these fea-tures are reflected by the process referred to in Fig. 9.As common in the literature, we now denote by x thenumber of units/monomers that a chain consists of. N gives the total number of units/monomers in the sample,whereas N denotes the total number of molecules of allsizes at a specific reaction time. Then the expression p = ( N − N ) /N describes the extent of reaction at thatspecific reaction time. As shown by Flory, the number N x of molecules that consist of x units/monomers is givenby [34–36] N x = N (1 − p ) p x − . (8)We now compare this distribution to the results fromour numerical solution of Eqs. (1), (2), and (7) for ¯ ψ = − . ψ ≥ N x , weextracted at each time the number of threads that consistof x connected pixels of ψ ≥ N is given by the totalnumber of threads (connected objects of ψ ≥
0) at thistime. In order to identify the total number of all units N , we counted the total number of pixels of ψ ≥ t = 800.The distributions were determined from nine separatenumerical samples of different initial fluctuations aroundthe spatially homogeneous state. Each sample was of size512 ×
512 grid points, the parameter values were chosenas given by the caption of Fig. 9. For small values of x , weexpect a significant deviation of N x from our numericalresults. The reason is that our process of phase separa-tion does not support objects of a diameter smaller than d . Also, diffusion can be important for small objects, be-fore the process of chain connections becomes dominant. x N x FIG. 12: (Color online). Time evolution of the chain/threadlength distribution N x , where x gives the number of units inthe chain/thread. The numerical data points were obtainedfrom nine independent numerical samples at the times indi-cated on the upper right. Corresponding parameter values aregiven by the caption of Fig. 9. Each point is an average over20 sizes. The lines follow from the Flory chain length distri-bution for step-growth polymerizations as given by Eq. (8).The necessary parameters were extracted from the numericalresults. Fig. 12 compares the results from our numerical solu-tions to the predictions of Eq. (8). The numerical datapoints are averaged over 20 object sizes. We chose thissize of averaging because it corresponds to the minimalsupported object size of π ( d/ . Clearly, we can see thedeviations outlined above. Considering, however, the factthat there is no adjustable parameter, the description byEq. (8) works reasonably well. It allows an approximateconnection to step-growth polymerization processes.In the second part of this section, we explain that aconcentration field interpretation of Eqs. (1), (2), and(7) can be seen as a phenomenological approach to qual-itatively describe micellar, vesicular, and bilayer mem-brane systems. For this purpose, we think, e.g., of sur-factant molecules in an aqueous solution. The surfactantmolecules consist of a hydrophilic head group and a hy-drophobic tail. We assume that the mass density of thehydrophilic and hydrophobic parts is identical. Then wecan characterize the pure hydrophilic components (waterand hydrophilic head groups) by a concentration ψ = − ψ = +1.The phase separation for high enough values of ϑ cor-responds to the separation of hydrophilic and hydropho-bic components. For example, in Fig. 1(a), the sphereswould correspond to micelles, and the stripes to bilayermembranes. The hydrophilic head groups are on averagelocated on the dark side ( ψ ≈ −
1) of the dark-to-brightboundaries. On the contrary, the hydrophobic tails pointinto the center of the bright regions ( ψ ≈ +1).Since hydrophilic and hydrophobic parts are connectedto each other, the width of the hydrophobic region can atmaximum be twice the average length of the hydrophobic1tails. We therefore roughly identify the thickness d ofthe bright regions in Fig. 1(a) as twice the length of thehydrophobic tails. The β -term in Eq. (7) determines thedegree of depletion of hydrophobic groups around theregions of hydrophilic head groups.It may be possible to relate processes similar to theones shown in Fig. 7 to the formation of a vesicle. In theearly stages of the numerical calculations, or throughoutthe calculations of section IV, the threads would haveto be interpreted as wormlike micelles. The higher cur-vature energy density at the ends of wormlike micellescan be mimicked by our energy functional, as, e.g., de-picted in the upper picture of Fig. 4. In the following,however, we will restrict ourselves to qualitative hydro-dynamic considerations. More precisely, we investigatethe flow behavior of vesicular structures in simple shearand pipe flows.For this purpose, we coupled Eqs. (1), (2), and (7) withthe Navier-Stokes equation, which leads to [2] ∂ψ∂t = − v · ∇ ψ + L ∆ δ F δψ , (9) ∂ v ∂t = −∇ p − v · ∇ v + η ∆ v − ( ψ − ¯ ψ ) ∇ δ F δψ . (10)Here, we have assumed incompressibility. This makes thecontinuity equation reduce to ∇ · v = 0. L is the sameparameter as in Eq. (2), and η denotes the viscosity. The(constant) mass density ρ has been scaled out.We set the parameter values to ϑ = 5, α = 3, d =15 . β = 0 . b = 10 .
2, and η = 100. In the case ofsimple shear flow, we set L = 10, whereas we chose L =20 in the case of pipe flow. The numerical calculationgrid was of size 256 ×
128 sites with a lattice distance of dx = 1, and the size of the time steps was dt = 0 . − g = (0 , g ) to Eq. (10), where g < ∇ · v = 0 was solvedthrough a spectral method [38]. Third and first order up-wind schemes were used to discretize the advective termin Eqs. (9) and (10), respectively.In order to start our investigations, we provided anannulus of ψ = +1 and thickness d in an environment of ψ = −
1. We let it relax to a spherical vesicle-like objectby numerically iterating Eq. (2) for a time interval ∆ t = FIG. 13: (Color online). Simple shear flow of a vesicle char-acterized by Eqs. (9) and (10). An approximately sphericalvesicle was provided as an initial condition. At time t = 5 theshear was started by setting the horizontal velocity compo-nent at the upper boundary to a nonzero value, pointing tothe left, the effect of which is clearly visible at time t = 27 . t = 200) and twice (time t = 455). Pink color: ψ > .
8, red color: ψ >
1. Parametervalues were set to ϑ = 5, α = 3, d = 15 . β = 0 . b = 10 . η = 100, and L = 10. The grid size was 256 ×
128 sites with alattice distance of dx = 1, and the size of the time steps was dt = 0 .
40. These vesicles were then used as an initial conditionfor the numerical investigations of Eqs. (9) and (10). Theobjects that we describe by this sort of concentration fieldapproach have to be classified as fluid vesicles. Thatis, viscous flow of the constituents is possible within thevesicle membranes. No elastic shear forces occur withinthe local membrane planes, which would be the case forelastic membranes.We give an example for the obtained behavior of avesicle under shear flow in Fig. 13. The first picture showsthe vesicle at the time when the shear flow is started bysetting the horizontal velocity component to a nonzerovalue at the upper boundary. In the second picture, theeffect of this shear deformation starts to become visible.The bottom row of Fig. 13 depicts the vesicle as it isstrongly deformed by the flow field of steady shear. Weshow its state after having crossed the periodic left-rightboundary once and twice.The long axis of the stretched vesicle approaches asteady orientation with respect to the flow direction.This behavior is well known from numerical studies,e.g. solvent multiparticle collision dynamics combinedwith coarse-grained membrane models [39, 40], as wellas experimental observations [41, 42]. It is referred toas the tank-treading regime, because it is usually accom-panied by a rotational motion of the vesicle membranearound its interior.We must remark in this context that our characteri-zation in its current implementation does not reproducean essential feature of the experimental systems. Realvesicles are practically impermeable for the surroundingsolvent on experimental time scales. Generally, in fielddescriptions, this can be achieved by setting the mem-2brane velocity equal to the local fluid flow velocity [43].In our case, this corresponds to letting L → .
01 with respect tothe vertical component.Apart from the tank-treading motion, at least twoother regimes were identified numerically and experimen-tally [39–42, 44]. If the inclination angle of the long axisof the elongated vesicle is not steady but varies in time,the motion is called tumbling. In addition, a more ir-regular dynamic mode of motion was identified: it in-cludes tumbling, features dynamic shape deformationsaway from the purely elongated state, and was namedswinging mode.Mainly two parameters were found to determine themode of motion. First, this is the reduced volume. Itmeasures the deviation from a spherical shape by relatingthe surface to the volume of the vesicle. The secondparameter is the ratio between the fluid viscosities insideand outside the vesicle. Imposing impermeability of themembrane is crucial to control these parameters. It willtherefore be an important step for future work on ourdescription.We close this section by pointing out our qualitativeresults for the pipe flow geometry. Fig. 15 shows an ini-tially spherical vesicle that has a size of the order of thechannel diameter. When the flow starts, the vesicle is de- x y (a) x y (b)
40 65 90 115 140 165 190 215 x y (c) FIG. 14: (Color online). Flow field v corresponding to thecase of simple shear in Fig. 13 at time t = 455. The red colorindicates the flow field within the membrane ( ψ > . .
01 w.r.t. the vertical one. FIG. 15: (Color online). Time series of a vesicle exposedto pipe flow. First, the vesicle folds into a parachute shape,following the Poiseuille flow profile. Then it tends to elongatealong the flow direction. Parameter values and color codeas given by the caption of Fig. 13, except for L = 20 and g = − . FIG. 16: (Color online). Two vesicles in a pipe geometry. Astime proceeds, the rear vesicle grows on cost of the front vesi-cle. Parameter values and color code as given by the captionof Fig. 13, except for L = 20 and g = − . formed according to the profile of the Poiseuille-like flow.It folds into a parachute-like shape. At longer times, wefind a transition to a shape that is elongated in the flowdirection (compare Ref. [45]).Finally, we put two vesicles into the channel, as de-picted by Fig. 16. Initially, they show the same behavioras the single vesicle in Fig. 15. However, as time pro-ceeds, the front vesicle looses material that is collectedby the rear one. Consequently, the first vesicle shrinks,whereas the second one grows. It would be interesting toknow whether a corresponding process can be observedin an experimental system. VI. CONCLUSIONS
We have introduced a simple nonlocal energy densityin Eq. (1) that was discretized for numerical investiga-tions as explained in the appendix. Starting from thedisordered (spatially homogeneous) state, this potentialdescribes the formation of non-periodic localized thread-like structures. These structures are metastable, and fi- nally we observe a macroscopic phase separation. Thescaling behavior of the latter coarse-graining process issimilar to the case of the Cahn-Hilliard model. We couldomit the macroscopic phase separation by adding expres-sion (7) to the energy density. It stabilizes the localized,non-periodic thread-like stripe structures. In this case,we observe a connecting process of different threads toreduce the number of high-energetic end points of thethreads. Our analysis revealed a regime during whichthe number of end points of the threads decays exponen-tially. Finally, we discussed two possible applications ofour model in the field of soft matter. On the one hand,this was a density field description of polymers, wherefor the moment we concentrated on the chain length dis-tribution function during the formation process. On theother hand, we focussed on the hydrodynamic proper-ties of vesicular structures. Simple shear and pipe flowgeometries were qualitatively addressed, where our de-scription reproduces generic features of the flow behav-ior.As pointed out, several aspects remain for future in-vestigations. We have indicated that the variation of thevarious parameter values involved in the description leadsto further qualitatively different textures. A completephase diagram needs to be established. Next, to allowmore quantitative studies, our numerical discretizationscheme should be improved. For example, a further de-coupling from the underlying grid directions might beachieved through an implementation on a hexagonal lat-tice. Furthermore, a generalization to three spatial di-mensions would be desirable. Lastly, we have followeda phenomenological approach so far. It would be inter-esting to investigate whether our functional can be re-lated to more microscopic approaches through appropri-ate coarse-graining, and which changes to the functionalform this would imply.
Acknowledgments
The author thanks Nigel Goldenfeld and Takao Ohtafor stimulating discussions on related fields, and NigelGoldenfeld for a stay in his group, where this work wasperformed. This work was supported by the DeutscheForschungsgemeinschaft through a research fellowship,which is gratefully acknowledged.
Appendix: Discretized numerical evaluation of thefunctional derivative and kinetic equation
In this appendix we describe our numerical evaluationof the functional derivative in Eq. (2). The form of thecorresponding energy density is given by Eq. (1).Due to the nonlocal contribution in Eq. (1), it is notstraightforward to perform the functional derivative ana-lytically. We use the following definition of the functional4derivative: δ F{ ψ ( r (cid:48) ) } δψ ( r ) = lim (cid:15) → (cid:15) (cid:104) F { ψ ( r (cid:48) ) + (cid:15) δ ( r (cid:48) − r ) }−F { ψ ( r (cid:48) ) } (cid:105) . (A.1)This is the correct starting point for the transformationfrom the analytic continuum picture to the discrete nu-merical calculation.The numerical lattice is of finite size and defines a dis-crete system. We use a rectangular domain on a squarelattice of N x × N y lattice sites with periodic boundaryconditions. The lattice sites are indexed by pairs ( i, j ),where i = 0 , . . . , N x − j = 0 , . . . , N y −
1. Conse-quently, the continuous functional form
F { ψ ( r ) } trans-forms to a simple multivariate function F [ ψ ij ] of N x × N y variables ψ ij .In analogy, the Dirac delta function δ ( r (cid:48) − r ) is nowgiven by a product of Kronecker deltas δ il δ jm , where ( i, j )and ( l, m ) correspond to the positions r and r (cid:48) , respec-tively. Therefore, the functional derivative in Eq. (A.1)transforms to a simple partial derivative, ∂ F [ ψ lm ] ∂ψ ij = lim (cid:15) → (cid:15) (cid:104) F [ ψ lm + (cid:15) δ il δ jm ] − F [ ψ lm ] (cid:105) . (A.2)This is the expression that we use to numerically evaluatethe functional derivatives.Furthermore, we discretize the gradient field in Eq. (1)as ∇ ψ | ij = (cid:18) ψ i +1 ,j − ψ ij dx , ψ i,j +1 − ψ ij dx (cid:19) . (A.3)We can motivate this choice by considering a quadraticgradient contribution ( ∇ ψ | r ) in the energy density, asit occurs, e.g., in the energy functional for the Cahn-Hilliard model [22]. On the one hand, analytically,the functional derivative of this term follows immedi-ately as − ∆ ψ | r . On the other hand, the discretiza-tion of the term ( ∇ ψ | r ) according to Eq. (A.3) reads (cid:2) ( ψ i +1 ,j − ψ ij ) + ( ψ i,j +1 − ψ ij ) (cid:3) /dx . We can per-form the functional derivative as described by Eq. (A.2).Equating both results at position | r ≡ | ij yields∆ ψ | ij = ψ i +1 ,j + ψ i − ,j + ψ i,j +1 + ψ i,j − − ψ ij dx . (A.4)As we can see, we correctly obtain the discretization ofthe Laplace operator as given by the stencil0 1 01 − . (A.5)This discretization of the Laplacian introducesanisotropy [46]. More sophisticated discretizationsthan expression (A.3) must be worked out for thegradient field when a fully isotropic discretization ofthe Laplacian operator should be obtained from the functional derivative. This is beyond the scope of theconsiderations in this paper.In order to write our results in a compact way, wedefine an abbreviation for the discretized nonlocal sub-scripts | r + d ∇ ψ (cid:107)∇ ψ (cid:107) | ij at site ( i, j ): | xp ij ,yp ij ≡ | ˆx · ( r + d ∇ ψ (cid:107)∇ ψ (cid:107) ) | ij , ˆy · ( r + d ∇ ψ (cid:107)∇ ψ (cid:107) ) | ij . (A.6)The pair ( xp ij , yp ij ) indexes the lattice site that the sub-script | r + d ∇ ψ (cid:107)∇ ψ (cid:107) points to (the “ p ” stands for “pointer”).Here, the gradients in the subscript are discretized ac-cording to Eq. (A.3).On the one hand, the value ψ ij contributes to the valueof the energy density at the grid site ( i, j ). This comesfrom the local ϑ -term and from the local term ∇ ψ | r inthe square brackets of Eq. (1). It is crucial to note that,on the other hand, ψ ij contributes to the value of theenergy density at all other grid sites ( l, m ) that point to( i, j ) through the subscript | xp lm ,yp lm ≡ | r (cid:48) + d ∇ ψ (cid:107)∇ ψ (cid:107) . Thelatter part comes from the nonlocal term ∇ ψ | r + d ∇ ψ (cid:107)∇ ψ (cid:107) inthe square brackets of Eq. (1).Consequently, at each time step, we must determinefor each lattice site ( i, j ) all other lattice sites ( l, m ) thatpoint to ( i, j ); i.e. all other lattice sites ( l, m ) satisfyingthe condition ( xp lm , yp lm ) = ( i, j ). In the supplementalmaterial [26] we present a scheme to perform this task atan overall computational cost of linear order, O ( N x N y ).We used dynamic list structures for this purpose.Following the procedure outlined in Eq. (A.2), we ob-tain the contributions to the functional derivative at site( i, j ). From the ϑ -term in Eq. (1), we directly find ϑ (cid:0) ψ ij − ψ ij (cid:1) . (A.7)Next, at site ( i, j ), the local term ∇ ψ | r leads to a con-tribution2 αdx (cid:88) l,m (cid:110) ( ψ l +1 ,m − ψ lm + ψ xp lm +1 ,yp lm − ψ xp lm ,yp lm ) × ( δ l +1 ,i δ mj − δ li δ mj )+ ( ψ l,m +1 − ψ lm + ψ xp lm ,yp lm +1 − ψ xp lm ,yp lm ) × ( δ li δ m +1 ,j − δ li δ mj ) (cid:111) . (A.8)For the nonlocal term ∇ ψ | r + d ∇ ψ (cid:107)∇ ψ (cid:107) , we obtain a con-tribution2 αdx (cid:88) l,m (cid:110) ( ψ l +1 ,m − ψ lm + ψ xp lm +1 ,yp lm − ψ xp lm ,yp lm ) × ( δ xp lm +1 ,i δ yp lm ,j − δ xp lm ,i δ yp lm ,j )+ ( ψ l,m +1 − ψ lm + ψ xp lm ,yp lm +1 − ψ xp lm ,yp lm ) × ( δ xp lm ,i δ yp lm +1 ,j − δ xp lm ,i δ yp lm ,j ) (cid:111) (A.9)5to the functional derivative at site ( i, j ).Finally, we have to add the contribution arising fromthe term with the coefficient β in Eq. (7). It turns outthat the gradients in the subscript | r − b ∇ ψ (cid:107)∇ ψ (cid:107) must be dis-cretized differently than in Eqs. (A.3) and (A.6). Oth-erwise the discrete nature of the numerical calculationleads to a systematic artificial drift.In analogy to Eq. (A.6), we define | xpβ ij ,ypβ ij ≡ | ˆx · ( r − b ∇ ψ (cid:107)∇ ψ (cid:107) ) | ij , ˆy · ( r − b ∇ ψ (cid:107)∇ ψ (cid:107) ) | ij . (A.10)Here, however, the gradients in the subscript are dis-cretized according to ∇ ψ | ij = (cid:18) ψ i +1 ,j − ψ i − ,j dx , ψ i,j +1 − ψ i,j − dx (cid:19) (A.11) instead of Eq. (A.3). For each site ( l, m ) that satisfiesthe condition ( xpβ lm , ypβ lm ) = ( i, j ) we obtain a contri-bution 2 β ( ψ ij + 1) (A.12)to the functional derivative at site ( i, j ).We present a schematic implementation of the aboveresults using a C-like pseudocode in the associated sup-plemental material [26]. [1] T. Ohta and K. Kawasaki, Macromolecules , 2621(1986).[2] M. Bahiana and Y. Oono, Phys. Rev. A , 6763 (1990).[3] M. Nonomura and T. Ohta, J. Phys. Soc. Jpn. , 927(2001).[4] K. Yamada and T. Ohta, J. Phys. Soc. Jpn. , 084801(2007).[5] K. Kawasaki and T. Ohta, Physica A , 223 (1986).[6] K. Yamada and T. Ohta, Europhys. Lett. , 614 (2006).[7] T. Ohta, Y. Enomoto, J. L. Harden, and M. Doi, Macro-molecules , 4928 (1993).[8] M. Doi, J. L. Harden, and T. Ohta, Macromolecules ,4935 (1993).[9] F. Drolet, P. Chen, and J. Vi˜nals, Macromolecules ,8603 (1999).[10] P. Chen and J. Vi˜nals, Macromolecules , 4183 (2002).[11] R. Tamate, K. Yamada, J. Vi˜nals, and T. Ohta, J. Phys.Soc. Jpn. , 034802 (2008).[12] K. R. Elder, M. Katakowski, M. Haataja, and M. Grant,Phys. Rev. Lett. , 245701 (2002).[13] K. R. Elder and M. Grant, Phys. Rev. E , 051605(2004).[14] J. Berry, M. Grant, and K. R. Elder, Phys. Rev. E ,031609 (2006).[15] P. Y. Chan, N. Goldenfeld, and J. Dantzig, Phys. Rev.E , 035701 (2009).[16] P. Y. Chan, G. Tsekenis, J. Dantzig, K. A. Dahmen, andN. Goldenfeld, Phys. Rev. Lett. , 15502 (2010).[17] B. P. Athreya, N. Goldenfeld, J. A. Dantzig, M. Green-wood, and N. Provatas, Phys. Rev. E , 056706 (2007).[18] N. Provatas, J. A. Dantzig, B. Athreya, P. Chan, P. Ste-fanovic, N. Goldenfeld, and K. R. Elder, JOM , 83(2007).[19] P. F. Tupper and M. Grant, Europhys. Lett. , 40007(2008).[20] K. R. Elder, N. Provatas, J. Berry, P. Stefanovic, andM. Grant, Phys. Rev. B , 064107 (2007).[21] S. Van Teeffelen, R. Backofen, A. Voigt, and H. L¨owen,Phys. Rev. E , 051404 (2009).[22] J. W. Cahn and J. E. Hilliard, J. Chem. Phys. , 258(1958).[23] U. M. B. Marconi and P. Tarazona, J. Chem. Phys. , 8032 (1999).[24] A. J. Archer and R. Evans, J. Chem. Phys. , 4246(2004).[25] P. Espa˜nol and H. L¨owen, J. Chem. Phys. , 244101(2009).[26] See Supplemental Material at [ URL will be inserted bypublisher ] for a scheme of numerical implementation.[27] J. Zhu, L.-Q. Chen, J. Shen, and V. Tikare, Phys. Rev.E , 3564 (1999).[28] J. Klein, Macromolecules , 105 (1986).[29] P. J. Mills, J. W. Mayer, E. J. Kramer, G. Hadziioan-nou, P. Lutz, C. Strazielle, P. Rempp, and A. J. Kovacs,Macromolecules , 513 (1987).[30] Y. Okumura and K. Ito, Advanced Materials , 485(2001).[31] T. Karino, Y. Okumura, C. Zhao, T. Kataoka, K. Ito,and M. Shibayama, Macromolecules , 6161 (2005).[32] P. G. de Gennes, Physica A , 231 (1999).[33] P. Munk, Introduction to macromolecular science (Wiley,New York, 1989).[34] J. M. G. Cowie,
Polymers: chemistry and physics of mod-ern materials (Chapman and Hall, New York, 1991), 2nded.[35] P. J. Flory, J. Am. Chem. Soc. , 1877 (1936).[36] P. J. Flory, Principles of polymer chemistry (Cornell Uni-versity Press, Ithaca, 1953).[37] C. A. J. Fletcher,
Computational techniques for fluid dy-namics: Specific techniques for different flow categories ,vol. 2 (Springer, Berlin, 1990).[38] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P.Flannery,
Numerical recipes in C (Cambridge UniversityPress, Cambridge, 1992).[39] H. Noguchi and G. Gompper, Phys. Rev. Lett. ,258102 (2004).[40] S. Meßlinger, B. Schmidt, H. Noguchi, and G. Gompper,Phys. Rev. E , 011901 (2009).[41] M. Abkarian, C. Lartigue, and A. Viallat, Phys. Rev.Lett. , 68103 (2002).[42] V. Kantsler and V. Steinberg, Phys. Rev. Lett. ,258101 (2005).[43] J. Beaucourt, F. Rioual, T. S´eon, T. Biben, and C. Mis-bah, Phys. Rev. E , 011906 (2004). [44] V. Kantsler and V. Steinberg, Phys. Rev. Lett. ,036001 (2006).[45] H. Noguchi and G. Gompper, Proc. Natl. Acad. Sci. USA , 14159 (2005). [46] M. Patra and M. Karttunen, Numer. Methods PartialDifferential Eq.22