Benchmark Computation of Morphological Complexity in the Functionalized Cahn-Hilliard Gradient Flow
Andrew Christlieb, Keith Promislow, Zengqiang Tan, Sulin Wang, Brian Wetton, Steven Wise
BBenchmark Computation of Morphological Complexity in theFunctionalized Cahn-Hilliard Gradient Flow
Andrew Christlieb
1, 2
Keith Promislow Zengqiang Tan
1, 3
Sulin Wang , ∗ Brian Wetton Steven Wise Abstract
Reductions of the self-consistent mean field theory model of amphiphilic molecule in solventleads to a singular family of Functionalized Cahn-Hilliard (FCH) energies. We modify the energy,removing singularities to stabilize the computation of the gradient flows and develop a series ofbenchmark problems that emulate the “morphological complexity” observed in experiments. Thesebenchmarks investigate the delicate balance between the rate of arrival of amphiphilic materialsonto an interface and least energy mechanism to accommodate the arriving mass. The result is atrichotomy of responses in which two-dimensional interfaces grow by a regularized motion againstcurvature, pearling bifurcations, or curve-splitting directly into networks of interfaces. We evaluatesecond order predictor-corrector time stepping schemes for spectral spatial discretization. Theschemes are based on backward differentiation that are either Fully Implicit, with PreconditionedSteepest Descent (PSD) solves for the nonlinear system, or linearly implicit with standard Implicit-Explicit (IMEX) or Scalar Auxiliary Variable (SAV) approaches to the nonlinearities. All schemesuse fixed local truncation error to generate adaptive time-stepping. Each scheme requires properpreconditioning to achieve robust performance that can enhance efficiency by several orders ofmagnitude. The nonlinear PSD scheme achieves the smallest global discretization error at fixedlocal truncation error, however the IMEX scheme is the most computationally efficient as measuredby the number of Fast Fourier Transform calls required to achieve a desired global error. Theperformance of the SAV scheme performance mirrors IMEX , at roughly half the computationalefficiency.
Keywords:
Phase Field Model; Benchmark Computations; Adaptive Time Stepping; Func-tionalized Cahn-Hilliard. Department of Computational Mathematics, Science and Engineering, Michigan State University, East Lansing,MI 48824, USA Department of Mathematics, Michigan State University, East Lansing, MI 48824, USA School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan 430074, China ∗ Corresponding author Department of Mathematics, University of British Columbia, Vancouver, BC V6T 1Z2, Canada Department of Mathematics, University of Tennessee, Knoxville, TN 37996, USAE-mail addresses: [email protected] (A. Christlieb), [email protected] (K. Promislow), [email protected] (Z.Tan), [email protected] (S. Wang*), [email protected] (B. Wetton), [email protected] (S. Wise)
Submitted to Elsevier in June 2020 a r X i v : . [ phy s i c s . c o m p - ph ] J un Introduction
We present a series of physically motivated computational benchmark problems addressing the evo-lution of a gradient flow that develops rich morphological structures with slightly different energieswhose faithful resolution requires significant computational accuracy. There has been considerablerecent attention to the development of energy stable computational schemes for gradient descentflows [32, 33, 26, 27, 13, 14, 20, 34]. Gradient flows are defined by the dissipation of a free energy,and it is essential that numerical schemes preserve that property. Energy stable schemes have thedesirable property that the energy, or a pseudo-energy, decreases at every time-step irrespective oftime-step size. We argue that energy decay should be a consequence of accuracy, indeed energydecay without accuracy can mask poor performance by leading to plausible but incorrect compu-tational outcomes. Accuracy should be balanced against computational cost, which motivates usto compare computational efficiency between schemes as measured by the minimal computationalcost required to achieve a desired global discretization error.Meaningful assessment of computational efficiency can be achieved from gradient flows that se-lect from significantly different outcomes with subtle energy differences arising from strong nonlinearinteractions. For motivation we emulate the “morphological complexity” experiments conductedby [22]. Dispersing amphiphilic diblock polymers in solvent and then allowing the mixture to relax,they observed the formation of a wide variety of structures whose evolution depended sensitivelyupon the polymer chain properties, see Figure 1 and [2, 1]. Reductions of the self-consistent meanfield theory for a model of amphiphilic molecules in solvent leads to a singular family of Function-alized Cahn-Hilliard (FCH) energies, [6]. We modify these energies, mollifying the singularitiesto produce a family of computationally stable but highly nonlinear FCH gradient flows similar tothose studied earlier, [17, 11, 12]. We present a series of benchmark problems that recover theonset of morphological complexity. These benchmarks reveal a delicate balance between the rate ofarrival of amphiphilic materials onto an interface and the gradient flow’s selection of a least energymechanism to redistribute them after absorption. This rate-based selection mechanism yields atrichotomy of responses in which two-dimensional interfaces either grow by a regularized motionagainst curvature, under-go pearling bifurcations, or directly curve-split into networks of interfaces.We present three numerical schemes, each combining second-order in time backward differentiationwith a predictor-corrector approach and spectral spatial discretization. The FCH energy is com-putationally stiff. Each of the second order methods considered here balances implicit and explicitterms and their efficiency is very sensitive to the choice of the implicit terms, with improvementsof several orders of magnitude possible when well pre-conditioned. These methods include animplicit-explicit (BDF2-IMEX) method and a scalar auxiliary variable approach (BDF2-SAV). Thelatter features provably unconditional pseudo-energy stability properties. Both of these schemes2re linear in their implicit stage. We compare these with a fully implicit, second order, backwarddifferentiation scheme based upon a preconditioned steepest descent with approximate line search(BDF2-PSD) for the nonlinear solve. We will drop the ’BDF2’ term in the sequel for brevity.The FCH gradient flows possess many distinct, emergent timescales that preclude a fixed time-stepping approach. For each scheme a specified target local truncation error is used to generatean adaptive time-stepping procedure. The first set of benchmarks, the sub-critical, critical, andsuper-critical, use a relatively well conditioned form of the FCH energy, and vary the initial data totrigger bifurcations that lead to morphological complexity that engenders many possible outcomesin the super-critical case. A proper resolution of the time evolution requires considerable accuracy.The second set of benchmarks build significant stiffness into the FCH energy by increasing thesecond derivative of the well at the background state, mimicking the singular nature of the FCHenergy. This adds a small foot to the left minima of the well, see Figure 2, hence these benchmarksare called Foot 1 and Foot 2. The stiffness increases the ratio of the absorption to the redistributiontime-scales and also induces morphological complexity without change to initial data.Each of the second order schemes we consider required a judicious choice of implicit terms orpreconditioner. This choice was typically based upon the linearization of a residual about a spatiallyconstant equilibrium solution. The linearly implicit IMEX and SAV accommodated the increase instiffness for the Foot 1 and Foot 2 benchmarks without significant adjustment. The nonlinear solvein the PSD scheme required optimization of internal parameters, in particular an error toleranceassociated to the iterative nonlinear solver, to converge. Moreover the efficiency of the PSD schemedecreased in comparison to the two linear implicit schemes with increasing numerical stiffness.We conduct grid refinement studies to verify that each benchmark has an adequate spatialresolution and develop highly accurate solutions for each benchmark by an extensive computationwith a very small local truncation error. Once spatially resolved, all three schemes yield concordantresults for sufficiently small specified local truncation error. We adjust the local truncation errorrestriction and use short runs to tune performance parameters in each scheme for each benchmark,and record the accuracy and cost of each optimized scheme. At given local truncation error wefind that the PSD approach is generically the most accurate and SAV is the least accurate, asmeasured by global error at the final time. However at fixed local truncation error the IMEX andSAV schemes required less computational effort than the PSD, with the IMEX scheme consistentlyrequiring about half the computational effort of the SAV scheme. For these benchmarks a global L discretization error of 3 × − was found to be a harbinger of accuracy, and within this constraintwe viewed the local truncation error as an internal parameter to be adapted for each scheme tooptimize global performance. For the sub-critical, critical, and super-critical benchmarks, all threeschemes achieved this global accuracy with comparable efficiency although at quite different valuesof the local truncation error. As presented in Figure 15 achieving this accuracy for the super-critical3enchmark requires 1 . × , × , and 3 × FFT calls for IMEX, PSD and SAV respectively. Asthe global error target was further tightened, the PSD scheme encountered increased computationaleffort, first increasing rapidly and then saturating. Conversely the computational effort of the IMEXand SAV schemes increases linearly with global discretization error related by a fixed factor of two.For the more strongly nonlinear benchmarks the efficiency of the linear-implicit schemes continuedits linear relationship to global discretization error, and to each-other by the same factor of two.The efficiency of PSD deteriorated in comparison to the linear-implicit methods with increasingnonlinear strength, as depicted in Figure 16.The SAV scheme is specifically designed to be energy stable with respect to an associatedpseudo-energy. This property either assumes fixed time-stepping, which is impractical for theFCH gradient flows, or an adaptive time stepping based upon modifications by factors of twoand nesting. This latter strategy was implemented for the super-critical benchmark and found toprovide no benefit for accuracy while increasing computational cost by a factor of two to three. Inall cases all convergent schemes preserved the energy decay property of the gradient flow.In section 2 we sketch the derivation of a singular FCH model from a random phase approxima-tion of self-consistent mean field theory, outline the regularization of the singular model to arrive atthe family of FCH models studied here-in, and present the initial data and associated benchmarkproblems. This derivation is relevant as it illuminates the incorporation of the well-stiffness in theFoot 1 and Foot 2 benchmarks that is the initial motivation for this computational study. In section3 we present the second-order adaptive numerical schemes which we use to resolve the benchmarkproblems and highlight the sensitivity of efficiency to choice of implicit terms. In section 4 wepresent an overview of the simulations of each of the five benchmark problems for a fixed localtruncation error, showing where the schemes agree and disagree. In section 5 we contrast the per-formance of the schemes, particularly with respect to accuracy in the far-field of the domain, energydecay, evaluation of the precise critical value for onset of defects, and comparison of time-steppingperformance and computational efficiency. We summarize the performance in section 6.
Self consistent mean-field (SCMF) approach derives density functional models that approximate thebulk interactions of collections of polymers represented by molecular units, [16]. When applied toamphiphilic diblock polymers suspended in a solvent the reduction yields a free energy for the threedensity components φ i for i = A, B, S , which represents the hydrophilic head A and hydrophobictail B of the diblock polymer and the solvent S , respectively. Considering a suspension of n s solventmolecules and n P polymer diblocks each comprised of N A and N B monomers of molecule A and B respectively, [8, 30] used the self-consistent mean field reduction to derive the free energy toa continuum phase-field model. More specifically, for k = { A, B, S } , they introduced the mean4ensities φ A = n P N A | Ω | , φ B = n P N B | Ω | , φ S = n S | Ω | , (2.1)and derived a bilinear approximation to the SCMF free energy expressed in terms of the variancefrom the mean φ k = φ k − φ k , F (2)UD ( φ ) = (cid:88) ij (cid:90) Ω a ij (cid:113) φ i φ j ( D − φ i )( D − φ j ) + b ij (cid:113) φ i φ j + χ ij φ i φ j + δ ij c ij φ i |∇ φ i | dx. (2.2)Here a = ( a ij ), b = ( b ij ), c = ( c ij ) denote material parameters and δ ij is the usual Kroneckerdelta function. Their derivation is similar to [9], with both approaches incorporating long-rangeinteraction terms through D := ( − ∆) , the square-root of the negative Laplacian subject to pe-riodic boundary conditions. The long-range terms describe entropic effects of chain folding andvolume exclusion derived from the interactions of the polymer chains with effective mean fields. Asimilar energy was proposed as a model of a microemulsions of oil, water, and surfactant by [29]who argued directly, and somewhat phenomenologically, from a Landau theory for a scalar density.This scalar model was extended to a more nonlinear one by [18] and [19], who proposed a densitydependence on the coefficients. Uneyama and Doi also proposed a nonlinear extension, [31], fortheir vector model in which the average density φ k was replaced with the local density φ k . Thisextrapolation yielded a family of models that include the Ohta-Kawasaki free energies. In [6] thenonlinear extrapolation approach was modified, first through a shift in dependent variables to thespatially averaged mean-zero density ψ k := D − φ k, , and then by an extrapolation step in whichthe average density φ k was replaced with a slowly varying average density, φ k → φ k (1 + ψ k ) . (2.3)In addition, [6] reduced the three-component model to a scalar field similar to [19] by requiringa point-wise incompressibility, ψ A + ψ B + ψ S = 0, and replacing the global constraint on theA- and B-polymer fractions with the point-wise constraint, φ A /N A = φ B /N B . This yields theparameterization ψ A = ψ B = u/m f , ψ S = − u, in terms of u that takes values in ( − m f , . The resulting model depends upon N P := N A + N B ,the polymer fractions α A = N A /N P and α B = 1 − α A , the polymer-solvent molecular mole fraction m f := n P N P /n S , and the dimensionless parameter ε = lL N / P (cid:28) l of the diblock polymer into a mean-square end-to-end length a single ideal diblock polymerchain expressed as a ratio of the domain length L . The amphiphilicity of the diblock molecules isexpressed in terms of a weighted Flory-Huggins parameter χ w := − α A α B χ AB + α A χ AS + α B χ BS , (2.4)5here for k, m ∈ { A, B, S } the Flory-Huggins parameters χ km > k -monomer and an s -monomer. The value of χ w depends upon thecomposition of the polymer diblock chain, but not on its length. We assume strong hydrophobicityof the tail groups, which implies that χ BS (cid:29) χ AB and is reflected in the strong-hydrophobicityhypothesis ( H sh ) χ w > . (2.5)With these reductions and notation, [6] reduced the Uneyama-Doi bilinear energy (2.2) to thesingular functionalized Cahn-Hilliard (S-FCH) form F S − FCH ( u ) = (cid:90) Ω (cid:0) ε ∆ u − W (cid:48) S ( u ) (cid:1) + P ( u ) dx, (2.6)where the singular potential W S is given by W (cid:48) S ( u ) = 24 E ( u ) − m f N P (ln | − u | + χ w u ) + C . (2.7)For notational convenience we introduced the bulk exclusion term E ( u ) = m f ln (cid:12)(cid:12)(cid:12)(cid:12) um f (cid:12)(cid:12)(cid:12)(cid:12) , (2.8)which records the energy of small densities of polymer. Significantly E tends to −∞ as u tends to − m f and is O ( m f ) elsewhere. The strong hydrophobicity assumption ( H sh ) guarantees that thepotential W (cid:48) S has three zeros b l < b M < b r in the interval ( − m f , W S is defined as theprimitive of W (cid:48) S that has a double zero at u = b l . Consequently W S has two, generically unequal-depth, local minima at u = b l , b r with − m f < b l < < b r < W S ( b l ) = 0 > W S ( b r ). Thefirst derivatives of the well W S are singular at the endpoints u = − m f and u = 1, correspondingto pure solvent and pure polymer phases. The perturbative potential P takes the form P ( u ) := 9 α A α B u E (cid:48) ( u ) − ( W (cid:48) S ( u )) . (2.9)The constant C does not impact the value of the energy for u ∈ H (Ω). It is chosen to minimizethe the perturbative potential P . We draw motivation for the benchmark simulations from the complexity observed in the experi-ments conducted in [22]. In that study the authors prepared well-stirred dispersions of amphiphilicdiblock of Polyethylene oxide (PEO) - Polybutadiene (PB) in water, and allowed the mixture torelax and come to quasi-equilibrium. The weight fraction of polymer was fixed at 1%, and theyconsidered a long and a short polymer chain, characterized by a fixed molecular length of the hy-drophobic PB, with N PB (= N B ) taken as 45 and 170. They varied the aspect ratio α A = N A /N B ,6 esult: network formation in two-compo-nent dilute aqueous solutions of PB-PEO.Previous experiments with relatively lowmolecular weight macromolecular surfac-tants produced three basic structural ele-ments: spheres (S), cylinders (C), andbilayers (B). These zero-, one-, and two-dimensional structures are readily dis-persed in water at low concentration asspherical micelles, wormlike micelles, andvesicles, respectively, mimicking the well-established states of aggregation createdwhen low molecular weight amphiphilesare mixed with water ( ). Here we dem-onstrate the formation of “ Y-junctions, ” which assemble into a dense, three-dimen-sional network (N) accompanied by macro-scopic phase separation. Isolated fragmentsof the network have been evaluated withcryogenic transmission electron microsco-py (cryo-TEM), from which we deduce thatthe Y-junction is preferred thermodynami-cally above a critical surfactant molecularweight at compositions between those as-sociated with the B and C morphologies.These findings are consistent with theoret-ical predictions recently developed ( , )and demonstrated ( , ) with traditionalthree-component (surfactant/oil/water) mi-croemulsions, signaling a fundamentaltransition in self-assembly behavior as theamphiphile molecular weight is increasedbeyond a critical value.Two sets of PB-PEO diblock copolymers(Scheme 1), each containing constant molec-ular weight poly( -butadiene) blocks( M n ,PB ! w PEO ), were syn-thesized with anionic polymerization tech-niques described elsewhere ( ). Fifteencompounds with degree of polymerization N PB !
46 and 0.30 ! w PEO ! N PB !
170 and 0.24 ! w PEO ! M w / M n "
100 to 300nm) films of vitrified aqueous solutions.The results for the 1% mixtures are sum-marized in Fig. 1. With N PB !
46, we findthe classic sequence of dispersed structures (B, C, and S) separated by mixed morphol-ogy regimes (B $ C and C $ S) withincreasing PEO content, in agreement withan earlier report ( ). Increasing the size ofthe hydrophobic block by nearly four times( N PB ! w PEO .However, the most dramatic structuralchanges we found were at compositions be-tween the B and C regions. Decreasing thelength of the PEO block leads to the forma-tion of Y-junctions in the nominally cylindri-cal micelles. Even for w PEO ! w PEO to 0.39 induces striking morphologicalchanges (Fig. 2A). Branched cylindricalloops, branched linear wormlike micelles,and Y-junctions terminated by enlargedspherical caps characterize this representativeimage. Macroscopically, this solution isopaque (milky) but showed no tendency to phase separate even after sitting quiescentlyfor several months.Unlike any of the other specimens exam-ined in this study, dilute mixtures of the w PEO ! O OR N PB N PEO -1H
Scheme 1. Fig. 1.
Morphology dia-gram for PB-PEO in wa-ter (1 wt%) as a func-tion of molecular sizeand composition, where N PB and w PEO are thedegree of polymeriza-tion and weight fractionof the PB and PEOblocks, respectively. Fourbasic structural motifs—bilayers (B), Y-junctions(Y), cylinders (C), andspheres (S)—have beenidentified by cryo-TEM,as illustrated in the mi-crographs, ( A to C ). At N PB ! w PEO from the cylindri-cal condition ( w PEO ! w PEO ! w PEO ! Y (Fig. 2A), isfollowed by network (N)formation and macro-scopic phase separationat w PEO ! N PB ! ) and current experimental results, respectively. Bars (A to C), 100 nm. R E P O R T S on F eb r ua r y , . sc i en c e m ag . o r g D o w n l oaded f r o m examples thr o ugh o ut this paper and elsewhere ; h o w - ever, n o ne are as co nspi c u o us as th o se depi c ted in thisblend. These beadlike de fo rmati o ns occ ur with a c har - a c teristi c peri o di c ity, whi c h appears t o be damped inthe l o ng c entral p o rti o ns of the c ylinders. Sh o rt c ylinderswith o ne and tw o beads c an be seen in Figure 8 A ( sh o rtand l o ng arr o w, respe c tively ) . Cylinders with o ne, tw o ,and three undulati o ns are marked in Figure 8 A . Figure8B sh o w bran c hed p o rti o ns of a c ylindri c al mi c elle withquantized undulati o ns apparently di c tating the arm lengths. Multiple undulating bran c hes in Figure 8Chighlight the l oc alizati o n of su c h f eatures near the j un c ti o ns and ends. These aggregates are stable, asheating the sample t o ° C fo r a f ew days did n o tpr o du c e any n o ti c eable c hange in the assembled m o r - ph o l o gies. A co mparis o n of OB9 - / OB1 - o wn inFigure 4B and OB9 - / OB1 - o ws a transiti o n f r o m mainly spheri c al mi c elles t o mainly c ylindri c al mi c elles o ver a very narr o w co mp o si - ti o n range. B l e n d s a r o un d w PEO * . Re c ently, we dis co verednetw o rk fo rmati o n in aque o us dispersi o ns of OB9 - ( w PEO ) N PB ) ) ( see Figure 9 A) . A s a part of the present study, we attempted t o mimi c thisnetw o rk stru c ture by blending OB9 dibl oc k co p o lymerswith greater ( w PEO > ) and lesser ( w PEO < )co mp o siti o ns. F o r example, we blended equim o lar mix - tures of OB9 - -
6, resulting in an average co mp o siti o n of 〈 w PEO ) 〉 , identi c al t o OB9 -
4. T o o ursurprise this premixed blend sel f- assembles int o ap o tp o urri of deli c ate l oo king o b j e c ts with bilayer, c ylin - dri c al, and co mplex j un c ti o n stru c tural elements, f re - quently mixed within individual m o ieties. These un - usual stru c tures are evident in all the c ry o- TEM imagestaken f r o m this mixture ; Figure 9C displays many of the prevalent f eatures. Perhaps the m o st striking is the oc t o pus (o r j elly f ish )- like mi c elles, whi c h are co mp o sed of a f lat bilayer with pr o truding c ylindri c al mi c ellesal o ng the edges. These oc t o pus - like entities are co mm o n,alth o ugh they occ ur with a variable number of c ylindri -c al arms. Several examples are sh o wn in Figure 10 co ntaining 4, 5, 6, 7, 8, 9, 10, and 14 arms atta c hed t o the f lat c entral p o rti o n. In all these oc t o pus - like ag - gregates the c ylindri c al arms are symmetri c ally dis - tributed al o ng the c ir c um f eren c e of the c entral f latbilayer as is evident in Figures 9C and 10. O cc asi o nally,these o b j e c ts appear t o be fo lded o n a side with the c ylindri c al arms pr o truding f r o m a hemispheri c al bi - layer c ap. Obvi o usly, the co n f inement c reated by the F i g u r e . Cry o- TEM images f r o m mixture OB9 - / OB1 - c ting undulati o ns and distended spheri c al end c aps o nw o rmlike c ylindri c al mi c elles. Sh o rt c ylinders with an undula - ti o n ( sh o rt arr o w ) and tw o undulati o ns ( l o ng arr o w ) and c ylinder ends with o ne, tw o , and three beads are marked co rresp o ndingly in (A) . Panels B and C sh o w bran c hing withquantized undulati o ns f ixing the segment length between j un c ti o ns. Tw o types of hyperb o li c ( saddle ) sur f a c es c hara c ter - ize these m o rph o l o gies ( see Figure 13 ) . F i g u r e . Cry o- TEM mi c r o graphs f r o m three dispersi o ns with identi c al co mp o siti o ns, 〈 w PEO ) 〉 ; the m o le c ular weightdistributi o n br o adens f r o m A t o B t o C. (A) A netw o rk f ragment f r o m OB9 -
4, a single co mp o nent dispersi o n. This pi c ture isrepr o du c ed f r o m re f ( B ) Blend OB9 - / OB9 - A bim o dal distributi o n of co mp o nent co mp o siti o ns ( w PEO ) ) breaksthe netw o rk. ( C ) A br o ader distribti o n ( w PEO ) - / OB9 - ) pr o du c es a variety of m o rph o l o gies in c luding vesi c les,w o rmlike mi c elles, and a new type of hybrid parti c le re f erred t o as an oc t o pus. Tw o of these o b j e c ts, co mprised of c ylindri c al armsradiating f r o m a single bilayer, are evident in this image, o ne with 11 and o ne with 4 arms. Ma c r o m o l ec ul e s, Vo l. 37, N o . 4, 2004 PEO - PB Mi c ellar Dispersi o ns network material spanning the vitrified film,which screens individual features. Neverthe-less, the edges of the object captured in Fig.2B reveal the same general structural ele-ments evident in the smaller fragment pre-sented in Fig. 2C. Decreasing the size of thePEO block (i.e., from w PEO ! N PB !
46 diblock copolymer solu-tions between the B and C states (Fig. 1)failed to uncover such phase behavior ( ).Increasing diblock copolymer molecularweight has another notable consequence. Theequilibrium solubility in water drops expo-nentially with the degree of polymerization,resulting in extremely slow exchange dynam-ics between individual micelles ( ). Hence,fragmentation of the network phase by stir-ring or sonication produces a nonergodic dis-persion of particles. Until buoyancy forcescompact and fuse these particles (smaller par- ticles should survive longer under the actionof thermal motion), each is subject to a lo-calized free-energy optimization.We have exploited this property to pro-duce small, isolated micelles in the 1 wt% w PEO ! w PEO ! ). Because thelocal movement of PB-PEO molecules isunhindered (the glass transition for the un-entangled PB blocks is about – ° C) ( ),certain rearrangements of the tubular struc-ture are feasible subject to the well-estab-lished rules governing microphase separa-tion of monodisperse block copolymers( ). Optimization of the particle free- energywill pit the overall surface tension ( propor-tional to the micelle surface area) againstchain stretching inside (PB) and outside(PEO) the tubular structure ( ). The unifor-mity in core diameters (34 nm on the basis ofthe cryo-TEM images) found within eachmicelle and across all micelles is evidencethat these rules are operative. Accounting forthe recorded shapes and topologies shouldprove a challenge to self-consistent field the-orists working with block copolymers. Evenif we ignore “ isomerizations ” involving topo-logical transitions (i.e., breaking or collaps-ing loops) ( ), Y-junctions can be createdby sprouting a spherical cap on a cylindricalsection with the required polymer drawnfrom loop sections. The prevalence of spher-ically capped Y-junctions in Fig. 2, B and C, Fig. 2.
Cryo-TEM images obtained from 1 wt% aqueous solutions of ( A ) w PEO ! B and C ) w PEO ! Y in Fig. 1 ( w PEO ! w PEO ! Fig. 3. ( A to N ) Representa-tive cryo-TEM images of net-work phase fragments ob-tained after agitating the w PEO ! R E P O R T S on F eb r ua r y , . sc i en c e m ag . o r g D o w n l oaded f r o m Figure 1: (left)
Experimentally observed bifurcation diagram for the morphology of blends of Polyethyleneoxide (PEO) - Polybutadiene (PB) amphiphilic diblock in water. The horizontal axis, w PEO , is the weightfraction of PEO as a percent of the total diblock weight, and the vertical axis denotes the molecular weightsof the PB component of the diblock, fixed at N PB = 45 or 170 (vertical axis). Morphological Complexity is observed for N PB = 170 but not for the shorter N PB = 45 chains. (right) Experimental images fromthe morphological complexity regime showing (top) network structures and (bottom) a mixture of end capsand Y -junction morphology corresponding to regions marked N and C Y in the bifurcation diagram. FromFigures 1 and 2AC of [22], Reprinted with permission from AAAS. characterized by the weight fraction, w P EO , of the amphiphilic PEO component. They recovereda bifurcation diagram, presented in Figure 1 (left), which shows that for the short chains the well-mixed dispersions largely formed codimension one spherical bilayer interfaces, codimension twosolid tubes, or codimension three solid spherical micelles, with some overlap depending upon theaspect ratio. However for α A ∈ (0 . , .
5) the suspensions of long chains form structures that areloaded with defects, such as the network structures and endcaps depicted in Figure 1 (right - topand bottom).The self-assembly of spatially extended morphologies from a relatively dilute suspension canbe viewed as an arrival and a redistribution process. The dispersed amphiphilic molecules aregenerically too dilute to self assemble, but may diffuse until they arrive at localized structure, insertthemselves and lowering their contribution to the system energy by isolating their hydrophobic tailfrom contact with the solvent. Within the FCH model we study, the rate of arrival determines thefinal outcome of this growth phase. The selection mechanism for the end state is delicate, with7any possible outcomes separated by slightly different final energies. This landscape affords anexcellent diagnostic to benchmark the performance of computational tools.To regularize the benchmark problems we make several changes to the initial configurationand to the model. In particular we replace the well-stirred initial dispersion, typically modeledwith random initial data, with a fixed bilayer interface configuration with an asymmetric shapeand an added constant background of amphiphilic diblock that emulates the dispersed reservoir.The asymmetry in the shape seeds the motion against curvature, which in a benchmark problemis best not left to random fluctuations as would be the case for a perfectly circular initial shape.For computational reproducibility we smooth the well, replacing the singular W S , especially itsprinciple singularity at u = − m f , with a smoother well W q with parameter q that regulates issecond derivative. As we shown in Figure 2, the most significant impact on the energy induced bythe increase in the diblock polymer length is an increase in the value of W (cid:48)(cid:48) S in the solvent phase.This raises the energy of a dilute, spatially constant suspension of the amphiphilic diblock molecule,as is consistent with the exposure of a long hydrophobic tail to solvent, [7].Intuitively, both a high density of dispersed diblock polymers or a high energy associated to anisolated diblock molecule correspond to a high rate of absorption of the dispersed polymers ontothe bilayer interface. This rate is a key quantity controlling defect formation. When the arrival rateis slow, the bilayer interface can grow in size to accommodate the new mass. The growth process isadiabatic and has been studied rigorously, [5], deriving a motion against curvature, regularized bya higher order Willmore term that includes surface diffusion. If the rate of arrival increases beyonda critical threshold, then defects, such as pearling, endcaps, and loop formation are observed. Atmoderate rates, a pearling bifurcation can be triggered, the onset of which is well understood withinthe context of the FCH gradient flow, [10]. The pearling can be transient, subsiding as the dilutesuspension of amphiphilic material is consumed. The pearling can also be lead to the formation ofend-cap type defects, essentially micelles that remain connected to the underlying structure fromwhich they emerged. The endcaps form most readily at points of high curvature of the bilayerinterface. The stem of the endcap can grow, forming a long trailing bilayer-type stem and mayultimately reconnect with the initial structure, forming a loop. At yet higher arrival rates thebilayer interface itself may undergo curve splitting – directly forming closed loops and networkstructures. The rich array of possible outcomes, and the wide variety of end-states of the gradientflow, provide an excellent diagnostic of the accuracy of the proposed schemes.Our benchmark problems stimulate the rate of arrival in two ways. In the first set of benchmarkswe fix a smooth bilayer profile and a version of the non-singular well W q with modest derivatives,determined by q = 0, and vary the initial background density of amphiphilic molecules. At alow background level, the sub-critical benchmark, the initial bilayer interface absorbs amphiphilicmaterial and increases its length, however the rate is sufficiently slow that there is no generation8f defects. In the super-critical benchmark the background level is raised to induce the formationof several defects that coalesce and merge over time. In the critical benchmark the aspect ratioparameter η is adjusted to trigger a more dramatic pearling transient within the bilayer interface.Accurate simulations in this benchmark approach but just fail to form an endcap before relaxingback to a smooth bilayer profile as the reservoir of dispersed diblock molecules is depleted. In thesecond set of benchmarks, we return to the initial data and systems parameters of the sub-criticalcase, but increase the value of W (cid:48)(cid:48) q ( b − ) by adjusting the value of q within the well. This increasesthe rate of absorption without adjusting the total amount of material absorbed, and induces defectformation. The FCH free energy with the non-singular well takes the form E FCH ( u ) := (cid:90) Ω (cid:16) ε ∆ u − W (cid:48) q ( u ) (cid:17) − (cid:18) ε η |∇ u | + η W q ( u ) (cid:19) dx. (2.10)The FCH equation is given by the H − gradient flow of E FCH u t = ∆ δ E FCH δu , (2.11)which takes the explicit form u t = ∆ (cid:104) ( ε ∆ − W (cid:48)(cid:48) q ( u ))( ε ∆ u − W (cid:48) q ( u )) − (cid:0) − ε η ∆ u + η W (cid:48) q ( u ) (cid:1) (cid:105) (2.12)and the regularized double well potential is given by W q ( u ) := (cid:20) ( u − b − ) ε (cid:18) − sech( u − b − ε ) (cid:19)(cid:21) (cid:20) ( u − b + ) γ (cid:16) u − b + − b − (cid:17)(cid:21) . (2.13)We fix b ± = ± γ = 0 .
3. The parameter q controls W (cid:48)(cid:48) q ( b − ), asdepicted in Figure 2 (bottom), adding a small foot-like extension to the left-most local minima.To compare the singular and the regularized wells it is convenient to exploit a rescaling of theFCH-SCMF energy that leaves the associated gradient flow invariant: ε → ε √ ν , W S → W S ν , P → Pν , t → ν t. The rescaling of ε is equivalent to a change in domain size L → √ νL. These rescalings simplifythe process of adjusting the regularized well to fit the singular well under changes in the polymerchain length. We take each monomer to have equal weight, equal to the molecular weight of thesolvent. Correspondingly the weight fraction of PEO, w PEO , equals the molar fraction, α A , andthe polymer weight fraction within the solvent reduces to the molar fraction of polymer, m f = n P N P n s = 1100 . W S -1 0 1-505101520 D u2 W S -1 0 1-0.2-0.100.10.2 W q q = 0q = 0.2q = 0.5 -1 0 10246810 D u2 W q W q ''(b - ) = 1.7 W q ''(b - ) = 5.1 W q ''(b - ) = 10.2 q = 0q = 0.2q = 0.5 Figure 2: (top) Graph of scaled, translated singular well W S and its second derivative as recovered byreduction of SCMF for N B = 45 (red) and N B = 170 (blue-dotted). (bottom) Graphs of the regularizedwell, W q and its second derivative for q = 0 , . , . For the benchmark problems that model the long and short chain polymers we take α A = 0 . , m f =0 .
01, and ξ w = 3. For the short-chain polymer benchmark we take N P = 45 and C = 0 . N P = 170 and C = 3 .
0, and rescale the well W S by a factorof ν = 4 .
4. For these parameters W S has a first and singular second derivative at u = − m f = − . W (cid:48) S , is located just to the right of this singularity,at b l = − . b l = − .
01 + 10 − for the long-chain polymers.We rescale and translate u → u − b l ) b r − b l − , so that the left and right wells occur at b ± = ± W q . The resulting W S is presented inFigure 2 (top) and compared to the regularized well W q used in the benchmark simulations. WhileFigure 2 (top-right) suggests that the second derivative of W S is larger for the short chains, in factthe curves cross just to the right of u = −
1, so that for the short chains W (cid:48)(cid:48) S ( −
1) = 607 while forthe long chains W (cid:48)(cid:48) S ( −
1) = 2022 , a ratio of 3 . . The values of η and η for the sub-critical benchmark were determined from a least-square fitof P for the long-chain data to the form of the functionalization terms used in E . For the other10enchmarks the values are reported in Table 1. For the critical case, the value of η and d weretuned to enhance the strength of the pearling transient. Table 1:
Parameters for Benchmark Cases.
Case \ Param q η η d ε γ α m (q) N MassSub-critical 0 1.45 ε ε ε ε ε ε ε ε ε ε To construct the initial data, we fix Ω = [0 , L ] with L = 4 π, N = 256 corresponding to a mesh spac-ing ∆ x ≈ . (cid:8)(cid:0) x ( t ) , y ( t ) (cid:1) , t ≤ t ≤ t (cid:9) ,we construct a region Γ R of uniform width R about Γ, with outer and inner boundaries Γ ± definedby Γ ± = (cid:26)(cid:16) x ( t ) ± y (cid:48) ( t ) s ( t ) R, y ( t ) ∓ x (cid:48) ( t ) s ( t ) R (cid:17) (cid:12)(cid:12)(cid:12) t ≤ t ≤ t (cid:27) , (2.14)where s is arc-length of Γ . We construct the piece-wise constant function φ Γ to be 1 inside Γ R and − F : L (Ω) → C ∞ per (Ω) defined as F [ φ Γ ]( x, y ) = (cid:88) k ,k ∈ I N ˆ φ Γ ( k , k ) exp (cid:0) − λ ( k + k ) (cid:1) exp (cid:18) πiL ( xk + yk ) (cid:19) , where ˆ φ Γ is the discrete Fourier transform of φ Γ interpolated to the N × N mesh with spacing h = L/N, and λ = 7 . × − . With the choice R = 0 . F [ φ Γ ] per unitlength of Γ approximates the mass of an exact bilayer dressing of Γ. For a fixed curve Γ we define φ := F [ φ Γ ], which is fully spatially periodic. For higher resolution simulations, such as N = 512or N = 1024 the corresponding φ N is constructed by sampling φ on the finer grid.For each of the benchmark cases we take the initial data in the form u = φ N + ε d α m (0) , (2.15)where d ∈ R is a parameter that varies in the benchmarks and α m (0) = W (cid:48)(cid:48) q ( b − ) (cid:12)(cid:12) q=0 . The curveΓ is fixed, defined through polar variables as Γ = { ( ρ ( θ ) cos( θ ) + L , ρ ( θ ) sin( θ ) + L ) (cid:12)(cid:12) θ ∈ [0 , π ) } ,where ρ ( θ ) = 3 − ε (cid:0) θ − π
11 ) (cid:1) − ε cos (cid:0) θ − π (cid:1) . u corresponding to N = 256 with this choice of Γ is shown in Figure 3 (right) ford = 0 . The curve Γ is chosen to break any symmetry with the periodic domain and to seed thecurvature growth of the bilayer interface. The mass, m , of the initial data, defined via the relation m := 12 (cid:90) Ω ( u + 1) dx, is reported in Table 1. N = 256N = 512N = 10249.2 9.25 9.3 9.350.71
Figure 3: (left) A 1 D cross-section of the profile φ and finer mesh realizations φ and φ . (right)The initial data u constructed from (2.15) with width R = 0 . u . We use the 2nd order backward differentiation formula (BDF2) to produce the IMEX, PSD, andSAV schemes, and use the solution by the 3rd Adams-Moulton (AM3) scheme as a predictor tocontrol the local error to resolve the benchmark problems described in Section 2.2.
To solve the Cauchy problem: y (cid:48) = F ( y ) , y ( t ) = y , and let us denote the temporal step size k n := t n − t n − , the classical uniform time-step BDF2 scheme is a 2nd order L-stable implicit linearmultistep method 3 y n +1 − y n + y n − k n = F ( y n +1 ) . We assume the variable step size BDF2 scheme has the form ay n +1 + by n + cy n − = F ( y n +1 ) . (3.1)12aylor expanding both sides at t = t n +1 and comparing the coefficients, we identify a = 1 k n +1 + 1 k n +1 + k n b = − k n +1 − k n c = 1 k n − k n +1 + k n . The classical uniform time-step AM3 scheme is a 3rd order (not A-stable) implicit linear multistepmethod reads y n +1 = y n + k n +1 (cid:104) F ( y n +1 ) + 4 F ( y n ) − F ( y n − ) (cid:105) . We assume the general variable step size AM3 scheme has the form y n +1 = y n + (cid:104) ω F ( y n +1 ) + ω F ( y n ) + ω F ( y n − ) (cid:105) . To identify the coefficients { ω i } i =1 , we make the approximation y ( t n +1 ) − y ( t n ) = (cid:90) t n +1 t n F ( y ( t )) dt ≈ (cid:90) t n +1 t n P ( t ) dt, where the quadratic polynomial P ( t ) is the interpolant of F ( y ( t )) at t n − , t n and t n +1 . Introducingthe time-step ratio τ := k n +1 k n , the general variable step size AM3 scheme can be written as y n +1 = y n + k n +1 (cid:104) τ τ F ( y n +1 ) + (3 + τ ) F ( y n ) − τ τ F ( y n − ) (cid:105) . (3.2)Further investigations can be found in [21]. The FCH gradient flow undergoes bifurcations that trigger hidden timescales. As these events occurat unpredictable times, an adaptive approach to time-stepping is required to balance accuracy andefficiency. The FCH equation (2.11) has the form u t = F ( u ), where F ( u ) = ∆ δ E FCH δu . We set a localtruncation error tolerance, σ tol , and the minimal and maximal time-step values k min and k max . Given u := u , t and some final time T , we fix the temporal step size k := k min . Compute thefirst time-step u , t of the FCH equation (2.11) via the 1st order accurate scheme: the locally 2ndorder BDF1-PSD, BDF1-IMEX and BDF1-SAV schemes for the PSD , IMEX and SAV schemesseparately. For successive time steps we assume the index n ∈ N + , with u n − , u n and a giventemporal step-size k n . The adaptive algorithm, based upon [28], initially sets k = k n and thenproceeds as follows. Step 1
Compute the 2nd order accurate main approximation ˜ u n +1 using a scheme with k n , k .13 tep 2 Compute the time step ratio τ = kk n and a 3rd order accurate approximation u p through u p = u n + k (cid:104) τ τ F (˜ u n +1 ) + (3 + τ ) F ( u n ) − τ τ F ( u n − ) (cid:105) , (3.3) Step 3
Calculate the relative error approximation e n +1 := (cid:107) ˜ u n +1 − u p (cid:107) L (cid:107) u p (cid:107) L , Step 4 If e n +1 ≤ σ tol or k = k min , thensave the solution u n +1 = ˜ u n +1 , the temporal step size k n +1 = k , the current time t n +1 = t n + k n +1 ,and finally increase n by 1.Recalculate the time step k = max { k min , min { A dp ( e n +1 , k ) , k max }} , where A dp ( e, k ) := ρ s (cid:16) σ tol e (cid:17) / k. Goto
Step 1 .We take the safety coefficient ρ s = 0 .
9, and k min = 10 − for all simulations. For the IMEX andSAV schemes k max is taken to be ∞ , while for the PSD scheme, the optimal value of k max dependsupon q , as shown in the Table 3. As discussed in [21], to ensure zero-stability for the variable stepsize BDF2 in (3.1), A dp ( e, k ) needs to be bounded from above by (cid:0) √ (cid:1) k . Numerical explorationwith this bound on A dp showed it afforded no significant impact on the benchmark problems. The FCH equation (2.12) can be rewritten as u t = ∆ (cid:104) L IMEX u + N IMEX ( u ) + ε η ∆ u − η W (cid:48) q ( u ) (cid:105) , (3.4)where we introduce the linear positive operator L IMEX := ε ∆ − α m ε ∆ + α m , (3.5)obtained by linearizing δ E FCH δu in (2.12) about u = b − and dropping the small, negative η and η terms. Here α m = W (cid:48)(cid:48) q ( b − ) depends strongly on q. The term N IMEX has zero linearization about u = b − and is given by N IMEX ( u ) := ε (cid:0) α m − W (cid:48)(cid:48) q ( u ) (cid:1) ∆ u + ε ∆ (cid:0) α m u − W (cid:48) q ( u ) (cid:1) + W (cid:48)(cid:48) q ( u ) W (cid:48) q ( u ) − α m u. The resulting 2nd order semi-implicit IMEX scheme is chosen to stabilize the spatially constantbackground state u ≡ b − . To this end we take the dominant linear terms implicit and the remainder,including smaller linear terms involving η and η , explicit, au n +1 + bu n + cu n − = ∆ (cid:104) L IMEX u n +1 + N IMEX ( u ∗ ,n +1 ) + ε η ∆ u ∗ ,n +1 − η W (cid:48) q ( u ∗ ,n +1 ) (cid:105) , (3.6)14here u ∗ ,n +1 := u n + k n +1 k n ( u n − u n − )gives a locally second order extrapolated approximation at time level t n +1 to make sure the schemeis numerically consistent. Now we can isolate and solve u n +1 in (3.6) from (cid:0) a − ∆ L IMEX (cid:1) u n +1 = − bu n − cu n − + ∆ (cid:104) N IMEX ( u ∗ ,n +1 ) + ε η ∆ u ∗ ,n +1 − η W (cid:48) q ( u ∗ ,n +1 ) (cid:105) . (3.7) To PSD attempts to solve (2.12) via a general implicit set up au n +1 + bu n + cu n − = ∆ δ E FCH δu (cid:12)(cid:12)(cid:12) n +1 . (3.8)Given u n − and u n , the equation for u n +1 is defined in terms of the zero of the negative of theresidual, R ( u n +1 ; u n , u n − ) := Π δ E FCH δu (cid:12)(cid:12)(cid:12) n +1 − ∆ − ( au n +1 ) − ∆ − ( bu n + cu n − ) = 0 , (3.9)where Π denotes the linear zero-mass orthogonal projection operator. The preconditioned steepestdescent (PSD ) method solves this nonlinear system (3.9) iteratively through a series of linearsystems. We introduce L PSD , the linearization of (3.9) about the spatially constant state u n +1 ≡ b − and dropping the small η and η terms, L PSD := ε ∆ − α m ε ∆ + α m − a ∆ − . This strictly positive, self-adjoint operator is well-defined on mass-less functions, and preconditionsthe iterative scheme. The search direction d n +1 s at u n +1 s is defined as d n +1 s := −L − R ( u n +1 s , u n , u n − ) . The solution u n +1 is defined as the limit of the sequence { u n +1 s } ∞ s =0 , constructed through theapproximate line search (ALS) recurrence relation u n +10 := u n + k n +1 k n ( u n − u n − ) , (3.10) u n +1 s +1 = u n +1 s + λd n +1 s , s = 0 , , , . . . (3.11)This solver is referred to as the PSD with approximate line search (ALS) algorithm, see [15, 3]. For aprescribed iterative stopping tolerance i tol , the ALS scheme is terminated once (cid:107) d n +1 s (cid:107) L (cid:107) u n +1 s +1 (cid:107) L < i tol . Theparameter λ in (3.11) is the search-step-size. Numerical investigations show that the optimal valueof λ is somewhat sensitive to the value of α m = α m (q) and k . This dependence was determinedby minimizing the average number of PSD iterations for a fixed k over the first 50 temporal steps15 able 2: Dependence of optimal value of search-step-size λ on temporal step size k . k ≤ − · − − · − − · − − · − . . λ for different values of q and k are reported in Table 2. Thevalues used in the simulations were determined by linear interpolation.The iterative stopping tolerance, i tol , impacts the accuracy and computational cost of the PSDscheme. Numerical optimization finds that an optimal choice of i tol is sensitive to both the wellstiffness, q and the local truncation error, σ tol . We determine this relation through the ratio i tol = ν (q) σ tol , and determine an optimal value of ν (q) . This requires balance, as overly small values of i tol leadto excessive iterations that do not improve the scheme’s accuracy. On the other hand i tol must besmall enough to ensure that numerical error from the iterative solver does not pollute the adaptivetime-stepping and does not impede the convergence of the iterative solver at subsequent time-steps. Instructively, the iterative convergence rate was found to depend upon the upper limit, k max ,imposed on the adaptive time-stepping algorithm. This lead to a coupled numerical optimizationstudy, presented in Table 3 which shows the sensitively of iterations numbers upon k max for thethree values of q, and the optimal value of ν . The iteration counts increase considerably with q,while ν decreases exponentially with q. If the upper bound k max is removed then the iterationcount may increase considerably, with associated increase in computational effort. The tuning of k max and ν with q was the most unpredictable element of the optimization process for any of thethree schemes. Table 3:
Dependence of PSD iteration count on ν and k max . Iteration count/1000 ν (q) Value of k max optimal k max .5 The SAV scheme Computational schemes based upon the SAV formulation have been applied to the FCH gradientflow, see [34]. The version presented here is a slight variation. We rewrite the FCH energy functional E FCH ( u ) in (2.10) in the form: E FCH ( u ) = (cid:90) Ω (cid:20) ε u ) − (cid:16) η ζ (cid:17) ε |∇ u | + G ( u ) (cid:21) dx, (3.12)where ζ > G ( u ) := − ε ∆ u (cid:0) W (cid:48) q ( u ) + ζu (cid:1) + 12 ( W (cid:48) q ( u )) − η W q ( u ) . (3.13)The choice of principle linear operator for the SAV scheme is a bit less intuitive than for theIMEX or PSD schemes. We introduce L SAV = ε ∆ + ε ( η + 2 ζ ) ∆ = L + L , (3.14)where the sub-operators are parameter dependent L ( β , β ) = ε ∆ − β α m ε ∆ + β α m , (3.15) L ( β , β ) = ε ( η + 2 ζ ) ∆ + β α m ε ∆ − β α m , (3.16)where α m = α m (q) and the constants β , β ≥ L defines the principle linear implicit terms in the SAV scheme. The default choice for theseparameters is β = 0 and β = 3 . Introducing the auxiliary energy E ( u ) = (cid:90) Ω G ( u ) dx, the FCH energy (3.12) takes the form E FCH ( u ) = 12 ( u, L SAV u ) L (Ω) + E ( u ) . (3.17)For fixed time-steps the SAV scheme is known to be pseudo-energy stable for an auxiliary energy,if the functional E ( u ) can be shown to be uniformly bounded from below over H (Ω), [26]. Thisis achieved by choice of ζ = ζ (q) . Specifically E ( u ) ≥ (cid:90) Ω (cid:16) W (cid:48)(cid:48) q ( u ) + ζ (cid:17) |∇ u | dx + (cid:90) Ω (cid:104)
12 ( W (cid:48) q ( u )) − η W q ( u ) (cid:105) dx, and choosing ζ larger than negative of the minimum of W (cid:48)(cid:48) q , see Figure 2, then the first integral ispositive and we estimate E ( u ) ≥ | Ω | min u (cid:16)
12 ( W (cid:48) q ( u )) − η W q ( u ) (cid:17) > − D , D > η and the form of W q .For the energy splitting approach, we introduce the scalar auxiliary variable r = r ( t ) := (cid:112) E ( u ) + D , then the FCH equation can be rewritten as ∂u∂t = ∆ µ, µ := L SAV u + r (cid:112) E ( u ) + D V [ u ] , (3.18) drdt = 12 (cid:112) E ( u ) + D (cid:90) Ω V [ u ] ∂u∂t dx, (3.19)where V [ u ] = δ E /δu = G (cid:48) ( u ). The SAV scheme takes the form au n +1 + bu n + cu n − = ∆ µ n +1 , µ n +1 = L u n +1 + L u ∗ ,n +1 + r n +1 (cid:112) E ( u ∗ ,n +1 ) + D V [ u ∗ ,n +1 ] , (3.20) ar n +1 + br n + cr n − = (cid:90) Ω V [ u ∗ ,n +1 ]2 (cid:112) E ( u ∗ ,n +1 ) + D (cid:0) au n +1 + bu n + cu n − (cid:1) dx, (3.21)where u ∗ ,n +1 can be chosen as any explicit approximation of u ( t n +1 ) with an error of O ( k ), forinstance, u ∗ ,n +1 = u n + k n +1 k n ( u n − u n − ) . We remark that the r n +1 variable in (3.21) also contribute to the implicit equation for u n +1 . Thefull resolution of u n +1 from (3.20) − (3.21) is presented in [34, 28], but is driven by the inversionof the operator L − a ∆ − . With a fixed time-step k , the SAV scheme is unconditionally energystable for the auxiliary energy E aux ( u n +1 , u n , r n , r n − ) := E FCH ( u n +1 ) + (cid:16) u n +1 − u n , L SAV (cid:16) u n +1 − u n (cid:17)(cid:17) L (Ω) + (cid:0) r n − r n − (cid:1) . That is, for any choice of time-step k >
0, the auxiliary energy has the property E aux ( u n +2 , u n +1 , r n +1 , r n ) ≤ E aux ( u n +1 , u n , r n , r n − ) , n ≥ . See [28] for details on the energy stability properties of SAV schemes.The stabilization parameters makes L a strictly positive operator and play an essential role inthe convergence, accuracy, and efficiency of the SAV scheme. The operator L reduces to L IMEX for the choice β = 2 and β = 1. Figure 4 shows FFT counts for simulations of IMEX and SAVusing the dominant implicit term based on L . Overall the schemes preform well if β + β = 3, withperformance deteriorating dramatically for smaller values and slowly for larger values of this sum.Indeed values of β + β < β = 3 − β , showing that the performance is optimal so long as neither β nor β for the IMEX scheme ( = 3- ) FFT c a ll s Sub-critical: T = 100~150Sup-ritical: T = 250~300Critical: T = 99~149Foot 1: T = 50~100Foot 2: T = 30~40 for the SAV scheme ( =0) FFT c a ll s Sub-criticalSup-criticalCriticalFoot 1Foot 2 Foot 1 Foot 2
Figure 4:
Total FFT calls verses stabilization parameters β and β for IMEX and SAV for each of the 5benchmark simulations. (left) Total FFT calls verses β for IMEX with β = 3 − β . The choice β = 2used in IMEX simulations herein is indicated with black arrow. (right) Total FFT calls as function of β for β = 0 for the SAV scheme. are too small. The choice β = 2 is the default for IMEX . The right panel shows performance ofthe SAV scheme with β = 0 for a range of value of β for each of the five benchmark problems,suggesting that the choice β = 3 yields the optimal computational efficiency over this class ofoperators. We present an overview of the Benchmark simulations for local truncation error σ tol = 10 − , forwhich the PSD scheme is accurate while the IMEX and SAV schemes are borderline accurate.Generically we find that a global L (Ω) discretization error of 3 × − is generically sufficient toensure that that scheme is quantitatively accurate, with the correct numbers, types, and placementsof defects. Table 4: L Error between PSD and IMEX or SAV simulations for each benchmark at final time.
Benchmark PSD/IMEX PSD/SAV T Sub-Critical 8.80E-02 1.46E-01 250Critical 2.64E-02 3.81E-02 250Super-Critical 7.18E-01 7.27E-01 250Foot 1 2.83E-02 6.00E-02 50Foot 2 3.98E-03 9.90E-03 50
The sub-critical benchmark has a low level of dispersed diblock polymer material, controlled bythe parameter d in (2.15), while the relatively mild concavity of W q at u = b − , controlled by α m (0) = W (cid:48)(cid:48) q ( b − ) (cid:12)(cid:12) q=0 , leads to a gentle absorption rate. The bilayer interface profile does not19 igure 5: Simulation of sub-critical benchmark 1 with q = 0, σ tol = 10 − and N = 256 at times T = 10(left)and T = 250(right). All three schemes agree to within L error 1 . × − as reported in Table 4. Red numberon colorbar indicates max u . pearl and remains a simple closed curve from initial data to its final equilibrium shape. As shownin [5], gentle absorption drives motion against curvature, regularized by surface diffusion, whichrelaxes to a curvature driven flow as the background material is depleted. All three schemes arein quantitative agreement, as can be verified by the contour plot comparison in Figure 8 (left) andthe data of Table 4. Figure 6:
Simulation of the critical benchmark with q = 0, σ tol = 10 − and N = 256 at times T = 15(left)and T = 21(right). All three schemes agree to within L error 3 . × − as reported in Table 4. Red numberon colorbar indicates max u . For the critical case the value of η and d are tuned to create a strongly pearled interface anda long pearled transient, lasting roughly from T = 4 to T = 21. The bilayer interface pearlstransiently, forming 21 pearls, whose discrete count generates a thresholding effect that slows the20bsorption of the dispersed amphiphilic polymer as the interface must generate new pearls tolengthen. During the 21-pearl transient period the pearled bilayer interface undergoes a “bicyclechain” meander in which adjacent pearls move in opposite directions, either in towards the center orout towards the boundary of the domain, as can be seen in Figure 6 (left). At time T = 21 the pearlshave reduced in size, and two extra pearls form at the points of highest curvature. The formationof the additional pearls facilitates an absorption of mass. As the background level of amphiphilicmaterial is depleted the rate of absorption slows and the the interface returns to an unpearled state,similar to that depicted in Figure 5 (right) that is able to move freely in a curvature driven motion.No endcap defects are formed in the critical benchmark, and each of the computational schemesare in quantitative agreement. Figure 7: (Top) Simulation of the q = 0 super-critical benchmark with σ tol = 10 − and N = 256 at time T = 50 . The top row presents the PSD simulation (left) and the SAV simulation (right). (Bottom) The PSD(left) and SAV (right) simulations at T = 250. The IMEX and SAV simulations are very similar, and butboth disagree with the PSD simulation by an O (1) in the L norm, as reported in Table 4. Red number oncolorbar indicates max u . PSDIMEXSAV
PSDIMEXSAV p Figure 8:
Contour curves from each of the simulations of each of the three schemes with σ tol = 10 − and N = 256. The level set u = − .
12 for (left) the sub-critical simulation at T = 10 and (right) the super-criticalbenchmark at T = 50. The sub-critical and super-critical benchmarks differ only in the level of the background mate-rial, controlled by the parameter d in (2.15). The elevated value of this parameter in the super-critical benchmark increases the rate of arrival of mass to the interface, exceeding the interface’scapacity to absorb the arriving mass via a curve lengthening flow or by pearl generation. The in-terface undergoes defect generation. For the super-critical benchmark with σ t ol = 10 − the outputfrom the three schemes disagree at leading order. For the PSD scheme the bilayer interface absorbsmaterial from the background and pearls locally at points of high curvature, and then ejects 8 end-cap defects, five of which intersect back with the underlying interface, forming closed loops. Twoof the loops subsequently merge to form an extended loop which grows into a cisternal structurecharacterized by two long parallel interfaces. The IMEX and SAV simulations differ from the PSD ,but agree with each-other. They also produce 8 endcap defects initially, however only four of themsubsequently form closed loops. Two of these loops merge, forming a cisternal structure, howeverthere are two small endcaps in the IMEX and SAV simulations, in contrast to the one small endcapin the PSD simulation. At longer times the cisternal region grows, consuming structures and attime T = 250 it leaves one loop, one long endcap, and one short endcap in all simulations – howeverin the SAV and IMEX simulations the distance between cisternal region and small loop is signifi-cantly longer than in PSD simulation. Figure 8 shows the levels sets corresponding to u = − . σ tol = 10 − and N = 256, showing theiragreement in the sub-critical benchmark and their disparity in the super-critical benchmark. Inthe super-critical benchmark the higher rate of absorption driven by the higher initial backgroundlevel of u produces dynamic choices associated to endcap formation that require greater accuracythan the linearly implicit schemes can achieve at σ tol = 10 − . If σ tol is reduced to 10 − , thenPSD simulations does not change quantitatively, while the SAV and IMEX simulations move into22uantitative agreement with the PSD scheme. Figure 9:
Value of u − b − at center point (solid) and corner point (dashed) of computational domain forthe sub-critical (left) and super-critical (right) benchmarks with σ tol = 10 − . Horizontal axis is log of time. The value of u in the far field, away from the interfacial structure, is asymptotically constant atequilibrium and has been shown to be a key bifurcation parameter for the onset of pearling, [24, 25].Faithful resolution of this value is essential to an accurate simulation. Figure 9 traces the evolutionof the value of u − b − at the domain center (solid lines) and domain corner (dashed lines) for thethree simulation strategies. For the sub-critical simulation no defects are formed and the far-fieldvalues of u relax to a tight range of equilibrium values over the time frame T = 75 ∼ T = 200 , , , , and 825 for the PSD scheme. Conversely the linear implicit IMEX and SAVschemes the background levels are in close agreement, recording excursions T = 175 , , , and875 , which differ in both timing and in number of events from the more accurate PSD simulation. The Foot 1 and sub-critical benchmarks, are identical in initial data and parameters with theexception of the value of the concavity of the well W q , controlled by the parameter q. For Foot1 we take q = 0 . α m (q) = W (cid:48)(cid:48) q ( b − ), as depicted in Figure 2. Thisadjustment raises the energy associated to small, spatially uniform values of u , thereby increasingthe rate of absorption of material from the bulk. Although the total amount of material in thebackground is the same in both benchmarks, the increased absorption rate in the Foot 1 benchmarkleads to defect formation. This result is captured by all three schemes with quantitative accuracy, asshown in Table 4. In Figure 10 (left) the pearling and defect formation are visible in the lower-rightof the bilayer interface already at time T = 1 .
5. At time T = 50 the simulations produce six closedloops place roughly symmetrically around the bilayer interface. This structure is quasi-stable, but23 igure 10: Simulation of Foot 1 benchmark with q = 0 . σ tol = 10 − and N = 256 at times T = 1 . T = 50(right). All three schemes agree to within L error 6 × − as reported in Table 4. Red numberon colorbar indicates max u . eventually evolves onto a double-sheeted bubble similar to that depicted in the right-most panel ofthe top row of Figure 18. Figure 11:
Simulation of the Foot 2 benchmark with q = 0 . σ tol = 10 − at times T = 1 (left) and T = 50(right) for N = 512. All three schemes agree to within L error 1 × − as reported in Table 4. Red numberon colorbar indicates max u . The Foot 2 and sub-critical benchmarks have an identical setup with the exception of the valueof q, which is taken as q = 0 . α m = W (cid:48)(cid:48) q ( b − ) (cid:12)(cid:12) q=0 . significantly increases the energy penalty associated to dispersedamphiphilic material. As a consequence its rate of absorption onto the bilayer interface increases,inducing a curve-splitting bifurcation in which the bilayer interface splits directly in two, as shownin Figure 11 (left) at T = 1. All three schemes agree qualitatively on the 512 ×
512 mesh, producingfour loops and two double loop. Grid refinement showed that the N = 256 grid was insufficient toproduce accurate results. Further grid refinement to N = 1024 yields quantitative agreement with24he N = 512 simulations. The large value of W (cid:48)(cid:48) q ( b − ) for q = 0 . (cid:113) W (cid:48)(cid:48) q ( b − ) (cid:46) ε ,which is significantly greater for q = 0 .
5, necessitating the higher spatial resolution.
Table 5: L Grid Refinement Error. N
256 / 512 512 / 1024Sub-Critical 6.218E-04Super-Critical 2.589e-04Critical 3.827E-04Foot 1 8.502E-02Foot 2 1.008 5.762E-04
Table 6: L Temporal Convergence Error and Rates. The error is determined by comparison to PSD witha fixed temporal step size k = 10 − and i tol = 10 − . IMEX PSD SAVfixed k L Error Rate L Error Rate L Error Rate2 × − × − × − × − − × − − × − − × − − × − − × − − × − u − b − evaluated at the domain center (solid) anddomain corner (dashed), are presented for the Foot 2 benchmark in Figure 12 (left). It has severalnotable differences from the sub-critical benchmark presented in Figure 9 (left). The most salientdistinction is that the large value of α m (0 .
5) greatly increases the temporal rate of absorption ofamphiphilic material from the background. For the Foot 2 benchmark the background state beginsto achieve its equilibrium value at T = 1 and is fully equilibrated around T = 7 ∼ . This is roughly10-15 times faster than the relaxation for the q = 0 sub-critical benchmark, depicted in Figure 9.The only difference in the two benchmarks is the increase in concavity of the well W q near u = b − ,reflected in Figure 2(bottom-right). 25 Computational Accuracy and Efficiency
The three schemes presented here are 2nd order accurate, as verified by the convergence studypresented in Table 6. Nevertheless, the performance of the schemes is not equally accurate norefficient, particularly as the nonlinear stiffness parameter q is increased. We discuss the relation ofaccuracy to energy decay, global discretization error, and computational efficiency.
A major feature of gradient schemes is the decay of the overall system energy. Much attention hasbeen given to the construction of gradient stable schemes for which energy decay is unconditionalwith respect to the temporal step-size. However in gradient flows that generate a rich variety ofstructures issues of accuracy move to the forefront and energy decay ideally becomes a consequenceof accuracy. For the super-critical benchmark, the various competing outcomes have only marginallydifferent energies and considerable accuracy is required for a scheme to differentiate between theavailable options. As shown in Figure 13 (left), with σ tol = 10 − for each of the 5 benchmarksthe energy decay behavior is very similar and decays uniformly. There are however importantdifferences. As the middle inset shows, for the super-critical benchmark the energy trace for theIMEX and SAV simulations are almost indistinguishable, while the PSD simulation shows visibledeviations in system energy, at roughly a 1% relative error. In these simulations the PSD schemeis more accurate and the linear-implicit schemes produce output that is erroneous to leading orderin the L norm. The differences in energy decay, and solution u , are largely erased for IMEX andSAV when σ tol is reduced to 10 − . The second inset shows detail of the Foot 1 benchmark. In thiscase the PSD and IMEX schemes have quantitative agreement, while the SAV scheme produces asimulation whose energy proceeds downhill too quickly at T = 650 . This leads to an output thatis time-shifted forward by about 200 units. Again, these errors in energy decay are eliminatedby reducing σ tol to 10 − . These features emphasize that system energy can be a poor proxy foraccuracy, and that energy decay is generally a minor benchmark for a gradient flow.The time-stepping profiles for the IMEX and SAV schemes are remarkably similar, and differ inimportant ways from that of the PSD scheme. As shown in Figure 12 (right), the PSD genericallytakes the largest time step-sizes, and typically hits the maximum step-size ceiling k max shortlyafter the resolution of the initial transient. This ceiling is required to insure the convergence ofthe nonlinear iterative scheme and to optimize its performance as measured by FFT per time unit.This value is smaller for the stiffer Foot 2 benchmark than for the super-critical benchmark. Indeedthe time-step profile for PSD is largely equivalent for the super-critical and the Foot 2 benchmarks,until it hits the lower value of k max for the Foot 2 benchmark. This is in contrast to the IMEXand SAV profiles which are different for the two benchmark problems, but largely agree with each26 -4 -2 -7 -6 -5 -4 -3 -2 -1 PSD: Super-criticalIMEX: Super-criticalSAV: Super-criticalPSD: Foot 2IMEX: Foot 2SAV: Foot 2
Figure 12: (left) Value of u − b − at center point (solid) and corner point (dashed) of the computationaldomain for the q = 0 . σ tol = 10 − and N = 512 . (right) Evolution of the adaptivetemporal step-size on a log-log scale for each of the three schemes for the super-critical benchmark (solid)and the q = 0 . other. Each of the schemes has swings in step size of roughly one order of magnitude during thevarious defect generation and merging events that occur after the initial transient. The step sizesfor the IMEX and SAV schemes are generically smaller than those for PSD , by as much as twoorders of magnitude for the stiffer Foot 2 benchmark. However this is offset by the growing numberof iterations required for solving the stiffer nonlinear system in this problem. -0.8-0.6-0.4-0.2 Sub: PSDSub: IMEXSub: SAVCri: PSDCri: IMEXCri: SAVSup: PSDSup: IMEXSup: SAVFoot 1: PSDFoot 1: IMEXFoot 1: SAVFoot 2: PSDFoot 2: IMEXFoot 2: SAV -0.79-0.78-0.77-0.76-0.75-0.74-0.73-0.72 Sup: PSD, = 10 -5 Sup: IMEX, = 10 -5 Sup: IMEX, = 10 -6 Sup: SAV, = 10 -5 Sup: SAV, = 10 -6
500 1000 1500-0.72-0.715-0.71-0.705-0.7-0.695
Foot 1: PSD, =10 -5 Foot 1: IMEX, =10 -5 Foot 1: SAV, =10 -5 Foot 1: SAV, =2*10 -6 Figure 13: (left) System energy verses time on a semilog-x scale for each of the five benchmark problemsfor each scheme with σ tol = 10 − on a semilog-x scale. The boxed insets for the super-critical (middle) andfoot 1 (right) benchmarks show more detail and include results for IMEX and SAV with σ tol = 10 − . An excellent proxy for accuracy is to determine the lowest (critical) value of the backgroundlevel, as measured by the initial data parameter d in (2.15), at which the a defect is generatedwithin the flow. The onset of a defect is easily detected through the maximum value of u , as themaximum value of the bilayer profile is determined by the largest zero of W q , for (2.13) about u = 0 . W q , with maximum values close to u = 0 . . Tracking the temporalevolution of max u yields a strong dichotomy. We fixed the parameters as in the critical benchmarkproblem but slightly adjusted the value of d to modify the amount of amphiphilic material inthe bulk. The critical d value, reported in Table 7 depends upon the local truncation error, butconverges to a common value of d = 0 . σ tol . Indeed the correct critical valueis achieved by the PSD scheme at σ tol = 10 − but requires reduction of σ tol to 10 − for IMEX andSAV . The time evolution of max( u ( · , t )) for the critical benchmark parameters but seven differentvalue of d is presented in Figure 14. t m a x ( u ( ,t )) d = 0.7522, max(u(100)) = 0.3566d = 0.7523, max(u(100)) = 0.3566d = 0.7524, max(u(100)) = 0.3566d = 0.7525, max(u(100)) = 0.3566d = 0.7526, max(u(100)) = 0.7445d = 0.7527, max(u(100)) = 0.7423d = 0.7528, max(u(100)) = 0.7407 Figure 14:
Running value of max u from the PSDscheme for the critical benchmark problem with d =0 . , ..., . . σ tol = 10 − .The onset of a defect occurs at the critical value d = 0 . . Table 7:
The dependence of the crit-ical value of d in (2.15) upon σ tol foreach scheme. σ tol PSD IMEX SAV10 − − − − The definitive measure of accuracy is to compute the global discretization error of a simulation asmeasured against a known highly accurate answer. To produce these highly-accurate solutions weconducted a spatial grid refinement study for each benchmark problem and each computationalscheme. For all but the stiffest Foot 2 benchmark increasing the grid from N = 256 to N = 512produced consistent results, with solution differences reported in Table 5. We present results onlyto the accuracy determined within this grid refinement study. Specifically the highly accuratesimulations were calculated with the PSD scheme with σ tol = 10 − for q = 0, and with σ tol =3 × − for q = 0 .
2. For q = 0 .
5, the IMEX scheme with σ tol = 3 × − was used. Theoutput of these simulations are taken as the highly accurate simulation against which others arecompared. For all three schemes, sufficient refinement of σ tol lead to a global error that waswithin the anticipated accuracy of the scheme. Indeed our computations find that a global L (Ω)discretization error of 3 × − is generically sufficient to ensure that that scheme is quantitativelyaccurate, with the correct numbers, types, and placements of defects.28 -9 -8 -7 -6 -5 -4 -3 -2tol -4 -3 -2 -1 G l oba l L e rr o r FFT c a ll s PSDIMEXSAV -4 -3 -2 -1 Global L error10 FFT c a ll s PSDIMEXSAV when L error = 3e-2PSD: max = 1.1e-3IMEX: max = 3.1e-6SAV: max = 1.4e-6 Figure 15: (left: blue y-axis and lines) Global L error verses σ tol at T = 250 for the super-critical bench-mark as measured by comparison to the most accurate solution. (left: red y-axis and lines) Computationalcost verses σ tol as measured by total number of FFT calls. (right) Computational cost verses global L error,plotted parametrically in σ tol . We measure the computational efficiency of the three schemes in two ways. First as globaldiscretization error, G tr , verses σ tol , and then more meaningfully as global discretization errorverses FTT calls. This latter is euphemistically referred to as the dollars-per-digit metric. Thefirst result, presented in Figure 15 (left), shows the decay in global L error with decreasing σ tol .The blue curves, corresponding to the left (blue) vertical axis, show that IMEX , PSD , andSAV all improve in global accuracy with decreasing σ tol . For the super-critical benchmark thelinear-implicit IMEX and SAV schemes are inaccurate for σ tol > × − and then have globaldiscretization errors that decay linearly on a log-log plot, corresponding to a global discretizationerror roughly proportional the the 2 / σ tol < × − , but its global accuracy at first improves sub-linearly with σ tol on thelog-log scale before setting into the 2 / − / σ tol acting as a parameterization of the curve. This is shown in Figure 15(right). In this plot, the lowest curve attains the desired global discretization error with the leastcomputation cost.Setting G tr = 3 × − as an acceptable upper limit, we find that the three schemes achieve thisglobal tolerance at comparable computational costs that correspond to disparate local truncationerrors. The IMEX scheme is the most efficient, hitting the global accuracy mark with 1 . × FFT29 -9 -8 -7 -6 -5 -4 -3tol -4 -3 -2 -1 G l oba l L e rr o r FFT c a ll s PSDIMEXSAV -5 -4 -3 -2 -1 Global L error10 FFT c a ll s Foot 1: PSDFoot 1: IMEXFoot 1: SAVFoot 2: PSDFoot 2: IMEXFoot 2: SAV
Foot 1: when L error = 3e-2PSD: max = 7.2e-5IMEX: max = 1.1e-5SAV: max = 3.4e-6 Foot 2: when L error = 3e-2PSD: max = 0.07IMEX: max = 3.2e-4SAV: max = 7.6e-5 Figure 16: (left: blue y-axis and lines) Global L error verses σ tol at T = 50 for the q = 0 . σ tol as measured by total number of FFT calls. calls at σ tol = 3 × − , while PSD does so with 2 × FFT calls at a much lower σ tol = 10 − , andSAV at 2 . × FFT calls and σ tol = 10 − . However the efficiency of the PSD code decays relativelyrapidly with global error above this acceptable upper limit, recovering at very small global error.The overall result is a large interval in which the linear-implicit schemes outperform PSD . For thestiffer Foot 1 benchmark simulations with q = 0 . . ∼ . / N = 512 . Conversely, the nonlinear-implicitPSD requires significantly more effort with increasing q as the iteration count in the nonlinear solverincreases significantly. The minimal cost for SAV to achieve the acceptable global discretizationerror is roughly twice that of IMEX , while for Foot 1 PSD is most efficient at a lower global errorof 9 × − , for which it requires 10 times the effort of IMEX and for Foot 2 PSD requires 40times the computational effort of IMEX . It is worth noting that PSD is more comparatively moreefficient at lower global error; indeed it requires only 3 and 7 times the computational effort ofIMEX to achieve an error of 9 × − for the Foot 1 and Foot 2 benchmarks, respectively.In Figure 17 the temporal trace of the global error is plotted for local truncation errors of σ tol = 10 − , − , and 10 − . In all cases the PSD is the most accurate, generically by an order30f magnitude at the same local truncation error. However the accuracy for PSD increases onlymodestly with decreasing σ tol while SAV and IMEX schemes have more significant improvements.For the sub-critical benchmark the global error accumulates slowly in each of the schemes asthe shape of the interface evolves and inaccuracies in its location accumulate. For the super-critical benchmark the error has peaks at each of the major defect merging events that occur at t = 50 , , , . These peaks reflect the impact of slight timing errors in the defect mergingevents and in the spatial structure of the merging transient. Each scheme shows about a half-orderof magnitude loss of accuracy during the merging that is recovered afterwards. This holds exceptfor the SAV and IMEX schemes with σ tol = 10 − which are both insufficiently accurate to capturethe correct sequencing of the defect evolution. -3 -2 -1 PSD: = 10 -5 PSD: = 10 -6 PSD: = 10 -7 IMEX: = 10 -5 IMEX: = 10 -6 IMEX: = 10 -7 SAV: = 10 -5 SAV: = 10 -6 SAV: = 10 -7 -3 -2 -1 PSD: = 10 -5 PSD: = 10 -6 PSD: = 10 -7 IMEX: = 10 -5 IMEX: = 10 -6 IMEX: = 10 -7 SAV: = 10 -5 SAV: = 10 -6 SAV: = 10 -7 Figure 17:
Time evolution of the global L error between output of the three schemes and the highlyaccurate solution for σ tol = 10 − , − , and 10 − for (top) the sub-critical benchmark and (bottom) thesuper-critical benchmark. We have demonstrated that the morphological complexity that develops within the gradient flows ofthe FCH energy requires accuracy for faithful representation. More significant than energy decay is31nergy accuracy, indeed small errors in evaluation of energy of different configurations can generatealternate temporal evolutions resulting in errors that grow to become leading order. The impact ofthis is magnified as the nonlinear stiffness in the model is increased. The nonlinear solve required inthe more strongly implicit PSD approach tends to raise the overall accuracy of the scheme, and forless-stiff forms of the model this compensates for the increased computational effort required for theiterative solver. The result is that the linear-implicit and nonlinear-implicit models are comparable.However for the more nonlinearly stiff versions of the model, the linear-implicit schemes require notuning and experience only modest decline in efficiency, while the nonlinear-implicit PSD requirestuning of the error tolerance and maximum time-step parameters to optimize its performance.Despite this tuning the efficiency of the PSD scheme falls behind the linear-implicit schemes bya factor that is comparable to the increase in stiffness, as measured by the left-well concavity α m (q) = W (cid:48)(cid:48) q ( b − ) . Within the linear-implicit schemes the performance of the IMEX and the SAV schemes arealmost indistinguishable. Their global accuracy as a function of local truncation error are almostidentical. The only discrepancy lies in the computational effort which is routinely a factor of twolarger for the SAV scheme. This is likely a result of the extra steps required to resolve the largerSAV system of equations. Beyond the guarantee of the decay of the associated pseudo-energy, itis difficult to identify a feature in the SAV scheme in which it improves upon the simpler IMEXapproach. Far and away the most important feature of the linear-implicit schemes is selecting aproper linear term for the implicit step. Given the theoretical understanding of the importance ofthe background state, that is of the value of u away from non-trivial structures, it is reasonable andefficient to use the linearization about the spatially constant state u ≡ b − . We generalized this to thefamily of operators presented in (3.15), and found that and choice of β + β ≈ β = 2 and β = 1 corresponds to the linearization about the spatiallyconstant background state. These constant coefficient linear operators are trivially inverted in thespatially periodic setting considered herein. It certainly may not be the case that such a convenientand efficient linear-implicit operator is available in all systems.As a final demonstration of the complexity possible within the FCH gradient flow, we presenta series of computations that show a putative equilibrium state resulting from the gradient flow ofthe initial data from the super-critical benchmark, see Figure 18. The only variation is in the valueof the parameter η , which represents the aspect ratio of the amphiphilic molecule. The decreasingvalues of η correspond to the increasing values of w PEO in the horizontal axis of the experimentalbifurcation diagram presented in Figure 1. Perhaps the fundamental contribution of this numericalstudy the suggestions that the shapes of the final structures produced in these casting problems arenot uniquely determined by the properties and densities of the molecules they are composed of, butalso depend upon the history of the morphology. Once defects are induced by transient flow, they32ecome an integral part of the energy landscape and can snag the gradient flow at a rich varietyof local minima. These gradient flow transients form an intriguing phylogenesis, whose resolutionrequires significant accuracy. network material spanning the vitrified film,which screens individual features. Neverthe-less, the edges of the object captured in Fig.2B reveal the same general structural ele-ments evident in the smaller fragment pre-sented in Fig. 2C. Decreasing the size of thePEO block (i.e., from w PEO ! N PB !
46 diblock copolymer solu-tions between the B and C states (Fig. 1)failed to uncover such phase behavior ( ).Increasing diblock copolymer molecularweight has another notable consequence. Theequilibrium solubility in water drops expo-nentially with the degree of polymerization,resulting in extremely slow exchange dynam-ics between individual micelles ( ). Hence,fragmentation of the network phase by stir-ring or sonication produces a nonergodic dis-persion of particles. Until buoyancy forcescompact and fuse these particles (smaller par- ticles should survive longer under the actionof thermal motion), each is subject to a lo-calized free-energy optimization.We have exploited this property to pro-duce small, isolated micelles in the 1 wt% w PEO ! w PEO ! ). Because thelocal movement of PB-PEO molecules isunhindered (the glass transition for the un-entangled PB blocks is about – ° C) ( ),certain rearrangements of the tubular struc-ture are feasible subject to the well-estab-lished rules governing microphase separa-tion of monodisperse block copolymers( ). Optimization of the particle free- energywill pit the overall surface tension ( propor-tional to the micelle surface area) againstchain stretching inside (PB) and outside(PEO) the tubular structure ( ). The unifor-mity in core diameters (34 nm on the basis ofthe cryo-TEM images) found within eachmicelle and across all micelles is evidencethat these rules are operative. Accounting forthe recorded shapes and topologies shouldprove a challenge to self-consistent field the-orists working with block copolymers. Evenif we ignore “ isomerizations ” involving topo-logical transitions (i.e., breaking or collaps-ing loops) ( ), Y-junctions can be createdby sprouting a spherical cap on a cylindricalsection with the required polymer drawnfrom loop sections. The prevalence of spher-ically capped Y-junctions in Fig. 2, B and C, Fig. 2.
Cryo-TEM images obtained from 1 wt% aqueous solutions of ( A ) w PEO ! B and C ) w PEO ! Y in Fig. 1 ( w PEO ! w PEO ! Fig. 3. ( A to N ) Representa-tive cryo-TEM images of net-work phase fragments ob-tained after agitating the w PEO ! R E P O R T S on F eb r ua r y , . sc i en c e m ag . o r g D o w n l oaded f r o m network material spanning the vitrified film,which screens individual features. Neverthe-less, the edges of the object captured in Fig.2B reveal the same general structural ele-ments evident in the smaller fragment pre-sented in Fig. 2C. Decreasing the size of thePEO block (i.e., from w PEO ! N PB !
46 diblock copolymer solu-tions between the B and C states (Fig. 1)failed to uncover such phase behavior ( ).Increasing diblock copolymer molecularweight has another notable consequence. Theequilibrium solubility in water drops expo-nentially with the degree of polymerization,resulting in extremely slow exchange dynam-ics between individual micelles ( ). Hence,fragmentation of the network phase by stir-ring or sonication produces a nonergodic dis-persion of particles. Until buoyancy forcescompact and fuse these particles (smaller par- ticles should survive longer under the actionof thermal motion), each is subject to a lo-calized free-energy optimization.We have exploited this property to pro-duce small, isolated micelles in the 1 wt% w PEO ! w PEO ! ). Because thelocal movement of PB-PEO molecules isunhindered (the glass transition for the un-entangled PB blocks is about – ° C) ( ),certain rearrangements of the tubular struc-ture are feasible subject to the well-estab-lished rules governing microphase separa-tion of monodisperse block copolymers( ). Optimization of the particle free- energywill pit the overall surface tension ( propor-tional to the micelle surface area) againstchain stretching inside (PB) and outside(PEO) the tubular structure ( ). The unifor-mity in core diameters (34 nm on the basis ofthe cryo-TEM images) found within eachmicelle and across all micelles is evidencethat these rules are operative. Accounting forthe recorded shapes and topologies shouldprove a challenge to self-consistent field the-orists working with block copolymers. Evenif we ignore “ isomerizations ” involving topo-logical transitions (i.e., breaking or collaps-ing loops) ( ), Y-junctions can be createdby sprouting a spherical cap on a cylindricalsection with the required polymer drawnfrom loop sections. The prevalence of spher-ically capped Y-junctions in Fig. 2, B and C, Fig. 2.
Cryo-TEM images obtained from 1 wt% aqueous solutions of ( A ) w PEO ! B and C ) w PEO ! Y in Fig. 1 ( w PEO ! w PEO ! Fig. 3. ( A to N ) Representa-tive cryo-TEM images of net-work phase fragments ob-tained after agitating the w PEO ! R E P O R T S on F eb r ua r y , . sc i en c e m ag . o r g D o w n l oaded f r o m Th e cor e dim en sion s of t h ese m icelles a r e iden t ica l t ot h ose det er m in ed for t h e n ea t copolym er solu t ion s (n otsh own for br evit y). H owever , pr em ixed solu t ion s of OB9-4/OB9-12 self-a ssem ble in t o lon g cylin dr ica l m icelleswit h fr equ en t ly pla ced Y-ju n ct ion s a lon g wit h a fewsph er ica l m icelles, a n in t er m edia t e m or ph ology t h a t isr ou gh ly con sist en t wit h t h e sin gle-com pon en t beh a viora t a com pa r a ble a ver a ge com posit ion 〈 w PEO 〉 ) w P E O *) a n d t h ose n ot n ea r w P E O *. B le n d s o u t s i d e w P E O *. In F igu r e 6 we pr esen t cr yo-TE M m icr ogr a ph s obt a in ed fr om t h r ee pr em ixed blen dsof OB9-4/OB9-12 a t differ en t m ixin g r a t ios. N ea t dis-per sion s of OB9-4 a n d OB9-12 for m n et wor k a n dsph er ica l m icelles, r espect ively. An equ im ola r m ixt u r eof t h es e d iblock cop olym er s s elf-a s s em bles in t o s h or tcylin dr ica l m icelles t h a t coexist wit h sph er ica l m icelles(F igu r e 6A), qu a lit a t ively m im ick in g t h e in t er m edia t em or ph ology iden t ified in F igu r e 2 bet ween t h e con st it u -en t diblock copolym er s. In cr ea sin g t h e OB9-4 con t en t(2:1, F igu r e 6B) r esu lt s in a decr ea se in t h e fr a ct ion ofm a t er ia l t h a t for m s sph er ica l m icelles a n d t h e develop- F ig u re 4.
Cr yo-TE M im a ges obt a in ed fr om post m ixed (A) a n dpr em ixed (B) sa m ples of t h e OB9-1/OB1-3 blen d. Th is com -pa r ison illu st r a t es t h e n on er godic n a t u r e of polym er ic m icellesfor m ed fr om diblock copolym er s wit h differ en t m olecu la rweigh t s. Th e sm a ll a n d la r ge m icelles in (A) fa iled t o m ix a ft ersever a l m on t h s.
Ta b le 3. Mi c e lla r D i m e n s i o n s o f Va ri o u s Mo rp h o lo g i e s m or ph ology M ncore (g/m ol) a d (n m ) b sph er es 2500 (
40 18. ( (
200 46. ( (
204 38. ( cylin der s 2500 (
40 14. ( (
200 34. ( (
204 25. ( bila yer s 2500 (
40 8. ( (
200 21. ( (
204 15. ( a E r r or in cor e m olecu la r weigh t is obt a in ed fr om gel per m ea t ionch r om a t ogr a ph y (equ ipped wit h ligh t sca t t er in g device) a n a lysis. b Micella r dim en sion , d , for sph er e a n d cylin der is t h e dia m et era n d for bila yer is t h e wa ll t h ick n ess. F i g u re 5.
Cr yo-TE M im a ges fr om post m ixed (A) a n d pr e-m ixed (B) sa m ples of OB9-4/OB9-12 (2:1) blen d illu st r a t in gn on er godicit y of polym er ic m icelles wit h iden t ica l cor e m olec-u la r weigh t polym er s.
M acrom olecu les, V ol. 37, N o. 4, 2004
P E O - P B Micella r Disper sion s
The core dimensions of these micelles are identical tothose determined for the neat copolymer solutions (notshown for brevity). However, premixed solutions of OB9-4/OB9-12 self-assemble into long cylindrical micelleswith frequently placed Y-junctions along with a fewspherical micelles, an intermediate morphology that isroughly consistent with the single-component behaviorat a comparable average composition 〈 w PEO 〉 ) w PEO *) and those not near w PEO *. Ble n ds ou tside w PEO *. In Figure 6 we present cryo-TEM micrographs obtained from three premixed blendsof OB9-4/OB9-12 at different mixing ratios. Neat dis-persions of OB9-4 and OB9-12 form network andspherical micelles, respectively. An equimolar mixtureof these diblock copolymers self-assembles into shortcylindrical micelles that coexist with spherical micelles(Figure 6A), qualitatively mimicking the intermediatemorphology identified in Figure 2 between the constitu-ent diblock copolymers. Increasing the OB9-4 content(2:1, Figure 6B) results in a decrease in the fraction ofmaterial that forms spherical micelles and the develop- Figure 4.
Cryo-TEM images obtained from postmixed (A) andpremixed (B) samples of the OB9-1/OB1-3 blend. This com-parison illustrates the nonergodic nature of polymeric micellesformed from diblock copolymers with different molecularweights. The small and large micelles in (A) failed to mix afterseveral months.
Table 3. Micellar Dimensions of Various Morphologies morphology M ncore (g/mol) a d (nm) b spheres 2500 (
40 18. ( (
200 46. ( (
204 38. ( cylinders 2500 (
40 14. ( (
200 34. ( (
204 25. ( bilayers 2500 (
40 8. ( (
200 21. ( (
204 15. ( a Error in core molecular weight is obtained from gel permeationchromatography (equipped with light scattering device) analysis. b Micellar dimension, d , for sphere and cylinder is the diameterand for bilayer is the wall thickness. Figure 5.
Cryo-TEM images from postmixed (A) and pre-mixed (B) samples of OB9-4/OB9-12 (2:1) blend illustratingnonergodicity of polymeric micelles with identical core molec-ular weight polymers.
Macrom olecules, Vol. 37, N o. 4, 2004
PEO - PB Micellar Dispersions e x a m p l e s t h r ou g h ou t t h i s p a p e r a n d e l s e w h e r e ; h ow -e v e r , n on e a r e a s con s p i cu ou s a s t h os e d e p i ct e d i n t h i sb l e n d . T h e s e b e a d l i k e d e for m a t i on s occu r w i t h a ch a r -a ct e r i s t i c p e r i od i ci t y , w h i ch a p p e a r s t o b e d a m p e d i nt h e lon g ce n t r a l p or t ion s of t h e cylin d e r s . S h or t cylin d e r sw i t h on e a n d t w o b e a d s ca n b e s e e n i n F i g u r e 8 A (s h or ta n d l on g a r r ow , r e s p e ct i v e l y ). C y l i n d e r s w i t h on e , t w o,a n d t h r e e u n d u l a t i on s a r e m a r k e d i n F i g u r e 8 A. F i g u r e8 B s h ow b r a n ch e d p or t i on s of a cy l i n d r i ca l m i ce l l e w i t hq u a n t i z e d u n d u l a t i on s a p p a r e n t l y d i ct a t i n g t h e a r m l e n g t h s . M u l t i p l e u n d u l a t i n g b r a n ch e s i n F i g u r e 8 Ch i g h l i g h t t h e l oca l i z a t i on of s u ch fe a t u r e s n e a r t h eju n ct i on s a n d e n d s . T h e s e a g g r e g a t e s a r e s t a b l e , a sh e a t i n g t h e s a m p l e t o 5 0 °C for a fe w d a y s d i d n otp r od u ce a n y n ot i ce a b l e ch a n g e i n t h e a s s e m b l e d m or -p h ol og i e s . A com p a r i s on of O B 9 -1 /O B 1 -3 s h ow n i nF i g u r e 4 B a n d O B 9 -1 /O B 1 -5 p r e s e n t e d i n F i g u r e 8s h ow s a t r a n s i t i on fr om m a i n l y s p h e r i ca l m i ce l l e s t om a in ly cylin d r ica l m ice lle s ove r a ve r y n a r r ow com p os i-t i on r a n g e .
B l e n d s a r o u n d w P E O *. R e ce n t l y , w e d i s cov e r e dn e t w or k for m a t i on i n a q u e ou s d i s p e r s i on s of O B 9 -4( w P E O ) N P B ) w P E O > w P E O < 〈 w P E O ) 〉 , i d e n t i ca l t o O B 9 -4 . T o ou rs u r p r i s e t h i s p r e m i x e d b l e n d s e l f-a s s e m b l e s i n t o ap ot p ou r r i of d e l i ca t e l ook i n g ob je ct s w i t h b i l a y e r , cy l i n -d r i ca l , a n d com p l e x ju n ct i on s t r u ct u r a l e l e m e n t s , fr e -q u e n t l y m i x e d w i t h i n i n d i v i d u a l m oi e t i e s . T h e s e u n -u s u a l s t r u ct u r e s a r e e vid e n t in a ll t h e cr yo-T E M im a ge st a k e n fr om t h i s m i x t u r e ; F i g u r e 9 C d i s p l a y s m a n y oft h e p r e va le n t fe a t u r e s . P e r h a p s t h e m os t s t r ik in g is t h eoct op u s (or je l l y fi s h )-l i k e m i ce l l e s , w h i ch a r e com p os e dof a fl a t b i l a y e r w i t h p r ot r u d i n g cy l i n d r i ca l m i ce l l e sa lon g t h e e d ge s . T h e s e oct op u s -lik e e n t it ie s a r e com m on ,a lt h ou gh t h e y occu r w it h a va r ia b le n u m b e r of cylin d r i-ca l a r m s . S e v e r a l e x a m p l e s a r e s h ow n i n F i g u r e 1 0con t a i n i n g 4 , 5 , 6 , 7 , 8 , 9 , 1 0 , a n d 1 4 a r m s a t t a ch e d t ot h e fl a t ce n t r a l p or t i on . I n a l l t h e s e oct op u s -l i k e a g -g r e g a t e s t h e cy l i n d r i ca l a r m s a r e s y m m e t r i ca l l y d i s -t r i b u t e d a l on g t h e ci r cu m fe r e n ce of t h e ce n t r a l fl a tb ila ye r a s is e vid e n t in F igu r e s 9 C a n d 1 0 . O cca s ion a lly,t h e s e ob je ct s a p p e a r t o b e fol d e d on a s i d e w i t h t h ecy l i n d r i ca l a r m s p r ot r u d i n g fr om a h e m i s p h e r i ca l b i -l a y e r ca p . O b v i ou s l y , t h e con fi n e m e n t cr e a t e d b y t h e F i g u r e 8 .
C r y o-T E M im a ge s fr om m ixt u r e O B 9 -1 /O B 1 -5d e p ict i n g u n d u la t i on s a n d d is t e n d e d s p h e r i ca l e n d ca p s onw or m lik e cylin d r ica l m ice lle s . S h or t cylin d e r s w it h a n u n d u la -t ion (s h or t a r r ow ) a n d t w o u n d u la t ion s (lon g a r r ow ) a n dcy li n d e r e n d s w it h on e , t w o, a n d t h r e e b e a d s a r e m a r k e dcor r e s p on d in g ly i n (A). P a n e l s B a n d C s h ow b r a n ch i n g w i t hq u a n t i ze d u n d u l a t i on s fi xi n g t h e s e g m e n t le n g t h b e t w e e nju n ct ion s . T w o t yp e s of h yp e r b olic (s a d d le ) s u r fa ce s ch a r a ct e r -i ze t h e s e m or p h ol ogi e s (s e e F i g u r e 1 3 ).
F i g u r e 9 .
C r y o-T E M m i cr og r a p h s fr om t h r e e d is p e r s i on s w i t h id e n t ica l com p os i t ion s , 〈 w P E O ) 〉 ; t h e m ol e cu l a r w e igh td is t r ib u t ion b r oa d e n s fr om A t o B t o C . (A) A n e t w or k fr a gm e n t fr om O B 9 -4 , a s in g le com p on e n t d i s p e r s ion . T h i s p ict u r e i sr e p r od u ce d fr om r e f 7 . (B ) B le n d O B 9 -1 1 /O B 9 -1 5 . A b im od a l d is t r ib u t ion of com p on e n t com p os it ion s ( w P E O ) w P E O ) M a cr om ol ecu l es , V ol . 3 7 , N o. 4 , 2 0 0 4
P E O - P B M i ce l l a r D i s p e r s i on s
Figure 18: (top-row) Approximate equilibrium states computed from super-critical benchmark initial dataand parameters except for values of η taken as 2 . ε, . ε, . ε, . ε, . ε (left-to-right). Final times are T = 5 K, K, K, K, K respectively. (bottom-row) Experimental comparisons showing (left-to-right)bubbles, bubbles with endcaps, bubbles and branched endcaps, long-branched filaments with endcaps, anddouble-sheeted bubbles (bubble inside of bubble). Figures 3I and 3N from [22], reprinted with Permissionfrom the AAAS. Figures 5A, 5B, and 9C, Reprinted (adapted) with permission from [23]. Copyright (2004)American Chemical Society. Acknowledgement
A. Christlieb acknowledges support from NSF grant DMS-1912183. K. Promislow acknowledgessupport from NSF grant DMS-1813203. Z. Tan recognizes support from the China ScholorshipCouncil under 201906160032. B. Wetton recognizes support from a Canadian NSERC grant. S.Wise recognizes support from NSF grant DMS-1719854.
References [1] S. Barnhill, N. Bell, J. Patterson, D. Olds, and N. Gianneschi, Phase diagrams of polynor-bornene amphiphilic block copolymers in solution,
Macromolecules (2015) 1152-1161.[2] A. Blanazs, S.P. Armes, and A. J. Ryan, Self-assembled block copolymer aggregates: from33icelles to vesicles and their biological applications, Macromolecular rapid communications , (2009) 267–277.[3] L. Chen, Xiaozhe Hu, and S.M. Wise, Convergence Analysis of the fast Sub-space Descent Method for Convex Optimization Problems, Math. Comp. (in press),https://doi.org/10.1090/mcom/3526.[4] Y. Chen, A. Doelman, K. Promislow, F. Veerman,
Robust stability of multicomponent mem-branes: the role of glycolipids , submitted (2019).[5] Y. Chen and K. Promislow, Regularized curve lengthening in the strong FCH gradient flow.arXiv preprint: https://arxiv.org/abs/1907.02196[6] Y. Chen and K. Promislow, Amphiphilic blends, derivation of a density functional theory.Preprint.[7] S.-H. Choi, T. Lodge, and F. Bates, Mechanism of molecular exchange in diblock copolymermicelles: hypersensitivity to core chain length,
Phys. Rev. Let. (2010) 047802.[8] R. Choksi and S. Ren, On the derivation of a density functional theory for microphase sepa-ration of diblock copolymers,
Journal of Statistical Physics
Physica D (1-2) (2005), 100-119.[10] A. Christlieb, N. Kraitzman, K. Promislow, Competition and Complexity in AmphiphilicPolymer Morphology,
Physica D
Proc. Roy. Soc. London, Series A , : 20120505 (20 pages).[12] S. Dai and K. Promislow, Competitive Geometric Evolution of Amphiphilic Interfaces, SIAMMath. Analysis (2015), 347-380.[13] Q. Du, L. Ju, X. Li, and Z. Qiao, Stabilized linear semi-implicit schemes for the nonlocalCahn-Hilliard equation, J. Computational Physics (2018) 39-54.[14] W. Feng, Z. Guan, J.S. Lowengrub, X. Wang, S.M. Wise, and Y. Chen, A uniquely solv-able, energy stable numerical scheme for the functionalized Cahn–Hilliard Equation and itsconvergence analysis,
J. Sci. Computing (2018) 1938-1967.[15] W. Feng, A.J. Salgado, C. Wang, and S.M. Wise, Preconditioned Steepest Descent Methodsfor some Nonlinear Elliptic Equations Involving p-Laplacian Terms, J. Comput. Phys. (2017) 45-67. 3416] G. Fredrickson, The equilibrium theory of inhomogeneous polymers, International Series ofMonographs on Physics 134, Oxford: Clarendon Press, (2006).[17] N. Gavish, G. Hayrapetyan, K. Promislow, L. Yang, Curvature driven flow of bi-layer inter-faces,
Physica D (2011) 675-693 .[18] G. Gompper and J. Goos, Fluctuating interfaces in microemulsion and sponge phases,
Phys.Rev E. (1994) 1325-1335.[19] G. Gompper and M.Schick, Correlation between structural and interfacial properties of am-phillic systems, Phys. Rev. Lett. (1990) 1116-1119.[20] U Gong, J. Zhao, and Q. Wang, Arbitrarily high-order unconditionally energy stableschemes for gradient flow models using the scalar auxiliary variable approach, arXiv preprintarXiv:1907.04254 (2019).[21] E. Hairer, G. Wanner, and S. P. Nørsett, Solving ordinary differential equations, I, Nonstiffproblems. Springer, Berlin, (1993).[22] S. Jain and F.S. Bates, On the origins of morphological complexity in block copolymer surfac-tants, Science , (5618):460–464 (2003).[23] S. Jain and F. Bates, Consequences of nonergodicity in aqueous binary peo-pb micellar dis-persions. Macromolecules , (2004) 1511–1523.[24] N. Kraitzman and K. Promislow, Pearling Bifurcations in the strong Functionalionalized Cahn-Hilliard Free Energy, SIAM Math Analysis (3) 3395-3426 (2018), DOI:10.1137/16M1108406[25] K. Promislow and Q. Wu, Existence of pearled patterns in the planar Functionalized Cahn-Hilliard equation, J. Differential Equations (2015) 3298-3343.[26] J. Shen, J. Xu, and J. Yang, The scalar auxiliary variable (BDF2-SAV) approach for gradientflows,
J. Computational Physics (2018) 407-416[27] J. Shen and J. Xu, Convergence and error analysis for the scalar auxiliary variable (BDF2-SAV)schemes to gradient flows,
SIAM Numerical Analysis (2018) 2895-2912.[28] J. Shen, J. Xu, and J. Yang, A New Class of Efficient and Robust Energy Stable Schemes forGradient Flows, SIAM Rev. (3), 474-506.[29] M. Teubner and R. Strey, Origin of scattering peaks in microemulsions, J. Chem. Phys. bf 87(1987) 3195-3200. 3530] T. Uneyama and M. Doi, Density Functional Theory for Block Copolymer Melts and Blends,
Macromolecules Macromolecules Phys. Rev. E , 066703 (2003).[33] S. Wise, Unconditionally stable finite difference, nonlinear multigrid simulation of the Cahn-Hilliard-Hele-Shaw system of equations, J. Sci. Computing44