Competition and complexity in amphiphilic polymer morphology

Physica D 400 (2019) 132144

Contents lists available at ScienceDirect

Physica D journal homepage: www.elsevier.com/locate/physd

Competition and complexity in amphiphilic polymer morphology ∗

Andrew Christlieb a , Noa Kraitzman b , , Keith Promislow c a b c

Department of Computational Mathematics, Science, and Engineering, Michigan State University, United States Department of Mathematics, University of Utah, United States Department of Mathematics, Michigan State University, United States

highlights • • • •

Analyze self-assembly of codimension one and two morphologies in amphiphilic castings. Derivation of regularized curve lengthening equations that drive morphological competition. Characterize morphological complexity as transient passage through pearling instability. Comparison of FCH flow with experiment, molecular dynamics, and self-consistent field theory.

article

info

Article history: Received 27 June 2018 Received in revised form 13 May 2019 Accepted 30 June 2019 Available online 22 July 2019 Communicated by Y. Nishiura Keywords: Geometric evolution Functionalized Cahn–Hilliard energy Amphiphilic interface Network formation Multiscale analysis Curvature driven flow

a b s t r a c t The strong functionalized Cahn–Hilliard equation models self assembly of amphiphilic polymers in solvent. It supports codimension one and two structures that each admit two classes of bifurcations: pearling, a short-wavelength in-plane modulation of interfacial width, and meandering, a longwavelength instability that induces a transition to curve-lengthening flow. These two potential instabilities afford distinctive routes to changes in codimension and creation of non-codimensional defects such as end caps and Y -junctions. Prior work has characterized the onset of pearling, showing that it couples strongly to the spatially constant, temporally dynamic, bulk value of the chemical potential. We present a multiscale analysis of the competitive evolution of codimension one and two structures of amphiphilic polymers within the H −1 gradient flow of the strong Functionalized Cahn– Hilliard equation. Specifically we show that structures of each codimension transition from a curve lengthening to a curve shortening flow as the chemical potential falls through a corresponding critical value. The differences in these critical values quantify the competition between the morphologies of differing codimension for the amphiphilic polymer mass. We present a bifurcation diagram for the morphological competition and compare our results quantitatively to simulations of the full system and qualitatively to simulations of self-consistent mean field models and laboratory experiments. In particular we propose that the experimentally observed onset of morphological complexity arises from a transient passage through pearling instability while the associated flow is in the curve lengthening regime. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Amphiphilic molecules are increasingly important in synthetic chemistry where they permit molecular level control of the self assembly of materials with desirable ionic and electronic conduction properties, [1]. A molecule is amphiphilic with respect to a solvent if it is comprised of two components, one of which has an energetically favorable or hydrophilic interaction with the solvent and the other with an energetically unfavorable or hydrophobic interaction. There is a growing literature for the construction and characterization of amphiphilic diblock polymers ∗ Corresponding author. E-mail addresses: [email protected] (A. Christlieb), [email protected] (N. Kraitzman), [email protected] (K. Promislow). https://doi.org/10.1016/j.physd.2019.06.010 0167-2789/© 2019 Elsevier B.V. All rights reserved.

comprised of polymer chains formed from adjustable lengths of hydrophilic and hydrophobic polymers covalently bonded together. Amphiphilic molecules are typically characterized by the Flory–Huggins parameters indicating the strength of the hydrophobic/hydrophilic interactions and by the aspect ratio of the two components, [2]. When immersed in solvent, amphiphilic polymers self assemble into a wide variety of structures with diverse morphologies that include codimension one bilayers (hollow vesicles), codimension two filaments (solid rods or cylinders), codimension three micelles (solid spheres), and various defect structures with no well defined codimension such as end-caps, ‘‘Y’’ junctions, and mixed morphologies. The bifurcation structure of these mixtures has been experimentally investigated for a variety of different

2

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

diblock structures. Two seminal experiments, presented in Fig. 1, have partially unfolded this bifurcation structure by studying the impact of varying the solvent blend to modify strength of the amphiphilic interaction, [3], and of varying the aspect ratio of amphiphilic polymer, [4]. Increasing the strength of the amphiphilic interaction and decreasing the aspect ratio of the minority phase produce similar results: a sequence of bifurcations in which codimension-one bilayers yield to codimension-two filaments which yield in turn to codimension-three micelles. In some experiments these codimensional structures are reported to coexist for a range of parameter values while others show regions of ‘‘morphological complexity’’ characterized by an abundance of defects. In this paper we analyze the coexistence, bifurcation, and longtime evolution of well-separated, defect free, codimension one and two structures within the context of the H −1 gradient flows of the strong scaling of the Functionalized Cahn–Hilliard free energy. We use multiscale analysis to derive the curvature driven evolution of codimension one and codimension two structures that are sufficiently far from self intersection. We show that the evolution of morphologies of distinct codimension couples through their exchange of amphiphilic molecules with the bulk. The bulk chemical potential is spatially constant and varies temporally on a slow time scale. Most significantly, codimension one and two structures switch between a regularized curve lengthening and a curve shortening evolution as the bulk chemical potential passes through critical values. This dichotomy presents a mechanism for morphological competition in which structures of one codimension grow at the expense of the other. This suggests that, in the absence of defects, the coexistence of structures with distinct codimensions is not generic however the resolution into structures of homogeneous codimension may require a substantial transient. This finding is supported by experimental results, see [5], which report that transient structures can persist for months. There are two mechanisms for a codimension one or two structure to develop an initial defect: self intersection and pearling. Pearling bifurcations are high-frequency tangential modulations of the width of the codimensional structure. In a companion paper, [6], we characterized the onset of the pearling bifurcation, showing that within the strong scaling of the functionalized Cahn–Hilliard free energy the onset of pearling is independent of the shape of the codimension one or two structure, but couples to the transient value of the bulk chemical potential. Self intersections can be non-local, arising when the initial intersection points are well separated in distance measured along the curve, as in a Fig. 8 intersection. They can also be local, as arises when a surface develops large curvatures, such as when the radius of a sphere tends to zero. In both cases the self intersections arise from the evolution of the underlying manifold that characterizes the structure. Curve-shortening flows render the manifold smaller and drive their curvatures towards constant values. Conversely, curve lengthening flows act like backward heat equations for the interfacial curvature and are well known to be locally ill-posed. We derive a regularized curve lengthening flow that includes a higher-order surface diffusion. The regularized curve lengthening is locally well posed evolution, but causes the underlying curve to grow and ‘‘meander’’ or buckle, and typically leads to finite-time non-local self intersections. In both cases the finite-time singularities can be arrested by the quenching element of the flow which slows the normal velocity as the far-field chemical potential approaches an equilibrium value. Depending upon initial conditions, our results can predict either a relaxation to an equilibrium state or provide the mechanisms for the generation of defect states. Our analysis is particularly relevant to the study of morphologies derived from casting processes in which an initial suspension

of amphiphilic molecules, reflecting a high bulk chemical potential, nucleates out structures of various codimension which initially grow, absorbing the amphiphilic molecules from the bulk suspension and lowering the bulk chemical potential. As the chemical potential falls it may trigger or inhibit pearling bifurcations in one or both codimensions, or trigger transitions from the regularized curve lengthening to curve shortening. We compare our asymptotic results to simulations of the full system, to simulations of self-consistent mean-field density functional models of amphiphilic polymer melts, and to experimental bifurcation diagrams. 1.1. The functionalized Cahn–Hilliard free energy The Functionalized Cahn–Hilliard (FCH) free energy models the free energy of a binary mixture of an amphiphilic molecule and a solvent. It supports stable network morphologies including codimension one bilayers and codimension two filaments as well as pearled morphologies and the defects such as end-caps and Y -junctions, [7–11]. The FCH free energy takes the form

∫

1(

F (u) := Ω

2

ε 2 ∆u − W ′ (u)

)2

− εp

(

ε 2 η1 2

|∇ u|2 + η2 W (u)

)

dx,

(1.1) where W is a smooth double-well potential with local minima at u = b± with b− < b+ . The two wells have unequal depths that are normalized so that W (b− ) = 0 > W (b+ ) and the left well is non-degenerate in the sense that α− := W ′′ (b− ) > 0. The value of α− is a key parameter that controls the rate of exchange of amphiphilic molecules between the bulk and the various morphologies. Here ε ≪ 1 is small parameter corresponding to the ratio of length of the amphiphilic molecule to the domain size, u = b− is associated to a bulk solvent phase, while the quantity u − b− > 0 is proportional to the density of the amphiphilic phase. The first term in the integrand of (1.1) is called the Willmore or the quadratic term, as it denotes the square of a variational derivative of a Cahn–Hilliard type free energy. The quadratic term is positive, and we refer to the class of u ∈ H 2 (Ω ) for which the residual of the quadratic term is small compared to ε as morphologies. The second grouping of terms in the integrand, multiplied by ε p , is called the functionalization terms. The strong and weak scalings of the FCH free energy correspond to the choice p = 1 and p = 2, respectively in (1.1) and represent two natural choices of distinguished limits between the residual of the quadratic term and the typical scaling of the functionalization terms. In the strong scaling of the FCH, the O(ε ) functionalization terms typically dominate the generically O(ε 2 ) residuals from the quadratic terms, in the weak scaling both terms balance at O(ε 2 ). The analysis of this paper focuses on the strong scaling of the FCH free energy for which the bifurcation analysis is more accessible. The functionalization parameters η1 and η2 characterize key properties of the amphiphilic molecules. Specifically η1 > 0 models the strength of the hydrophilic interaction, modeling the propensity of amphiphilic molecules to form monolayers by rewarding increases in interfacial area or curve length with a decrease in free energy. The parameter η2 ∈ R encodes the aspect ratio of the amphiphilic molecule, as discussed in Section 4.3. Equivalently these parameters are analogous to the surface and volume energies typical of models of charged solutes in confined domains, see [12] and particularly equation (67) of [13]. With these parameter choices the minus sign in front of the functionalization terms has great significance — it incorporates the propensity of the amphiphilic surfactant phase to drive the creation of interface. Indeed, experimental tuning of solvent quality identifies molecular level phase separation and self assembly in

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

3

Fig. 1. Morphological bifurcation diagrams for two classes of amphiphilic polymers in solvent. (left) Morphology of Polystyrene (PS)–Polyacrylic acid (PAA) diblocks as function of increasing water content in water–dioxane solvent blend (horizontal axis) and polymer %-weight fraction of the overall mixture (vertical axis). Increased water volume fraction drives bifurcation to lower codimensional morphologies without inducing pearling. From [3]. Reprinted with permission from AAAS. (right) Morphology of amphiphilic Polyethylene oxide (PEO)–Polybutadiene (PB) diblock suspension as function of the PEO weight fraction, wPEO (horizontal axis) for molecular weights of PB fixed at NPB = 45 and 170 (vertical axis). Defect loaded phases are observed for the higher value of NPB . Images of polymer morphologies are included inset, as well as schematic representations of codimension one, two, and three morphologies and ‘Y’ junctions. From [4]. Reprinted with permission from AAAS.

amphiphilic mixtures with the onset of negative values of surface tension in mesoscale agglomerates, [14] and [15]. Prior work on FCH gradient flows has focused on the weak scaling, corresponding to p = 2 in (1.1). In [9] the authors derived the geometric evolution of bilayers at the quenched meancurvature flow on the O(ε −1 ) time scale and as a surface area preserving Willmore flow on the O(ε −2 ) time scale. The geometric evolution of codimension two structures was derived in [16], obtaining a curvature driven competitive geometric evolution of the filament curve on the O(ε −1 ) time scale and a length-preserving Willmore type flow on the O(ε −2 ) time scale. Moreover, it was found that the codimension one and two structures can co-exist on the faster O(ε −1 ) time scale in the weak FCH, but compete on the longer O(ε −2 ) time scale. However, rigorous investigation of the pearling bifurcation in the weak FCH is complicated by its leading order coupling to the curvature of the underlying curves. Conversely in the strong scaling of the FCH, the pearling bifurcation is independent of morphology and was rigorously characterized in [6] for a wide class of codimension one and two structures. In the remainder of this paper we consider the strong scaling of the FCH free energy. Fixing Ω = [0, L]d ⊂ Rd for d = 2, 3, . . . and applying periodic boundary conditions to H 4 (Ω ), the first variation of F , also called the chemical potential µ, associated to a spatial distribution u ∈ H 4 (Ω ) takes the form

driven flow mediated through the common value of the spatially constant far-field chemical potential, µ1 defined in (2.17). We show in Section 2 that the nonlocal Mullins–Sekerka problem familiar to Cahn–Hilliard evolution is present but is unforced, and on the long time scales we consider the far-field chemical potential relaxes to a trivial, spatially-constant for both codimension one bilayer and codimension two filament morphologies. As a consequence, the geometric evolution is local. While spatially constant, the far-field chemical potential µ1 = µ1 (t), is temporally dynamic and is linearly proportional to the density of free amphiphilic molecules in the bulk. It serves as a key bifurcation parameter, triggering two potential types of instability for each codimension of morphology. Indeed, in [17] it is shown for the FCH free energy that the pearling and self intersection via geometric motion are the only possible mechanisms to generate defects in bilayers. In the companion paper, [6], a sharp condition for pearling stability is derived that relates the chemical potential to the parameter ηd and constants that depend implicitly on the form of the double-well W . Specifically the bilayers are stable with respect to the pearling bifurcation if and only if

δF µ := (u) = (ε 2 ∆ − W ′′ (u) + εη1 )(ε 2 ∆u − W ′ (u)) + εηd W ′ (u), δu

( 2  2 ) µ1 Sf + ηd ψf′ ,0,0 L + λf ,0,0 ψf ,0,0 L < 0,

(1.2) where ηd := η1 − η2 . The Functionalized Cahn–Hilliard equation is the associated H −1 gradient flow, ut = ∆µ(u),

(1.3)

supplemented with periodic boundary conditions on Ω . The choice of the H −1 gradient is a reflection of its status as the simplest local gradient that preserves mass. The mathematical focus of the paper is on the multiscale analysis of the evolution codimension one and two structures. On the O(ε −1 ) time scale we find that the H −1 gradient flow drives well separated filament and bilayer structures through a competitive, mean-curvature

 2 µ1 Sb + ηd λb,0 ψb,0 2 < 0,

(1.4)

and similarly filaments are pearling stable if and only if R

R

(1.5)

where λb,0 is the ground-state eigenvalue of the linear operator Lb,0 , defined in (2.7), with eigenfunction ψb,0 , and λf ,0,0 is the ground state eigenvalue of the linear operator Lf ,0,0 , defined in (3.16), with the corresponding eigenfunction ψf ,0,0 . The constants Sb , Sf are the shape factors of the bilayers and the filaments, respectively, defined by the relations

∫

Φb,1 W ′′′ (φb )ψb2,0 dz , ∫ ∞ Sf := 2π Φf ,1 W ′′′ (φf )ψf2,0,0 RdR. Sb :=

R

(1.6)

0

Here φb and φf are the bilayer and filament profiles, defined in (2.6) and (3.10), while Φb,1 and Φf ,1 , defined in (2.8) and

4

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

(3.17), encode the impact of a change in chemical potential on the shape of the bilayer and filament, respectively. For each codimension, if the shape factor is negative then pearling stability is favored by large (positive) values of µ1 , while if it is positive then pearling stability is favored for small (negative) values of µ1 . In [11] the existence of pearled codimension one circular and flat equilibrium was demonstrated in R2 for the strong FCH. 1.2. Summary of analytical results Our main analytical result is the derivation of curvature driven flow laws valid on slow time τ = ε t for codimension one and two structures embedded in Ω ⊂ R3 via a multiscale analysis. Bilayer and filament morphologies are defined as dressings of collections of admissible interfaces and curves, respectively. We consider an admissible codimension one interface Γb , see Definition 2.1, with O(1) area and codimension two curve Γf , see Definition 3.1, with O(ε −1 ) length. We assume the reaches of these two classes of surfaces, as defined in (2.2) and (3.4) respectively, are all disjoint and introduce the composite solution defined in (4.2). The evolution of these composite solutions under the gradient flow (1.3) is parameterized at leading order by the triple (Γb , Γf , µ1 ) according to their ε -scaled normal velocities Vb = νb (µ1 − µ∗b )H0 + ε kb ∆s H0 ,

(1.7)

Vf = − νf (µ1 − µf )κ⃗ + ε kf ∂ κ⃗ .

(1.8)

2 s

∗

[

]

Here we have introduced the mean curvature H0 and Laplace– Beltrami operator ∆s of Γb and vector curvature κ⃗ = (κ1 , κ2 ) and surface diffusion ∂s2 of Γf . On this slow time τ = ε t, the chemical potential satisfies

[ ∫ 2 ( ) α− =− mb νb µ1 − µ∗b H02 − ε kb |∇s H0 |2 ds dτ |Ω | Γb ] ∫ ( ) 2 2 ∗ κ | − ε|∂s κ⃗f | ds + 2 π mf ε νf µ1 − µf |⃗

dµ1

Γf

+ O(ε ), 2

(1.9)

where the constants νb > 0, kb > 0, and mb > 0, are defined in (2.51) while νf > 0, kf > 0, and mf > 0 defined in (3.77). While the surface diffusion terms are formally lower order, they are leading order singular perturbations that keep the resulting flow locally well-posed. The system (1.7)–(1.9) describes the competitive dynamics between codimension one and two morphologies. The key critical values

µ∗b = − µ∗f =

kb 2νb

kf

νf

(η1 + η2 ),

(1.10)

η1 ,

(1.11)

indicate the value of µ1 at which the rates of amphiphilic molecule insertion and ejection are balanced for bilayers and filaments, respectively. Specifically, if the far-field chemical potential lies above this number, then structures of the corresponding codimension will grow. Indeed, the rate of change of the area of bilayers is given by

∂τ |Γb | =

∫ Γb

∫ Vb H0 ds =

Γb

νb (µ1 − µ∗b )H02 − εkb |∇s H0 |2 ds, (1.12)

with the corresponding expression

∂τ |Γf | =

∫ Γf

Vf · κ⃗ ds =

∫ Γf

νb (µ1 − µ∗f )|⃗ κ |2 − ε kf |∂s κ⃗ |2 ds, (1.13)

for the length of the filaments. The competitive dynamics system provides the leading order evolution so long as the interfaces remain admissible with disjoint reaches and the µ1 -dependent pearling conditions (1.4)–(1.5) hold. The surface diffusion terms in (1.7)–(1.8) are relevant to mass balance only if the gradients of the curvatures become asymptotically large or if µ1 becomes asymptotically close to one of the critical values µ∗b or µ∗f . In particular, it follows from (1.9) that net growth of bilayers and filaments corresponds to a decrease in µ1 , while large curvature gradients enhance the ejection rates and increase the value of µ1 . Remark. For the codimension two filament term to contribute to the evolution of the chemical potential, µ1 , at leading order, we assume that their collective length |Γf | is O(ε −1 ). Our generic assumption on the codimension one phase is the surface area |Γb | = O(1), so that both bilayers and filaments occupy an O(ε ) volume fraction. This limits the applicability of the asymptotic results as the assumption of disjoint reaches becomes non-generic, and represents a significant caveat in the application of our analytical results. Since the geometric flow reduction does not apply to morphologies with defects, both a large number of short filaments or a small number of long filaments complicate the non-self intersection assumption. Our analysis applies to each disjoint component of filament morphology separately, and while ε ≪ 1 is small, it is viewed as fixed within the model and need not be vanishingly small. A more detailed analysis of mass scaling in the ε → 0 convergence issues for FCH models can be found in [19]. For µ1 > µ∗b we call the normal velocity (1.7) a regularized (codimension one) curve lengthening flow and a (codimension one) curve shortening flow if µ1 < µ∗b , with similar terminology for the codimension two flow based upon the sign of µ1 − µ∗f . When the structures have a homogeneous codimension, then in the absence of singularities in the curvature flow, Eq. (1.9) drives the chemical potential µ1 to the corresponding critical value, µ∗b or µ∗f , and the leading-order term in the geometric flow goes to zero, and the system is said to be ‘‘quenched’’. To illustrate the nature of the geometric flow, it is instructive to rewrite it as a corresponding evolution equation for the curvatures. For codimension one structures in two space dimensions, up to tangential reparameterization it takes the simple form

( ) ∂τ H0 = −(∆s + H02 )Vb = −(∆s + H02 ) νb (µ1 − µ∗b )H0 + ε kb ∆s H0 . (1.14) For the curve lengthening regime, the dominant term is a backward heat equation, with a fourth-order regularization and a H03 nonlinearity with a negative (stable) coefficient. For the curve shortening regime, both differential terms are stabilizing, but the cubic nonlinearity has a positive coefficient that supports finite time singularity which may be arrested by the fourth order regularization. In numerical simulations the curve lengthening flows, with µ1 > µ∗b,f show distinct regimes. For positive but O(ε ) values of µ1 − µ∗b the bilayer interfaces will bend and buckle at O(1) length scales, leading to shapes reminiscent of a meandering river. This regime is called a ‘‘meandering flow’’, and is studied rigorously in [20]. For O(1) positive values of µ1 − µ∗b the curve lengthening flows can lead to growth of high-curvature regions which self intersect on a τ = O(1) length scale. For filaments this can lead to the formation of many closed loops, see Figures 2 and 3 of [4] for experimental examples or Fig. 5 of this article for an example of meandering motion within the FCH gradient flow.

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

5

Fig. 2. Molecular dynamics simulations of amphiphilic molecules with distinct aspect ratios. Increasing the size of the hydrophilic head group, relative to a fixed tail leads to a preference for morphologies with increasing codimensions: (b) bilayers with edge caps (codim 1), (c) branched filaments (codim 2), (d) unbranched filaments (codim 2), and (e) micelles (codim 3). Morphologies with higher codimension have a higher density of free-floating amphiphilic molecules in the far-field (bulk), corresponding to a higher critical value of µ1 . Source: Reprinted from [18], with permission from Elsevier.

1.3. Competition and morphological complexity Our main scientific results are the conjectured mechanisms for the morphological bifurcations observed in the casting of amphiphilic suspensions. In particular we propose a mechanism for the onset of the so-called ‘‘morphological complexity’’ observed in the experimental casting processes presented in Fig. 1 (right). In a casting process amphiphilic molecules are dispersed (mixed) in a solvent and the mixture is allowed to relax, generically leading to self-assembly of structures with distinct codimension. For the shorter chains, NPB = 45 corresponding to the lower horizontal row of symbols, casts from molecules with increasing weight fraction of the amphiphilic PEO component (wPEO – horizontal axis) lead to a series of codimensional bifurcations in which the self-assembly prefers structures with increasing codimension. Indeed experiments produce first bilayers (codimension one) marked (B), then bilayers coexisting with cylinders (codimension two) marked (B + C), then cylinders, then cylinders coexisting with spheres (codimension three) marked (S), and finally spheres. A similar series of codimensional bifurcations are presented in Figs. 1(left) and 2. In the former the bifurcations are induced by lowering the percentage of water in the water– dioxane solvent which forms the basis for the casting. In the latter they are realized within coarse-grained molecular dynamics simulations of casting processes by increasing the size of the hydrophilic head group, and hence the aspect ratio, of a simple three-group amphiphilic molecule. We conjecture that the competitive dynamics implicit in the system (1.7)–(1.9) forms the basis for these codimensional bifurcations. In a casting process there is initially a relatively large density of dispersed amphiphilic molecules, corresponding to a high value of the far-field chemical potential µ1 . As various structures self-organize, amphiphilic molecules are removed from the far-field and the scalar value of µ1 falls. The critical values µ∗b and µ∗f , given in (1.10)–(1.11), gauge the relative ability of the corresponding bilayer and filament morphologies to absorb and retain amphiphilic molecules from the far-field (bulk) environment. The morphology with the lowest corresponding critical value will, in the absence of defects, lower the value of µ1 and drain the mass of the morphologies with higher critical values of µ1 . Increasing values of either η1 or η2 will drive µb greater than µf and trigger a competitive imbalance that favors codimension two filaments over codimension one bilayers. We argue in Section 4 that increasing the aspect ratio of the amphiphilic molecule corresponds to an increase in the value of η2 , while decreasing the percentage of water within the water–dioxane solvent blend

corresponds to an increase in the value of η1 . These produce shifts in µ∗b,f in agreement with bifurcation from codimension one to codimension two. The coexistence of codimension one and two structures for large parameter ranges are not supported by the analysis. However the time scale to reach equilibrium can be quite long, [5] suggest times on the order of months, and we propose that longer experimental trials may decrease the size of the regions of coexistence. We do not present an analysis of codimension three micelles within this work as they do not have a spatially extended direction that can accommodate incremental growth, rather simulations suggest that they swell, form dumbbell shapes, and then break into distinct micelles. This behavior is outside the scope of our analysis. For the longer NPB = 170 chains in Fig. 1(right), increasing wPEO one finds that the codimensional bifurcation structure is interrupted by the onset of so-called ‘‘morphological complexity’’. Specifically, the casting sequences yield bilayers, bilayers coexisting with branched filaments, strongly connected network morphologies, and Y -junction dominated filaments, before reverting to the familiar codimension two and codimension three structures. The term morphological complexity refers not only to the wide variety of possible outcomes, but also to the difficulty in controlling the outcomes, see [4] and [5]. Our second conjecture is that morphological complexity arises from the interplay between the pearling bifurcation, the competitive dynamics, and the evolving value of µ1 . Indeed, the criteria for pearling stability depends upon the value of µ1 , see (1.4) and (1.5), and as µ1 deceases during the casting it may trigger or inhibit pearling stability. In particular, in Section 4.4.2 we present regimes in which bilayers have a competitive advantage over filaments, but are transiently pearling unstable, while filaments are globally pearling stable. This cascade of bifurcations provides mechanisms to produce complex blends of defects and morphologies and affords a clear mechanism for hysteresis. In such an environment the ultimate outcome of a given casting process could depend sensitively on secondary effects such as the rate at which amphiphilic molecules are initially added to the dispersion or upon spatial inhomogeneities. The spatial complexity of the end states in this regime is born out both by experiments and by simulations of the FCH free energy, see Fig. 10 (center and right). We emphasize that the morphological complexity conjecture encompasses structures with codimensional defects that are outside of our analysis. Moreover, the pearling bifurcation cannot be robustly suppressed within the FCH gradient flows. The experiments and simulations exhibiting the simple codimensional bifurcation route do not display signs of pearling bifurcation.

6

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

In particular there is no mechanism within the FCH energy to explain why pearling would be expressed in longer polymers and inhibited in shorter but otherwise identical polymers. These limitations of the model are expanded upon in the discussion of Section 5. In [21] and [22] two-component extensions of the FCH are proposed which possess more detailed internal layer structure and afford precise mechanisms to robustly inhibit pearling bifurcations. In Sections 2 and 3, respectively, we apply a multiscale analysis to derive the long time-scale evolution of admissible codimension one and two structures under the FCH equation. In particular we extract the coupling of the bulk chemical potential on the curvature driven flow. In Section 4 we extend these results to composite morphologies, deriving the laws governing their competitive evolution. We present bifurcation diagrams that show the regions of pearling stability and curve shortening of each morphology, and compare them to simulations of the FCH equation, to self consistent mean field density functional theory simulations, and to the experimental bifurcation diagrams presented in Fig. 1. These are the first results that explicitly quantify the complexity of the transients associated to amphiphilic polymer blends and identify the role of the pearling bifurcation in the generation of complex network morphologies.

For an admissible codimension one interface Γb the change of variables x → ρ (s, z) given by (2.1) is a C 4 diffeomorphism on the reach of Γb , defined as the set

⏐ { } ⏐ Γb,ℓ := ρ (s, z) ∈ Rd ⏐s ∈ S , −ℓ/ε ≤ z ≤ ℓ/ε ⊂ Ω ,

(2.2)

with complement Γ˜ b,ℓ := Ω \Γb,ℓ . On the reach we may expand the Cartesian Laplacian in terms of the Laplace–Beltrami operator ∆s and the curvatures,

ε 2 ∆x = ∂z2 + ε H0 (s)∂z + ε 2 (zH1 ∂z + ∆s ) + O(ε 3 ),

(2.3)

where Hi (s) is related to the ith powers of the curvatures Hi = (−1)i

d−1 ∑

kij+1 ,

(2.4)

j=1

and, in particular, H0 is the mean curvature of Γb . See [9] for more details.

2. Geometric evolution of codimension one structures In this section we derive the geometric evolution of admissible codimension one interfaces, which we refer to as bilayers. These calculations are carried out in Rd for d ≥ 2 but we will restrict our attention to R3 for the analysis of filaments. We consider the local, mass-preserving H −1 gradient flow of the strong FCH given in (1.3). The multiscale analysis of this section follows closely from the calculations of [16], which considered the weak scaling of the FCH gradient flow. For brevity we present only the main calculations. It is well known that for the single-layer interfaces supported in Cahn–Hilliard type models, that is codimension one interfaces which separate distinct phases, the O(1) and O(ε −1 ) time scales yield Stefan and Mullins–Sekerka problems for the interfacial motion [23–25]. For single layers motion of the interface requires transport of materials on either side. For bilayers we derive these reduced flows, but they have trivial solutions, as the interfacial motion of an interface with the same material on either side does not require long-range transport but is facilitated by permeation. The result is a local geometric flow, driven by membrane curvatures and coupled to the bulk value of the chemical potential. 2.1. Admissible codimension one manifold and their dressings Given a smooth, closed (d − 1)–dimensional manifold Γb immersed in Ω ⊂ Rd , we define the local ‘‘whiskered’’ coordinates system in a neighborhood of Γb via the mapping x = ρ (s, z) := ζb (s) + εν (s)z ,

(i) The W 4,∞ (S ) norm of the 2nd Fundamental form of Γb and its principal curvatures are bounded by K . (ii) The whiskers of length 3ℓ < 1/K , in the unscaled distance, defined for each s0 ∈ S by, ws0 := {x : s(x) = s0 , |z(x)| < 3ℓ/ε}, neither intersect each-other nor ∂ Ω (except when considering periodic boundary conditions). (iii) The surface area, |Γb |, of Γb is bounded by K .

(2.1)

Definition 2.2. Given an admissible codimension one interface Γb ∈ GKb ,ℓ and f : R → R which tends to constant value f∞ at an exponential rate as z → ±∞, then we define the H 2 (Ω ) function fΓb (x) := f (z(x))χ (|z(x)|/ℓ) + f∞ (1 − χ (|z(x)|/ℓ)),

where χ : R → R is a fixed, smooth cut-off function which takes values one on [0, 1] and 0 on [2, ∞). We call fΓb ∈ L2 (Ω ) the dressing of Γb with f ∈ L2 (R), and by abuse of notation will drop the Γb subscript when doing so creates no confusion. Within the reach Γb,ℓ of an admissible Γb the quadratic term within the FCH, (1.1), can be re-written in the codimension one whiskered coordinates system (2.1). Setting the quadratic term equal to zero, and formally taking the leading order terms in ε leads to a second-order ODE in z. The bilayer profile φb , is defined to be the solution of this equation

∂z2 φb = W ′ (φb ),

(2.6)

which is homoclinic to the left well b− of W . We denote by Ub ∈ L2 (Ω ) the dressing of Γb by φb ∈ L2 (R), and introduce the associated linear operator Lb,0 := ∂z2 − W ′′ (φb ).

(2.7)

This is a Sturm–Liouville operator on L2 (R) and has a positive ground-state eigenvalue λb,0 > 0 with eigenfunction ψb,0 ≥ 0 and a translational eigenvalue λb,1 = 0 associated to the eigenfunction ψb,1 = φb′ . In addition, we define the functions

where ζb : S ↦ → R is a local parameterization of Γb and ν (s) is the outward unit normal to Γb . The variable z is often called the ε -scaled, signed distance to Γb , while the variables s = (s1 , . . . , sd−1 ) parameterize the tangential directions of Γb .

Φb,j := Lb,0 1,

Definition 2.1. For any K , ℓ > 0 the family, GKb ,ℓ , of admissible interfaces is comprised of closed (compact and without boundary), oriented d − 1 dimensional manifolds Γb embedded in Rd , which are far from self-intersection and with a smooth second fundamental form. More precisely,

2.2. Inner and outer expansions

d

(2.5)

−j

(2.8)

for j = 1, 2 which converge to a non-zero value at z = ±∞. We also define their Γb dressings, which are denoted by Φb,1 and Φb,2 by the abuse of notation mentioned above.

Assuming initial data arising from the dressing of an admissible initial codimension one interface Γb (t0 ) ∈ GKb ,ℓ , we describe the coupled geometric evolution of the interface and the far-field

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

chemical potential as a flow in time t. We consider formal, multiscale analysis of the density u and the chemical potential µ. In the far-field or bulk region, Γ˜ b,ℓ := Ω \Γb,ℓ , the outer solution u has the expansion u(x, t) = u0 (x, t) + ε u1 (x, t) + O(ε 2 ).

(2.9)

Within the reach or inner region, Γb,ℓ , we use the whiskered coordinates with the ε -scaled distance z to the interface. The standard assumption is that the inner solution u˜ is smooth in the tangential s-variables. However the leading order result is a curvature driven flow whose coefficient may switch sign. Flow against curvature is not locally well posed due to uncontrollable growth in high-frequency growth terms. To regularize this, we incorporate a two-scale tangential expansion, introducing the fast √ tangential variable S := s/ ε , so that the inner variable admits the expansion. u(x, t) = u˜ (s, S , z , τ ) = u˜ 0 (s, S , z , τ ) + ε u˜ 1 (s, S , z , τ ) + O(ε 2 ). (2.10)

7

the Laplacian and the inner solution to obtain an expansion for P = P0 + ε P1 + ε 2 P2 + · · ·, where P0 =∂z2 − W ′′ (u˜ 0 ),

(2.20)

P1 =H0 (s, S)∂z + ∆S − W (u˜ 0 )u˜ 1 + η1 , ′′′

(2.21)

P2 =zH1 (s, S)∂z + D2,S + ∆s − W ′′′ (u˜ 0 )u˜ 2 +

1 2

W (4) (u˜ 0 )u˜ 21 ,

and for Q = Q0 + ε Q1 + ε 2 Q2 + · · · Q0 =∂z2 u˜ 0 − W ′ (u˜ 0 ),

(2.23)

Q1 =H0 ∂z u˜ 0 + ∆S u˜ 0 + (∂z2 − W ′′ (u˜ 0 ))u˜ 1 , Q2 =zH1 ∂z u˜ 0 + D2,S u˜ 0 + ∆s u˜ 0 + (H0 ∂z + ∆S )u˜ 1

(2.24)

1 (2.25) 2 With these reductions we expand the inner chemical potential as

+ (∂z2 − W ′′ (u˜ 0 ))u˜ 2 − W ′′′ (u˜ 0 )u˜ 21 .

µ ˜ 0 =P0 Q0

(2.26)

µ ˜ 1 =P1 Q0 + P0 Q1 + ηd W (u˜ 0 ), ′

(2.27)

µ ˜ 2 =P0 Q2 + P1 Q1 + P2 Q0 + ηd W (u˜ 0 )u˜ 1 . ′′

The inclusion of the fast tangential variable promotes a formally lower order surface diffusion term to the leading order, where it regularizes the curve lengthening flow. The normal velocity V = V (s, S , t) of Γb is denoted by V (s, S , t) := −ε

∂z . ∂t

The second order form of the H gradient induces inner–outer matching conditions for the chemical potential,

µ± ˜ 0 (z , s, S , t), 0 (x, t) = lim µ z →±∞

(2.11)

∂ u˜ ∂ u˜ ∂ u˜ ∂τ = −ε −1 V (s) + . ∂t ∂z ∂τ ∂ t

(2.12)

At the interface we have the standard matching conditions lim u(x + hν, τ ) = lim u˜ (s, S , z , τ ).

(2.13)

z →±∞

which reduce to the relations

˜ u± 0 (x, τ ) = lim u0 (s, S , z , τ ),

(2.14)

± ˜ u± 1 (x, τ ) + z ∂ν u0 (x, τ ) = lim u1 (s, S , z , τ ).

(2.15)

z →±∞ z →±∞

where ∂ν is the derivative in the normal direction of Γb , and u± i denote the values of the limits of the left-hand side of (2.13) as h → 0± respectively. The chemical potential, µ, defined in (1.2), admits similar outer and inner expansions. The terms of its outer expansion are slaved to the density u through the outer relations,

µ0 =W ′′ (u0 )W ′ (u0 ),

(2.16)

µ1 =(W (u0 )u1 − η1 )W (u0 ) + (W (u0 )) u1 + ηd W (u0 ). ′′′

′

′′

2

′

(2.17)

For the inner expansion we introduce the nonlinear operators P and Q to rewrite the chemical potential (1.2) as

µ ˜ = P(u˜ )Q (u˜ ) + εηd W ′ (u˜ ).

(2.18)

In the multiscale tangential variables the Laplacian expansion (2.3) takes the form

ε ∆x = ∂ + ε (H0 (s, S)∂z + ∆S ) 2

(2.28)

−1

(2.29) ± µ± ˜ 1 (s, S , z , t) 1 (x, t) + z ∂ν µ0 (x, t) = lim µ

On the slow time τ = ε t, the time derivative of the inner density function u˜ , defined in (2.10), combined with the normal velocity, (2.11), takes the form

h→0±

(2.22)

2 z

+ ε2 (zH1 (s, S)∂z + zD2,S + ∆s ) + O(ε3 ),

(2.19)

where ∆S is the scaled Laplace–Beltrami operator and D2,S denotes a higher order elliptic term in S. Details on the D2,S term can be found in section 6 of [17] however its precise form is immaterial to our presentation. We combine the expansion of

z →±∞

(2.30) 1

± µ± z 2 ∂ν2 µ± ˜ 2 (s, S , z , t). 2 (x, t) + z ∂ν µ1 (x, t) + 0 (x, t) = lim µ

2

z →±∞

(2.31) 2.3. Time scale τ = ε t: quenched curvature driven flow We focus on the first relevant slow time-scale, τ = ε t, inserting the time derivative and chemical potential expansions into the FCH gradient flow, (1.3). As the interface Γb is codimension one, it separates the region Ω into two disjoint sets, Ω \Γb = Ω+ ∪ Ω− with the normal to Γb pointing towards Ω+ . In this bulk region we obtain the relations O(1) : 0 = ∆ W ′′ (u0 )W ′ (u0 ) ,

(

in Ω− ∪ Ω+ , (2.32)

)

O(ε ) : u0,τ = ∆ (W (u0 )u1 − η1 )W (u0 )

(

′′′

) +(W ′′ (u0 ))2 u1 + ηd W ′ (u0 ) ,

′

in Ω− ∪ Ω+ . (2.33)

The relevant solution to the O(1) relation is the spatially constant density, u0 = b− for which µ0 = 0. With this reduction and the fact that W ′ (b− ) = 0 and W ′′ (b− ) := α− > 0, the O(ε ) relation reduces to 0 = ∆ u1

in Ω− ∪ Ω+ ,

(2.34)

which is subject to interior layer matching and exterior boundary conditions derived in the sequel. In the inner region we supplement the expansion (2.10) with the form of the Laplacian in inner variables, given in (2.3). Collecting orders of ε we find O(ε −2 ) :

0 = ∂z2 µ ˜ 0,

in Γb,ℓ ,

(2.35)

where µ ˜ 0 is defined in (2.26). This relation is satisfied by u˜ 0 = Ub , which is consistent with our choice of initial data corresponding to the dressing of an admissible bilayer with the bilayer profile. Modulo this form the next orders in the expansion take the form, O(ε −1 ) : 0 = ∂z2 µ ˜ 1, O(1) :

− V (s, S)∂z u˜ 0 = ∂z2 µ ˜ 2 + (H0 ∂z + ∆S )µ ˜ 1,

in Γb,ℓ ,

(2.36)

in Γb,ℓ ,

(2.37)

8

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

where µ ˜ 1 and µ ˜ 2 are defined in (2.27) and (2.28) respectively. With the reduction u˜ 0 = Ub , the matching condition (2.30) reduces to the relations ∂z µ ˜ 1 ≈ ∂ν µ0 = 0 as z −→ ±∞. Applying these to (2.36) we find that µ ˜ 1 is independent of z, i.e., µ ˜1 = µ ˜ 1 (s, S , τ ). Similarly, we simplify the inner expression for µ ˜ 1 in (2.27) which yields the expression 2 ′ u˜ 1 = µ ˜ 1 Φb,2 − ηd L− b,0 W (Ub ),

(2.38)

where Φb,2 is defined in (2.8). Since µ ˜ 1 is independent of z, the relation (2.37) reduces to O(1) :

˜ 2 + ∆S µ ˜ 1, − V (s, S)∂z Ub = ∂z2 µ

in Γb,ℓ .

(2.39)

To obtain interfacial jump conditions for µ1 we introduce Uˆ b := Ub − b− > 0, and integrate (2.39) twice from 0 to z we obtain the relation

µ ˜ 2 (z) = µ ˜ 2 (0) − V (s, S)

z

∫

(

Uˆ b (t) dt + z ∂z µ ˜ 2 (0) + V (s, S)Uˆ b (0)

)

0

z

2

∆S µ ˜ 1 (s, S). (2.40) 2 Comparing to the jump condition (2.31), and recalling that µ0 = 0, we deduce that ∆s µ ˜ 1 = 0. Since Γb is closed, this implies that µ ˜1 = µ ˜ 1 (s, τ ) is constant in S. Using this information, we integrate (2.39) with respect to z over R. As Ub is homoclinic we obtain the relationship

+

lim ∂z µ ˜ 2 (z) − lim ∂z µ ˜ 2 (z) = 0,

z →∞

(2.41)

z →−∞

which when reported to the matching condition (2.31) yields the key outer interfacial jump relations Jµ2 K = 0,

J∂ν µ1 K = 0.

(2.42)

Coupling these boundary conditions with the elliptic problem (2.34) typically yields a Mullins–Sekerka problem for the longrange transport of material; however the homogeneous jump conditions imply that ∆µ1 = 0 in all of Ω , which subject to the exterior boundary conditions implies that the far-field chemical potential, µ1 , is spatially constant; however it remains a function of time. To extract the normal velocity we return the reduction u˜ 0 = Ub to (2.28), so that P0 reduces to Lb,0 and the chemical potential takes the form

µ ˜ 2 =L2b,0 u˜ 2 − Lb,0 Q˜ 2 + (H0 (s, S)∂z + ∆S − W ′′′ (Ub )u˜ 1 + η1 )(Lb,0 u˜ 1 + H0 (s, S)Ub′ ) + ηd W ′′ (Ub )u˜ 1 ,

V (s, S)Uˆ b (0) = lim ∂z µ ˜ 2 − ∂z µ ˜ 2 (0) = ∂ν µ1 − ∂z µ ˜ 2 (0) z →−∞

(2.44)

Using (2.44) to replace V (s)Uˆ b (0) in Eq. (2.40) yields

µ ˜ 2 (z) = µ ˜ 2 (0) − V (s, S)

z

∫

Uˆ b (t) dt .

(2.45)

L2b,0 u˜ 2 =Lb,0 Q˜ 1 − H0 ∂z − W ′′′ (Ub )u˜ 1 + η1

)(

− ∆S H0 (s, S)∂z Ub − ηd W ′′ (Ub )u˜ 1 ∫ z +µ ˜ 2 (0) − V (s, S) Uˆ b (t) dt . 0

(

H0 (Lb,0 u˜ 1 , ∂z2 φb )L2 (R) + (W ′′′ (φb )u˜ 1 , ∂z φb ) − η1 ∥φb′ ∥2L2

)

+ V ∥φˆ b ∥2L2 − ∆S H0 ∥φb′ ∥2L2 = 0.

(2.47)

From (2.6) it is easy to verify that Lb,0

(z 2

) φb′ = φb′′ ,

(2.48)

Lb,0 φb′′ = W ′′′ (φb )(φb′ )2 .

(2.49)

Using these relations and the form (2.38) of u˜ 1 , the coefficient of H0 in (2.47) reduces to (Lb,0 u˜ 1 , ∂z2 φb )L2 (R) + (W ′′′ (φb )u˜ 1 , ∂z φb ) − η1 ∥φb′ ∥2L2 z = (L2b,0 u˜ 1 , φb′ )L2 − η1 ∥φb′ ∥2L2 , 2 1 = µ1 mb + (η1 + η2 ) σb . 2 Returning this reduction to (2.47) and solving for the normal velocity we find

µ1 mb + 21 (η1 + η2 )σb

V (s, S) =

B1

H0 +

σb B1

∆S H0 ,

(2.50)

where here and above we have introduced the positive constants

∫

 2  

B1 := φˆ b 

φˆ b dz > 0,

mb := R

L2 (R)

 2 σb := φb′ L2 (R) . (2.51)

,

The sign of the coefficient of H0 is indeterminate, as η2 can be negative and moreover the bulk chemical potential µ1 varies temporally. To emphasize this fact we introduce the constants 1

σb

2

mb

µ∗b = − (η1 + η2 )

,

νb :=

mb B1

,

kb :=

σb B1

,

(2.52) and return the S variable to its original scaling, obtaining the regularized curvature driven flow V (s) = νb µ1 − µ∗b H0 + ε kb ∆s H0 .

(

)

(2.53)

To close the system and fully determine the normal velocity we evoke conservation of mass to specify the temporally varying value of the bulk external chemical potential, µ1 . The mass balance is determined by the interplay between the length of the interface Γb and the total mass of amphiphilic material. From the form of u˜ 0 and u˜ 1 in (2.38), we have the composite formulation 2 ′ 2 u(x, t) = Ub + ε (µ1 Φb,2 − ηd L− b,0 W (Ub )) + O(ε ),

(2.54)

which has the spatially constant far-field asymptotic value

0

Replacing µ ˜ 2 in (2.45) with its expression from (2.43) and solving for L2b,0 u˜ 2 yields an expression for u˜ 2

(

are either functions of z or of s and S, we factor out the functions of s and S, replace Ub with φb , and take the inner product of (2.46) with ∂z φb in L2 (R). Recalling that u˜ 1 , defined in (2.38), is even in z and the operator Lb,0 preserves symmetry, parity considerations reduce the solvability condition to

(2.43)

where we have introduced the quantity Q˜ 2 := Q2 − Lb,0 u˜ 2 . Integrating (2.39) from z = −∞ to z = 0, using the matching condition (2.31), and recalling that µ1 is constant, ∆S µ ˜ 1 = 0, and Uˆ b → 0 as z → ±∞ yields the expression

= −∂z µ2 (0).

Fixing the values of s and S, this equation has a solution u˜ 2 (s, S , ·) ∈ L2 (R) if and only if the right-hand side is perpendicular to ker Lb,0 , which is spanned by ∂z φb . This solvability condition is enforced by selecting the value of V (s, S). Since the terms in (2.46)

Lb,0 u˜ 1 + H0 ∂z Ub

u(x, t) = b− + ε

µ1 + O(ε 2 ) in Ω \Γb,ℓ , 2 α−

(2.55)

)

where α− := W ′′ (b− ) > 0. The gradient flow (1.3) conserves the total mass,

(2.46)

M :=

∫ Ω

u(x, t) − b− dx =

∫ Ω

u(x, 0) − b− dx.

(2.56)

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

Using the form of the composite solution, (2.54), we evaluate the integral over the reach and its complement, M=ε

∫ Γ˜ b,ℓ

µ1 dx + 2 α−

∫ Γb,ℓ

2 ′ 2 Uˆ b +ε (µ1 Φb,2 −ηd L− b,0 W (Ub )) dx + O(ε ).

(2.57) Since Γb is admissible, its area |Γb | ∼ O(1). Changing to whiskered coordinates in the localized integral yields M=ε

(

|Ω |

µ1 + 2 α−

∫ ∫ Γb

l/ε

−l/ε

φˆ b dz ds + O(ε 2 ).

(2.58)

|Γb | = γb,0 + εγb,1 + O(ε 2 ).

(2.59)

Evaluating the integrals in Eq. (2.58) and solving for µ1 yields the expression

) 2 ( α− ˆ − γb,0 mb + O(ε ), M |Ω |

(2.60)

where mb is defined in (2.51). On the other hand, the area of a smooth curve subject to normal velocity V evolves according to

∂|Γb | = ∂τ

∫ Γb

V (s)H0 (s) ds,

(2.61)

and for the normal velocity (2.53) this reduces to the result presented in (1.12). Taking the time derivative of (2.60), using (2.61) to eliminate ddτ γb,0 , the normal velocity (2.53) drives the bulk chemical potential according to dµ1 dτ

=−

2 α− mb |Ω |

(

νb (µ1 − µ∗b )

where N1 (s) and N2 (s) are orthogonal unit vectors which are also orthogonal to the tangent vector ζf′ (s), defined by

∂ Ni = −κi T, ∂s

i = 1, 2,

(3.2)

where

κ⃗ (s, t) := (κ1 , κ2 )t ,

(3.3)

is the normal curvature vector with respect to {N , N }. 1

)

Our choice of initial data implies that the mass can be rescaled as ˆ + O(ε 2 ). We also expand the surface area M = εM

µ1 =

9

∫ Γb

H02 (s) ds − ε kb

∫ Γb

) |∇s H0 |2 ds

+ O(ε 2 ).

(2.62)

The coupled system (2.53) and (2.62) prescribes the interfacial evolution for the dressing of an admissible codimension one interface with a shifted bulk value of u given by (2.55). 3. Geometric evolution of codimension two structures

f

Definition 3.1. For any K , ℓ > 0 the family, GK ,ℓ , of admissible curves is comprised of closed (compact and without boundary), oriented 1 dimensional curves Γf embedded in R3 , which are far from self intersection and with a smooth second fundamental form. More precisely, (i) The W 4,∞ (Sf ) norm of the 2nd Fundamental form of Γf and its principal curvatures are bounded by K . (ii) The whiskers of length 3ℓ < 1/K , in the unscaled distance, defined for each s0 ∈ Sf by, ws0 := {x : s(x) = s0 , |z(x)| < 3ℓ/ε}, neither intersect each-other nor ∂ Ω (except when considering periodic boundary conditions). (iii) The length, |Γf |, of Γf is bounded by K /ε . For an admissible codimension two curve Γf the change of variables x → ρ (s, z) given by (3.1) is a C 4 diffeomorphism on the reach of Γf , defined as the set

⏐ } { ⏐ Γf ,ℓ := ρ (s, z) ∈ R3 ⏐s ∈ S , −ℓ/ε ≤ |z| ≤ ℓ/ε ⊂ Ω .

ε 2 ∆x = ∆z − ε⃗ κ · ∇z + ε 2 (∂s2 − (z · κ⃗ )κ⃗ · ∇z ) + O(ε 3 ),

Given a smooth, closed, non-self intersecting one-dimensional manifold Γf immersed in Ω ⊂ R3 , and parameterized by the map s ∈ Sf ↦ → ζf (s) ∈ Ω , we may uniquely decompose points x near Γf as x = ρf (s, z1 , z2 ) = ζf (s) + ε (z1 N1 (s) + z2 N2 (s)) ,

(3.1)

(3.5)

where the lower order terms are immaterial for the analysis. The Jacobian of the change of variables (3.1) takes the form J = ε 2 − ε 2 z · κ⃗ .

(3.6)

If the underlying curve Γf evolves in time, then its normal velocity vector V = (V1 , V2 ) of Γf at a point s(t) takes the form

∂ z1 + ε z2 N 2 · ∂t ∂ z2 V2 := −ε + ε z1 N 1 · ∂t

∂N1 , ∂t 2 ∂N . ∂t

(3.7) (3.8)

The terms z2 N 2 · ∂∂Nt and z1 N 1 · ∂∂Nt reflect lower-order contributions to the normal velocity induced by the rotational motion of the normal vectors to Γf (t). See [16] for further details. 1

2

Definition 3.2. Given an admissible codimension two curve Γf ∈ f GK ,ℓ and a smooth function f : R+ → R which tends to a constant value f∞ at an O(1) exponential rate as R → ∞, we define fΓf ∈ H 2 (Ω ), called the dressing of Γf with f , according to the rule fΓf (x) := f (R(x))χ (R(x)/ℓ) + f∞ (1 − χ (R(x)/ℓ)),

3.1. Admissible codimension two curves and their dressings

(3.4)

where z := (z1 , z2 ). We introduce R(x) = |z(x)| which denotes the scaled distance of x to Γf . Within the reach the cartesian Laplacian admits the local form

V1 := −ε

In this section we derive the geometric evolution of admissible codimension two curves, called filaments, in R3 under the H −1 gradient flow (1.3). As remarked in the introduction, we make the assumption that the combined length of all the filament curves scales as O(ε −1 ) so that the combined mass of the codimension two structures is O(ε ). This gives a comparable mass to filament and bilayer structures, so they may contribute to the mass balance at the same order of magnitude. Codimension two structures are much less studied than codimension one structures, however our analysis leads to a qualitatively similar result: a surface diffusion regularized curvature-vector driven normal flow that may be curve lengthening or curve shortening depending upon the value of the spatially constant far-field chemical potential.

2

(3.9)

where χ : R → R is a fixed, smooth cut-off function taking values one on [0, 1], and zero on [2, ∞). By abuse of notation we will drop the Γf subscript when doing so creates no confusion. Within the reach Γf ,ℓ the Cartesian Laplacian reduces formally at leading order to the two-dimensional Laplacian in z, which may be written in turn in polar coordinates in R. We may eliminate the dominant terms in quadratic component of the (1.1) by taking

10

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

u at leading order to be the dressing of the codimension two profile φf (R), defined as the solution of 1

∂R2 φf + ∂R φf − W ′ (φf ) = 0,

(3.10) R subject to the boundary conditions ∂R φf (0) = 0 and φf → b− as R → ∞. We denote the dressing of Γf with φf by Uf . As in the codimension one case, we introduce the operator 1 ∂R − W ′′ (φf ), (3.11) R corresponding to the linearization of (3.10) about φf . The operator Lf ,0 is the radially symmetric reduction of the associated cylindrical Laplacian,

Lf ,0 := ∂R2 +

1

1

The time derivative of the inner density function u˜ , defined in (3.18), combined with the normal velocity, (3.7), takes the form

∂ u˜ ∂ u˜ ∂τ = −ε −1 V · ∇z u˜ + . ∂t ∂τ ∂ t

(3.20)

For a whisker identified by s ∈ Sf , with base point x = ρf (s, 0) ∈ Γf we choose vectors n, m ∈ span{N1 , N2 } in the normal plane of Γf at x, and choose θ so that n = cos(θ )N1 + sin(θ )N2 .

(3.21)

The usual directional derivative along n is denoted

∂n := n · ∇x = cos(θ )N1 · ∇x + sin(θ )N2 · ∇x ,

(3.22)

(3.12)

and for f ∈ C ∞ (Ω /Γf ) we introduce the n, m limit

This operator is self-adjoint in the usual R-weighted L2 (R+ ) inner product,

∂nj f m (x) := lim (n · ∇x )j f (x + hm, t) for all j ≥ 0,

Lf := ∂R2 +

⟨f , g ⟩R :=

R

∂R +

R2

∂θ2 − W ′′ (Uf ).

∞

∫

f (R)g(R)R dR.

(3.13)

0

Moreover, the translational eigenfunctions {φf (R) cos(θ ), φf (R) sin(θ )} lie in the kernel of Lf and their associated dressings of Γf agree with {∂z1 Uf , ∂z2 Uf } respectively, up to exponentially small terms. For each m ∈ N, we define the spaces ′

⏐ := {f (R) cos(mθ ) + g(R) sin(mθ ) ⏐ f , g ∈ C ∞ (0, ∞), m ∈ N}. (3.14)

These spaces are invariant under the operator Lf , and mutually orthogonal in L2 (Ω ). Moreover, on these spaces Lf reduces to Lf (f (R) cos(mθ ) + g(R) sin(mθ )) = cos(mθ )Lf ,m f + sin(mθ )Lf ,m g , (3.15) where Lf ,m :=

∂2 1 ∂ m2 + − 2 − W ′′ (Uf ). 2 ∂R R ∂R R

(3.23)

and the limit of the gradient

∇x f m (x) := lim ∇x f (x + hm, t),

(3.24)

h→0+

′

Zm

h→0+

where the limit exists. If f ∈ C 1 (Ω ) then the normal derivative of f will satisfy

∂n f −m = ∂n f m .

(3.25)

This motivates the following definition of the jump condition. Definition 3.3. Given a radial function f := f (R) localized on Γf , we define the jump of f across a given whisker by J∂n f m KΓf (x) := ∂n f m (x) − ∂n f −m (x)

(3.26)

which is zero when f has a smooth extension through Γf . With this notation we develop matching conditions

(3.16)

Each operator Lf ,m is self-adjoint in the R-weighted inner product, and the operator Lf ,1 has a 1-dimensional kernel spanned by its ground state ∂R φf > 0. For m > 1 we observe that ⟨Lf ,m f , f ⟩R < ⟨Lf ,1 f , f ⟩R and since Lf ,1 ≤ 0 we deduce that Lf ,m < 0 for m > 1, and is boundedly invertible. We denote the eigenfunctions and ∞ eigenvalues of Lf ,m by {ψf ,m,j }∞ j=0 and {λf ,m,j }j=0 , respectively. We address the kernel of Lf ,0 with the following assumption. Assumption 1. We assume that the operator Lf ,0 has no kernel and a one-dimensional positive eigenspace spanned by ψf ,0,0 . With this assumption we may define the functions

lim u(x + ε Rn, t) = lim u˜ (s, S , R, θ, t).

R→0+

(3.27)

R→∞

Expanding the left-hand side yields the following expression u(x + ε Rn) =un0 (x, t) + ε un1 (x, t) + R∂n un0 (x, t)

(

)

) ( 1 + ε 2 un2 (x, t) + R∂n un1 (x, t) + R2 ∂n2 un0 (x, t) 2

+ O(ε 3 ),

(3.28)

where uni denotes the limit of the left-hand side of (3.27) as ε R → 0+ . Equating orders of ε for the matching condition (3.27) yields un0 = lim u˜ 0 (s, S , R, θ, τ ),

(3.29)

for j = 1, 2 and their Γf dressings, also denoted Φf ,1 and Φf ,2 .

un1 + R∂n un0 = lim u˜ 1 (s, S , R, θ, τ ).

(3.30)

3.2. Inner and outer expansions

The chemical potential, defined in (1.2) admits similar inner and outer expansions. The outer expansion is identical to that for the codimension one case, see (2.16) and (2.17). To obtain the inner expression for the chemical potential we first note that in the multiscale tangential variables the Laplacian expression (3.5) takes the form

−j

Φf ,j := Lf ,0 1,

(3.17)

Considering initial data that is close to a filament dressing of f an admissible curve, Γf (0) ∈ GK ,ℓ , embedded in Ω ⊂ R3 . In the far-field, Γ˜ f ,ℓ , the outer solution u has the expansion u(x, t) = u0 (x, t) + ε u1 (x, t) + ε 2 u2 (x, t) + O(ε 3 ),

(3.18)

and within the reach Γf ,ℓ , we incorporate a two-scale tangential expansion, introducing the variable S = √sε , and the inner spatial expansion takes the form u(x, t) = u˜ (s, S , z, τ ) = u˜ 0 (s, S , z, τ )

+ ε u˜ 1 (s, S , z, τ ) + ε 2 u˜ 2 (s, S , z, τ ) + O(ε 3 ).

(3.19)

R→∞ R→∞

( ) ε 2 ∆x = ∆z −ε κ⃗ · ∇z + ∂S2 +ε 2 (∂s2 −(z·⃗ κ )κ⃗ ·∇z +2z · κ⃗ ∂S2 )+O(ε 3 ), (3.31) Introducing the nonlinear operators P and Q , the inner chemical potential is written as

µ ˜ (x, t) = P(u˜ )Q (u˜ ) + εηd W ′ (u˜ 0 )

(3.32)

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

where P admits the expansion P = P0 + ε P1 + ε 2 P2 + · · ·, with P0 = − ∆z + W (u˜ 0 ), ′′

(3.33)

P1 =⃗ κ · ∇z + ∂ + W (u˜ 0 )u˜ 1 − η1 , 2 S

′′′

(3.34)

P2 = − (z · κ⃗ )κ⃗ · ∇z + 2z · κ⃗ ∂ − ∂ + W (u˜ 0 )u˜ 2 + 2 S

2 s

′′′

1 2

W

(4)

˜ ˜ ,

(u0 )u21

(3.35)

Q0 = − ∆z u˜ 0 + W ′ (u˜ 0 ),

(3.36)

Q1 =⃗ κ · ∇z u˜ 0 + ∂S2 u˜ 0 + −∆z + W ′′ (u˜ 0 ) u˜ 1 ,

(

)

(3.37)

Q2 = − (z · κ⃗ )κ⃗ · ∇z u˜ 0 − 2z · κ⃗ ∂ ˜ − ∂ ˜ 2 S u0

2 s u0

(u0 )u21

µ ˜ 0 =P0 Q0

where µ1 is a spatial constant and Φf ,2 is defined in (3.17). With this simplification equation (3.48) becomes

where

ξ (s, S , R, θ ) :=

µ ˜ 2 =P0 Q2 + P1 Q1 + P2 Q0 + ηd W (u˜ 0 )u˜ 1 .

(3.41)

The relevant matching conditions for the chemical potential extend to second order in ε :

µ

n 0 (x

µ

, t) = lim µ ˜ 0 (s, S , R, θ, t), R→∞

(3.53) We observe that V · ∇z Uf = ∂R Uf (R) (V1 cos θ + V2 sin θ) ∈ Z1 ,

∂ −

1 R R

1

µn2 (x, t) + R∂n µn1 (x, t) + R2 ∂n2 µn0 (x, t) = lim µ ˜ 2 (s, S , R, θ, t), 2

R→∞

(3.44) 3.3. Time scale τ = ε t: quenched vector-curvature driven flow The analysis of the outer expansion of the chemical potential is identical to the codimension one case, and we find at leading order that u0 = b− , µ0 = 0, while at O(ε ) we obtain

∆x µ1 = 0 in Ω \Γf .

(3.45)

In the inner region we supplement the inner expansions (3.19) and (3.32) with the inner expression of the Laplacian (3.5). At leading order in ε we find 0 = ∆z µ ˜ 0,

in Γf ,ℓ ,

(3.46)

where µ ˜ 0 is defined in (3.39). This equation is consistent with the choice of initial data u˜ 0 = φf which implies that µ ˜ 0 = 0 via the matching conditions. With this reduction the subsequent orders become O(ε

−1

) : 0 = ∆z µ ˜ 1,

in Γf ,ℓ , (3.47)

O(1) :

( ) − V · ∇z u˜ 0 = ∆z µ ˜ 2 − κ⃗ · ∇z + ∂S2 µ ˜ 1,

in Γf ,ℓ . (3.48)

where µ ˜ 1 and µ ˜ 2 are defined in Eqs. (3.40) and (3.41), respectively. The combined system (3.45) and (3.47) couples through the matching condition (3.43). Since µ0 = 0 we deduce from (3.43) that µ ˜ 1 is bounded as R → ∞, and hence from (3.47) that µ ˜ 1 is constant in z. In particular µ ˜1 = µ ˜ 1 (s, S , τ ) ≈ µ1 . Since u˜ 0 = Uf , Eq. (3.40) for µ ˜ 1 reduces to a linear equation for u˜ 1 , L2f u˜ 1 = µ ˜ 1 − ηd W ′ (Uf ).

(3.49)

arrive at the system 1

˜ 1 (s, S), ∂R2 C + ∂R C = ∂S2 µ

, t) + R∂n µ0 (x, t) = lim µ ˜ 1 (s, S , R, θ, t), R→∞

(3.54)

while ∂ µ ˜ ∈ Z0 . We project (3.51) onto Zm where ∆z = ∂R2 + 2 S 1 m2 and R2

1

R 1

R 1

R2 1

±

(3.43)

O(ε −2 ) :

(Am (s, S , R) cos(mθ ) + Bm (s, S , R) sin(mθ )) .

m=2

(3.42) n 1 (x

∞ ∑

(3.40) ′′

(3.51)

To impose interfacial matching conditions for µ ˜ 2 we solve (3.51) by expanding µ ˜ 2 in (R, θ ) inner-polar coordinates associated to the spaces {Zm }∞ m=0 as

(3.39)

µ ˜ 1 =P1 Q0 + P0 Q1 + ηd W ′ (u˜ 0 ),

(3.50)

(3.52)

W ˜ ˜ . (3.38) 2 With these reductions we expand the inner chemical potential as

−

′′′

2 ′ u˜ 1 = µ1 Φf ,2 − ηd L− f W (Uf ),

µ ˜ 2 = A1 (s, S , R) cos θ + B1 (s, S , R) sin θ + C¯ (s, S , R) + ξ (s, S , R, θ ),

+ (κ⃗ · ∇z + ∂S2 )u˜ 1 + (−∆z + W ′′ (u˜ 0 ))u˜ 2 1

By Assumption 1 we know that ker Lf ⊂ ker Lf ,1 ⊂ Z1 , defined in (3.14), while the right-hand side of (3.49) lies in Z0 . Since the spaces Zm are mutually orthogonal, we may solve for u˜ 1 ,

O(1) : −V · ∇z Uf = ∆z µ ˜ 2 − ∂S2 µ ˜ 1.

and similarly Q = Q0 + ε Q1 + ε Q2 + · · ·, where 2

11

∂R2 A1 + ∂R A1 + ∂R2 B1 + ∂R B1 +

A1 = V1 (s, S)∂R Uf (R),

(3.55) (3.56)

B1 = V2 (s, S)∂R Uf (R), (3.57) R R2 plus homogeneous equations for {Am , Bm }∞ m=2 that have nonsingular solutions Am = am (s, S)Rm and Bm = bm (s, S)Rm . The Z0 equation has solution C (s, S , R) = C0 (s, S) +

R2

∂ 2µ ˜ 1 (s, S), (3.58) 4 S while the non-homogeneous equations, (3.56) and (3.57), have the solutions A1 (s, S , R) = a1 (s, S)R − a(R)V1 (s, S),

(3.59)

B1 (s, S , R) = b1 (s, S)R − a(R)V2 (s, S),

(3.60)

where a(R) is the solution of the non-homogeneous ordinary differential equation 1 ′ 1 a − 2 a = ∂R Uf (R), R R which enjoys the explicit formula a′′ +

a(R) =

1 R

(3.61)

R

∫

r Uˆ f (r) dr ,

(3.62)

0

where we have introduced Uˆ f := Uf − b− . In particular, a(R) → 0 as R → ∞. From the matching condition (3.44) we see that µ ˜ 2 grows at most linearly as R → ∞ and

∂µ ˜2 = ∂n µn1 = cos θ N1 · ∇x µn1 + sin θ N2 · ∇x µn1 , R→∞ ∂ R lim

(3.63)

where the second equality follows from the definition of the directional derivative along n, given in (3.22). Taking the R derivative of (3.52) and using the results above yields

∂µ ˜2 R ∂ξ = ∂S2 µ ˜ 1 + (a1 − a′ (R)V1 ) cos θ + (a2 − a′ (R)V2 ) sin θ + . (3.64) ∂R 2 ∂R

12

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

Projecting each of (3.63) and (3.64) onto Zm and matching terms, we conclude that ξ = 0 and ∂S2 µ ˜ 1 = 0, and in particular

∂µ ˜2 ∂µ ˜2 (s, S , R, θ , τ ) = − (s, S , R, θ + π, τ ). ∂R ∂R

J∂n µ

= 0,

(3.66)

for any choice of normal vector n. As the codimension two curve Γf has zero capacity, it follows from the zero-jump condition and (3.45) that µ1 has a harmonic extension to all of Ω , and hence is spatially constant. In particular ∇x µn1 = 0 for all choices of direction n. With these reductions µ ˜ 2 takes the form

µ ˜ 2 = C0 (s, S) − a(R) (V1 (s, S) cos θ + V2 (s, S) sin θ) .

(3.67)

To determine the normal velocity we substitute u˜ 0 = Uf into the expression (3.41) for µ ˜ 2 , so that P0 reduces to Lf and the chemical potential takes the form

µ ˜ 2 =L2f u˜ 2 − Lf Q˜ 2 + (κ⃗ · ∇z + ∂S2 + W ′′′ (u˜ 0 )u˜ 1 − η1 ) × (−Lf u˜ 1 + κ⃗ · ∇z Uf ) + ηd W ′′ (Uf )u˜ 1 ,

(3.68)

where we have introduced Q˜ 2 := Q2 − Lf u˜ 2 . To solve equation (3.68) for u˜ 2 we rewrite it in the form

˜ 2 − Q + Lf Q˜ 2 , L2f u˜ 2 = µ

Uˆ f2 RdR.

(3.79)

0

Substituting (3.76) and (3.78) into (3.75) we arrive at the expression for the normal velocity V(s, S) = −

2µ1 mf − η1 σf mf ,2

σf

κ⃗ −

m f ,2

∂S2 κ⃗ .

η 1 σf

µ∗f :=

2mf

,

νf :=

2mf mf ,2

,

kf =

Q := (κ ⃗ · ∇z + ∂ + W (u˜ 0 )u˜ 1 − η1 )(−Lf u˜ 1 + κ⃗ · ∇z Uf )

+ ηd W ′′ (Uf )u˜ 1 .

V(s, S) = − νf (µ1 − µ∗f )κ⃗ + ε kf ∂s2 κ⃗ .

]

[

u(x, t) = b− + ε

µ1 + O(ε 2 ) in Γ˜ f ,ℓ . 2 α−

(3.84)

M := Ω

u(x, t) − b− dx =

∫ ∫Ω

u(x, 0) − b− dx

∫ (u − b− ) dx +

= Ω \Γf ,ℓ

Γf ,ℓ

where the outer integral takes the value

∫

Q = Q0 + Q1 + Q0,2 ,

Ω \Γf ,ℓ

(3.71)

(u − b− ) dx = ε

∫

Q0 := −W ′′′ (Uf )u˜ 1 Lf u˜ 1 + ∂S2 Lf u˜ 1 + η1 Lf u˜ 1 + ηd W ′′ (Uf )u˜ 1 ,

(3.72)

Γf ,ℓ

(u − b− ) dx = ε 2

(3.73) (3.74)

Since the spaces Zm are orthogonal and Lf (Q˜ 2 )⊥ ker Lf , the solvability conditions take the form for i = 1, 2.

(3.75)

To evaluate these conditions we expand Q1 using the expression (3.50) for u˜ 1 ,

( ) ⟨Q1 , ∂zi Uf ⟩R = −2π mf µ1 κi − πσf η1 κi + ∂S2 κi ,

for i = 1, 2, (3.76)

σf :=

∞

∫

(Uf′ )2 R dR.

Γf

( R2

(3.86)

2 ′ Uˆ f + ε (µ1 Φf ,2 − ηd L− f W (φf ))

(3.77)

0

The R-weight inner product of µ ˜ 2 , given in (3.67), with ∂zi Uf yields (3.78)

(3.87)

= ε 2π |Γf |mf + O(ε |Γf |). 2

3

Assuming that |Γf | ∼ O(ε −1 ), as is commensurate with an O(1) amphiphilic mass, we expand

|Γf | = ε −1 γf ,−1 + γf ,0 + O(ε ),

(3.88)

to arrive at the total mass expansion M=ε

(

) µ1 |Ω | + 2π mf γf ,−1 + O(ε 2 ). 2 α−

(3.89)

Taking the τ = ε t time derivative of the total mass, (3.89), and dγ 1 solving for df ,− yields τ dγf ,−1

where we have introduced ∞

∫ ∫

) +O(ε 2 ) (1 − ε z · κ⃗ ) dz ds

Q1 := −⃗ κ · ∇z Lf u˜ 1 + W ′′′ (Uf )u˜ 1 κ⃗ · ∇z Uf − ∂S2 κ⃗ · ∇z Uf Q0,2 := (κ ⃗ · ∇z )2 Uf .

µ1 (|Ω | − |Γf ,ℓ |) + O(ε 2 ). 2 α−

Using (3.83) and the Jacobian, (3.6), we evaluate the inner integral

where Q0 ∈ Z0 , Q1 ∈ Z1 , Q0,2 ∈ Z0 + Z2 , are given by

− η1 κ⃗ · ∇z Uf ,

(u − b− ) dx, (3.85)

We decompose Q into its Zm components

⟨µ ˜ 2 , ∂zi Uf ⟩R = π Vi mf ,2 ,

in Γf ,ℓ , (3.83)

which has the far-field asymptotics,

(3.70)

0

(3.81)

(3.82)

2 ′ 2 u(x, t) = Uf + ε (µ1 Φf ,2 − ηd L− f W (Uf )) + O(ε )

ker Lf = span{∂R Uf cos θ , ∂R Uf sin θ} = ker Lf ,1 .

Uˆ f RdR,

,

The constant value of u1 is determined by the conservation of total mass, and is coupled to changes in length of the curve Γf . Combining the inner and outer expansions of u yields the composite expansion

For fixed values of s and S, the expression (3.69) can be solved for u˜ 2 if and only if the right-hand side is perpendicular to

mf :=

m f ,2

and return the S variable to its original scaling, we obtain the normal velocity

∫

∫

σf

The total mass of the system is given by

′′′

⟩ ⟨ µ ˜ 2 − Q1 , ∂zi Uf R = 0,

(3.80)

Introducing the quantities

(3.69)

where 2 S

∞

∫ mf ,2 :=

(3.65)

This latter result, substituted into (3.63) yields the no-jump condition across the curve Γf n 1 KΓf

where we have introduced

dτ

=−

|Ω | dµ1 . 2 2πα− mf dτ

(3.90)

On the other hand, any admissible codimension two curve evolving with normal velocity V satisfies d|Γf | dτ

∫ =− Γf

V · κ⃗ ds.

(3.91)

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

Combining this expression with (3.82), (3.88), and (3.90) yields dµ1

2π α− mf 2

=ε

dτ

|Ω | ∫

+ε kf

Γf

( ) ( −νf µ1 − µ∗f

∫ Γf

|⃗ κ |2 ds

) |∂s κ⃗ | ds + O(ε 2 ). 2

(3.92)

This system exhibits the same quenching behavior as the codimension one evolution, with the distinction being that the equilibrium far-field density for a codimension two curve takes the form lim u = b− + ε

τ →∞

µ∗f α− 2

+ O(ε 2 ).

(3.93)

13

pearling bifurcation, Pb and Pf , depend upon the functionalization parameters η1 and η2 through their difference ηd := η1 − η2 . The µ1 thresholds for the transition from curve shortening to regularized curve lengthening, µ∗b and µ∗f given in (2.52) and (3.81) respectively, have a more subtle dependence upon the functionalization parameters. These relations are summarized below: µ1 sign(Sb ) < Pb (ηd ) := −ηd

λb,0 ∥ψb,0 ∥22 L

|Sb |

,

Bilayers Pearling

(4.5)

Stable

µ1 sign(Sf ) < Pf (ηd ) := −ηd

∥ψf′ ,0,0 ∥22 + λf ,0,0 ∥ψf ,0,0 ∥22 L R

L R

|Sf |

,

Filaments Pearling

(4.6)

Stable

4. Competitive evolution of amphiphilic suspensions

µ1 < µ∗b (η1 , η2 ),

We assume that Ω ⊂ R3 and combine the geometric flow results derived in Sections 2 and 3 with the pearling stability results for bilayers and filaments presented in [6]. The goal is to derive an overall picture of the complexity of transients and bifurcation structure of the H −1 gradient flow of the strong scaling of the FCH system. 4.1. Competitive evolution of codimension one and two systems Fix Ω ⊂ R3 and let Γb and Γf be admissible codimension one and codimension two morphologies with disjoint reaches, Γb,ℓ and Γf ,ℓ . We emphasize that Γb and Γf may be comprised of multiple disjoint surfaces and curves. For a given value of the chemical potential, µ1 , the codimension one bilayer morphology ub , given in (2.54) and codimension two filament morphology uf , given in (3.83) satisfy identical far-field asymptotics lim uf = lim ub = b− − ε z →∞

R→∞

µ1 + O(ε 2 ). 2 α−

(4.1)

Consequently we may form the composite solution

µ1 = ub + uf − b− − ε 2 α−

(

ub,f

)

+ O(ε 2 ),

(4.2)

parameterized by Γb , Γf , and the common, slowly varying, chemical potential µ1 . Recalling the scalings (2.59) and (3.88) of the surface area and length of Γb and Γf respectively, the total mass of the composite solution satisfies M=ε

(

) µ1 | Ω | + m γ + 2 π m γ + O(ε 2 ), b b,0 f f ,−1 2 α−

(4.3)

where mb , the bilayer mass per unit area, is defined in (2.51), and 2π mf denotes the filament mass per unit length, defined in (3.77). ˆ + O(ε 2 ) and solving for µ1 yields relation Expanding M = ε M between the morphology size and the chemical potential µ1 ,

µ1 =

) 2 ( α− ˆ − mb γb,0 − 2π mf γf ,−1 . M |Ω |

(4.4)

Taking the time derivative of (4.4) and using the relations (2.61) and (3.91) to relate the growth of the curves to the normal velocities, yields an evolution equation for the chemical potential given in (1.9). Coupling this equation to the normal velocities for the bilayer and filament derived in (1.7) and (1.8) gives a closed system for the combined curve motion and far-field chemical potential. 4.2. Analysis of competitive geometric evolution and bifurcation The analysis presents a bifurcation diagram with four thresholds that delineate distinct behaviors. The µ1 thresholds for the

Bilayers Curve

(4.7)

Shortening

µ1 < µ∗f (η1 , η2 ),

Filament Curve Shortening.

(4.8)

The signs of the shape factors Sb and Sf , defined in (1.6), depend upon the choice of the double well, W , and impact not only the sign of the right-hand sides of (4.5) and (4.6) but also the direction of the inequalities, see Fig. 3 (right). The chemical potential and the geometric flows evolve on the same t = O(ε −1 ) timescale. Within the H −1 gradient flow, the pearling instability produces eigenvalues of size O(ε −1 ), [6] and hence will manifest itself on the t = O(ε ) timescale, essentially instantaneously on the time scale of the geometric flow and the chemical potential. We define the pearling instability region to be set of values (µ1 , η1 , η2 ) for which either codimension one or codimension two structures are pearling unstable. We investigate the variation of the pearling instability regions, and its relation to regions of curve lengthening flows, as functions of the well shape. For simplicity we present these regions in the µ1 –ηd plane, with the assumption that η1 = 0.15, unless specified otherwise. To parameterize the well shape we fix b± = ±1 and insert a one-parameter well tilt, ξ into the double-well potential,

W (u; ξ ) :=

1 2

(u − b− )2

(

1 2

(u − b+ )2 −

ξ 3

( u−

3b+ − b− 2

))

,

(4.9) where the parameter ξ controls the value of W at the right well u = b+ , see Fig. 3 (left). For the µ1 dynamics, as η1 and η2 are fixed parameters, the temporal evolution of µ1 traces a vertical line segment on the diagram, with µ1 decreasing if it is larger than both µ∗b and µ∗f and increasing if it is smaller than both. The µ1 -region bounded by the points {µ∗b , µ∗f } is attracting and forward invariant under the flow see Fig. 3 (center). Within this invariant region the direction of the flow depends upon the overall curvatures of the two classes of morphologies; however to leading order the total area/length of the codimensional structures with the larger equilibrium value decreases monotonically and the other increases monotonically, so long as the morphologies remain admissible. A simulation of the FCH gradient flow with one spherical bilayer and two circular filaments is presented in Fig. 4 for parameter values for which µ∗b < µ∗f . Motivated by this example, and to further illustrate the dynamics, we consider a composite morphology consisting of Nb spherical bilayers of radii R1 , . . . , RNb and Nf = O(ε −1 ) circular filaments of radii r1 , . . . , rNf . 1 For these special shapes the

14

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

Fig. 3. (left) Graph of the double well, W versus u for ξ = −0.9, −0.7, −0.5 (see (4.9)). (center) Diagram of the curve shortening regions in ηd –µ1 plane, for ξ = −0.9 and η1 = 0.15. Arrows indicate the direction of the temporal evolution of µ1 under the flow (1.9). The shaded region is forward invariant and globally attracting so long as the curves remain admissible. (right) Values of the bilayer and filament shape factors, Sb (blue) and Sf (red) as a function of the well tilt, ξ . The change in sign of Sb near ξ = −0.8 flips the direction of the inequality in (4.5) . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. Simulation of the FCH gradient flow for ε = 0.03, η1 = 0.15, η2 = 0.24, and ξ = −0.15. For these values µ∗b < 0 < µ∗f and after a short transient bilayers will grow while filaments shrink. (a) Initial data consisting of two circular filaments and a spherical bilayer inside the [−π, π]3 computational domain. (b) 2D slice along x–z plane (y = 0) of initial data. (c) 3D pose of the t = 200 (τ = 6) final computation stage and (d) corresponding 2D slice showing the larger bilayer and smaller filament radii.

competitive evolution (1.7)–(1.9) reduces to

( ) 2 R˙ i = νb µ1 − µ∗b ,

i = 1, . . . , Nb ,

(4.10)

r˙j = νf µ1 − µf

j = 1, . . . , ε −1 Nf ,

(4.11)

Ri ) ∗ 1

(

rj

,

⎛

⎞ Nf 2 ∑ α− 1 ⎝16π mb νb (µ1 − µ∗b )Nb + ε4π 2 mf νf (µ1 − µ∗f ) ⎠, µ ˙1 = − |Ω | rj j=1

(4.12) where the dot notation denotes differentiation with respect to τ = t /ε . The µ1 evolution depends upon the spherical bilayers only through their total number. Consistent with the discussion above, the bulk chemical potential µ1 decreases if µ1 > max{µ∗b , µ∗f } and increases if µ1 < min{µ∗b , µ∗f }. The radii shrink or grow depending upon the signs of µ1 − µ∗b and µ1 − µ∗f . If both Nb and Nf are positive, and µ∗b ̸ = µ∗f then the system has no equilibrium. Assuming for simplicity of presentation that µ∗b < µ∗f , then after a possible transient the system enters a regime in which all spheres are growing and all circular filaments are shrinking. The radii rj will then decrease to zero in finite time due to the inverse relation between r˙j and rj . If the zero radius filaments are removed from the system and the remaining filaments re-indexed, then after a transient the set of hoops will be empty (Nf = 0) and only spheres will remain. At this point µ1 will relax (quench) at an exponential rate from above to µ∗b as the spherical radii grow, albeit at an exponentially decreasing rate as µ1 quenches to µ∗b . The equilibrium will be a collection

of spherical bilayers of differing radii. We emphasize that in the growing regime spherical shapes are unstable to perturbation under the full flow (1.7). Manifestation of this instability requires sufficiently large values of νb (µ1 − µ∗b ) > 0 in relation to the coefficient ε kb of the surface diffusion term in (1.7), and is triggered more easily with increasing radius. However, if the non-spherical excursions are not so large as to induce self-intersection, then the shapes return to spherical as µ1 quenches to µ∗b . Fig. 5 presents a simulation of the full FCH gradient flow that illustrates this phenomenon, while a rigorous derivation of the transient instability regime of the curve lengthening flow in 1 + 2D for nearly circular bilayers is presented in [20]. In Fig. 6, the pearling bifurcation lines are added and the full stability diagram is shown for values of the well-tilt parameter ξ = −0.85, −0.7, −0.45, and −0.2. For ξ = −0.85, Sf is positive, Sb is slightly negative, and for positive values of µ1 filament pearling instability region contains the bilayer pearling instability region. As ξ is decreased from −0.7 to −0.45 the two pearling instability lines almost coincide, and for ξ = −0.2 they have crossed, with the filament pearling instability region now lying above the bilayer pearling instability region for µ1 > 0. In all figures the intersections of the pearling instability lines occur to the left of the crossing of the curve lengthening lines. This implies that increasing ηd from negative values will excite pearling bifurcations at smaller values of ηd than those for which the filament morphology becomes dynamically favored over the bilayer morphology. For the non-generic value of ηd at which the critical values µ∗b = µ∗f coincide, codimension one and two

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

15

Fig. 5. Simulation of the FCH gradient flow corresponding to a single circular bilayer in [−2π, 2π]2 ⊂ R2 in the regularized curve lengthening regime. The larger radii and value of µ1 − µ∗b , compared to simulations presented in Fig. 4, induce the interfacial meander. Yet larger initial values of µ1 − µ∗b lead to self-intersection and defect generation. (a) The initial data has µ1 > µ∗b , with sufficient excess to initiate shape perturbations as the circular bilayer grows. (b) The onset of the shape instability at time t = 484. (c) Quenching of the flow at t = 986 as the value of µ1 relaxes towards µb ∗ and the higher order surface diffusion returns the interface back to a larger, circular shape.

structures can co-exist on the t = O(ε −1 ) time-scales under consideration here. 4.3. Analytic bifurcation diagrams and comparison to simulations We compare the evolution of numerical simulations of the H −1 gradient flow of the FCH free energy, (1.3), to the corresponding bifurcation results and to simulations from a self-consistent mean field model. The parameters impacting stability are the functionalization parameters η1 and η2 , the time-dependent scaled chemical potential µ1 , and the shape of the double well, parameterized by ξ . For simplicity we fix η1 = 0.15 and vary η2 and ξ . Simulations of the strong FCH gradient flow, (1.1), were conducted in a domain Ω = [−π, π]3 ⊂ R3 , with initial data of the form (4.2) with Γb consisting of a single sphere of radius R1 = 1.1 with center at (0, 0, 0) and Γf comprised of two circular filaments of radius r1 = r2 = 1.6 oriented parallel to the ) x–y ( plane of Ω and with center point at 0, 0, ± 12 (π + R1 /2) . The initial value of µ1 varied with each simulation and is reported in Fig. 7. The system parameters are ε = 0.03, η1 = 0.15 and η2 = 0.24, hence ηd = −0.09, and four separate values of ξ := −0.15, −0.2, −0.25, −0.3. In addition, there was one simulation for ξ = −0.2 and η1 = η2 = 0.15. The simulations were conducted for t ∈ [0, 100], equivalently τ ∈ [0, 3]. Results for the five simulations are superimposed upon the corresponding µ1 − ηd bifurcation diagram and presented in Fig. 8. The initial and final values of µ1 for each simulation are indicated with a closed circle and closed square respectively in each of the four diagrams. For ξ = −0.15 the value of µ1 starts in a region of bilayer and filament pearling stability, filament curve shortening, and bilayer regularized curve lengthening, see Fig. 8(bottom-right). During the simulation the bilayer radius grew, the filaments shrunk, and neither pearled, the τ = 3 end state is presented in Fig. 7 (left/top-left). Two simulations were conducted for ξ = −0.2, for the simulation with ηd = −0.09 the initial value of µ1 lies at the boarder of the bilayer pearling region, and the initial stages of the simulation (t < 5) displayed the onset of pearling, however the value of µ1 decreased out of the pearling region as the filaments and bilayers grow and the pearling evanescent, restoring the unpearled bilayer structure. The end-state is presented in Fig. 7 (left/top-right), there is less shrinking of the circular filaments than in the case η2 = −0.15 and the filaments are thinner due to the stronger well tilt. The ξ = −0.2 simulation with ηd = 0 but an identical initial value of µ1 starts in the middle of the bilayer pearling region, the bilayer

pearled fully (end-state not shown) and the value of µ1 increased, see Fig. 8 (bottom-left). The simulations with ξ = −0.25 and ξ = −0.3 begin within the bilayer pearling region, see Fig. 8 (top-left and top-right) and rapidly pearled with the ξ = −0.25 simulation pearling around t = 5 and the ξ = −0.3 simulation fully pearled at the first output time of t = 1. The pearling leads to an increase in µ1 as the bilayer sheds net amphiphilic molecule mass to the bulk (far-field). The end states are depicted in Fig. 7 (left/bottom-left and bottom-right). The bilayer pearling instability was observed by Fraaije and Sevink, who developed a self-consistent mean field density functional model describing the free energy of amphiphilic diblock polymer surfactants embedded in solvent. Their model parameters are based upon poly(propylene oxide)-poly(ethylene oxide) diblock in an aqueous solution, see [2] and reference therein. They simulated spherical nanodroplets of 15% solvent and 85% polymer by volume. By decreasing the block ratio — the ratio of the length of the hydrophilic portion of the diblock chain to the length of the hydrophobic portion, they uncovered a series of bifurcations that lead from stable bilayers, to pearled bilayers to a continuous filament pattern decorated with Y -junctions and endcap defects, see Fig. 7 (center). This change in block ratio from 35% amphiphilic polymer down to 20% amphiphilic polymer increases the aspect ratio of the minority phase. Within the context of the FCH, the parameter η2 weighs the energy requirement of compressing the minority phase into the restricted core of a higher codimensional structure, corresponding to values of u for which the double well W is negative. In particular, bilayer profiles reside in the positive region of the double well, W (φb ) > 0, while the filament profile accesses the negative regions of W , consequently negative values of η2 lower the energy cost of filaments relative to bilayers. In this sense positive values of η2 penalize higher codimensional structures, which is analogous to amphiphilic molecule aspect ratio near one, while negative values of η2 favor the packing found in high codimension morphologies, analogous to a high aspect ratio amphiphilic molecule, see Fig. 7 (right). These preferences are born out in Fig. 6: for fixed values of µ1 , decreasing η2 under constant η1 , hence increasing ηd , results in a crossing of the pearling stability line, leading to the pearled morphology seen in Fig. 7 (center/b and c). Further increase in ηd leads to a defect laden filament structure, see Fig. 7 (center/d) consistent with slight crossing of the filament pearling stability line. This sequence, bilayer stability followed by bilayer pearling and then filament pearling for µ1 > 0 and increasing ηd is consistent with a value of ξ in the range [−0.4, −0.2] depicted in Fig. 6. The qualitative

16

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

Fig. 6. A µ1 –ηd bifurcation diagram versus µ1 (vertical axis) and ηd (horizontal axis) for ξ = −0.85, −0.7, −0.45 and −0.2 and η1 fixed at η1 = 0.15. Pearling instability holds to the right of the dotted lines (blue-bilayers, red-filaments) while regularized curve lengthening holds above the solid lines (blue-bilayers, redfilaments). Areas of pearling instability and curve lengthening are indicated in key regions. The value of µ1 is generically time dependent, and the gray shaded region is forward invariant under the flow . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

agreement between the pearled morphologies observed in the FCH free energy and the self-consistent density-functional model is striking. Both models present emergent pearling with smaller, round holes, compare Fig. 7 (left/bottom-left) and (center/b), while fully emerged pearling leads to larger, pentagonal shaped holes, compare Fig. 7 (left/bottom-right) and (center/c). 4.4. Comparison of analytical bifurcation diagrams to experiments We compare the bifurcation structure derived for codimension 1 and 2 composite morphologies within the FCH gradient flow to experimental results of Dicher and Eisenberg reprinted in Fig. 9 and of Jain and Bates reprinted in Fig. 1. 4.4.1. Bifurcations of Dicher and Eisenberg The bifurcation experiments of Dicher and Eisenberg, [3], depicted in Fig. 1 (left), characterize the end-states of casts of PS–PAA amphiphilic diblocks dispersed in a water–dioxane solvent blend. Dioxane is a good solvent for both PEO and PS, but PS, like its commercial relative styrofoam, is strongly hydrophobic in water. Increasing the water content from zero leads to endstate morphologies, in order: solvability, only micelles, micelles and rods, only rods, rods and vesicles, and only vesicles. Here micelles are codimension 3, while rods and vesicles refer to codimension 2 and 1 respectively. An increase in water content in the solvent is analogous to an increase in η1 — the energy release per unit of interface formation, under constant diblock aspect ratio, hence constant η2 . This morphological bifurcation sequence can be emulated within the FCH equation by fixing ξ =

−0.5 and η2 = −0.4 and allowing η1 to vary from 0 to 0.4. The resulting µ1 –η1 bifurcation diagram is presented in Fig. 9 (left). The curve shortening lines are presented within the view, the pearling stability lines are outside the view. For small values of η1 filaments are dynamically favored, while bilayers are dynamically favored for larger values of η1 . In particular, the horizontal line is color coded to depict a probable end-state morphology of initial data consisting of an admissible composite solution with µ1 = 0.05. Here red denotes a pure filament end-state, blue a pure bilayer end-state and the yellow corresponds to a region of longtime coexistence of the two morphologies due to the approximate equality of the critical values µ∗b = µ∗f for η1 = 0.22. However the bifurcation diagram lies within the pearling instability region of both bilayers and filaments, and as a result the FCH predicts that neither of these pure states would persist. This is a limitation in the FCH model, with the well choice considered here the intersection of the pearling bifurcation curves occurs at smaller values of η1 than the intersection of the curve shortening lines. Agreement with these experimental results requires a robust inhibition of the pearling mechanism. Indeed, the experimental dynamics for the PS–PAA polymers considered here are largely reversible, as shown in Fig. 9 (right), increasing and then decreasing the water content leads to a fully reversible sequence of morphological bifurcations. The presence of the pearling bifurcation generates complex morphologies with strong hysteresis, see Fig. 10 (center). Matching this class of morphological bifurcation diagrams requires a tuning mechanism within the well shape W that affords robust pearling inhibition. Such a mechanism is proposed within the context of multicomponent models in [21] and [22].

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

17

Fig. 7. (left) The t = 100 (τ = 3) end state of simulations of the FCH gradient flow for ε = 0.03, η1 = 0.15, η2 = 0.24 from initial data described in the text and well tilt ξ = −0.15, −0.20 (top row, left to right), ξ = −0.25, −0.3 (bottom row, left to right). The less negative values of ξ have initial values of µ1 that inhibit pearling, for the more negative two values the bilayer pearls. Images courtesy of Andrew Christlieb and Jaylan Jones. (center) Simulations of a mean-field density functional model of amphiphilic diblock copolymers with ratio of amphiphilic component of the diblock decreasing from (a) 35% (b) 30% (c) 25% (d) 20%. The minority solvent-hydrophilic phase is imaged. Reprinted (adapted) with permission from [2]. Copyright 2003 American Chemical Society. (right) Depiction of aspect ratio of lipids, a biological diblock with a short amphiphilic head and a long hydrophobic tail. The aspect ratio is defined as the lipid volume v divided by the head c 2011 WILEY-VCH Verlag GmbH & Co. KGaA, cross-sectional area a and tail length l. From [26], reprinted with permission from the Wiley publication. Copyright ⃝ Weinheim.

4.4.2. Bifurcations of Jain and Bates The experimental bifurcation diagram of Jain and Bates [4], depicted in Fig. 1 (right), shows the end-state morphology of dispersions of PEO-PB with different polymer lengths and different weight fractions of PEO, wPEO , and hence different aspect ratios of the overall amphiphilic diblock. For comparison Fig. 10 (left) depicts the end states of the FCH gradient flow corresponding to values of ξ = −0.2 and η1 = 0.15 for initial values of µ1 = 0.075 and various values of ηd arising from variation in η2 which models the changes in diblock aspect ratio. Small, negative values of η2 correspond to low values of wPEO while larger, positives values of η2 correspond to larger values of wPEO and to negative values of ηd . The horizontal line in Fig. 10 (left) is color coded to depict a probable end-state of the FCH evolution starting from an admissible composite solution with this value of µ1 . On the left where ηd < −0.5, an initial value of µ1 = 0.075, both codimension one and two morphologies are in their curve lengthening region and both would increase in surface area/length. However as µ1 is depleted, the suspension enters the curve shortening region for filaments, which will either vanish in finite time or pearl as µ1 crosses the red dotted line – the end result is a pure bilayer state – indicated by the blue color of the horizontal line for this value of ηd . We focus on the values of ηd in [−0.35, −0.05], for which the horizontal line is colored green, to indicate a region of morphological complexity. The initial data lies in the bilayer pearling region, but passes transiently through it to bilayer pearling stability region. Depending upon the form of the initial bilayer morphology they may either fully pearl and form filament networks, or persist as bilayers and then grow after the return to pearling stability. Filament networks formed from the pearling of a bilayer typically host many Y -junctions and end caps. The filament network will expand until µ1 crosses the red-solid filament curve shortening line, at this point any defect free filaments will shrink, although

the presence of any end-caps and Y -junctions in a particular component will render its evolution unclear. This uncertainty is reflected in the dotted nature of the bottom half of the green vertical line. The end result of the evolution is strongly dependent upon the form of the initial data, and will be very hysteretic in this regime. The uncertainty in the evolution is consistent with co-existing bilayers, filament networks, Y -junctured filaments, and end-cap defects, see Fig. 10 (center) for a depiction of the experimental morphology found in this regime and Fig. 10 (right) for a corresponding end-state of simulation of the FCH gradient flow from random initial data. For large values of ηd the µ1 flow remains in the bilayer pearling instability region, the bilayers either will not form or will pearl and transform to filament networks with the endresult being a Y -juncture dominated filament network. The last prediction, for ηd = −0.02 corresponds to the triple intersection (red circle) of the descending µ1 line, the filament pearling (red-dotted), and curve shortening (red-solid) lines. The arrival of µ1 to the filament curve shortening line from above is consistent with a stable filament phase, but the emergent pearling bifurcation signals a transition to end-cap and micelle formation, and corresponds to the possible coexistence of filament and micelle phases. This transition from filament to filament-micelle is reflected in the dotted nature of the yellow µ1 horizontal line for ηd > −0.02. A final transition to a pure micelle stage is plausible but is outside the scope of our analysis. The overall trend depicted in Fig. 10 (left) suggests increasing ηd at fixed η1 results in bilayers, bilayers mixed with defective filaments (end caps and Y -junctions), filaments, and filaments coexisting with micelles. This bifurcation sequence is in excellent qualitative agreement with the NPB = 170 bifurcation sequence depicted in Fig. 1 (right) as the PEO weight fraction wPEO decreases from high values to low values, corresponding to decreasing η2 and hence increasing ηd subject to constant values of η1 .

18

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

Fig. 8. µ1 –ηd bifurcation diagrams with η1 = 0.15 for the same values of ξ as in the end-states depicted in Fig. 7 (left). For the five simulations the initial and t = 100 final value of µ1 is indicated on the corresponding bifurcation diagram with a solid black circle and square, respectively. The simulations with ξ = −0.25, −0.3 and ξ = −0.2 with ηd = 0 lead to pearled morphologies.

Fig. 9. (left) A µ1 –η1 bifurcation diagram with ξ = −0.5 and η2 = −0.4 that shows the curve shortening lines for bilayers and filaments and vertical arrows showing generic evolution of µ1 from initial data starting at 0.05 for η1 running from 0 to 0.4. The color coding of the µ1 = 0.05 line indicates the final result of the end state with red denoting pure filament, yellow coexistence of filament and bilayer, and blue denoting pure bilayer; compare to Fig. 1 (left) for increasing values of water in solvent phase. (right) Complete reversibility in a PS–PAA system under change in water solvent concentration, implying robust inhibition of pearling instabilities . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Source: From [3]. Reprinted with permission from AAAS.

The bifurcation sequence depicted in Fig. 1 (right) for NPB = 45 corresponds to a much shorter, stiffer diblock polymer. In this regime the network and defect-laden filament phase are not observed, rather increasing wPEO weight fraction leads to the codimensional bifurcation sequence which leads from bilayers, to coexistence of bilayers and filaments, to filaments, to coexistence of filaments and micelles, and finally to micelles. While the general trend of the codimensional bifurcation sequence is supported by the competitive geometric motion and its bifurcations, as in

Fig. 9 (left), we reiterate that within the context of the scalar version of the FCH free energy presented herein, the pearling bifurcation cannot be fully suppressed. 5. Discussion We have presented a multiscale analysis of the H −1 gradient flow of the FCH free energy corresponding to initial data close to

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

19

Fig. 10. (left) A µ1 − ηd bifurcation diagram for ξ = −0.2 and fixed η1 = 0.15, showing the end states of initial data with µ1 = 0.04 and varying values of η2 , with color coding corresponding to probable end state: blue — bilayer, green — network and defect structure (morphological complexity), red — filament, dotted yellow — filament and micelle. Compare to experimental results of Fig. 1 (right) for NPB = 170, with increasing values of wPEO corresponding to decreasing values of η2 . (center) Experimentally observed network, end cap, and Y -junction morphologies corresponding to the CY phase of the bifurcation diagram from Fig. 1 (right). Scale bar is 200 nm, From [4]. Reprinted with permission from AAAS. (right) End state of simulation of FCH gradient flow corresponding to green arrow, coarsened from random initial data, courtesy of Zhengfu Xu . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

dressings of admissible codimension one (bilayer) and codimension two (filament) morphologies. We derived their curvature driven flow which couples to the evolving, spatially constant farfield chemical potential, µ1 . This flow is the basis of the morphological competition, which barring other bifurcations leads to an end state corresponding to the dynamically favored morphology with the lower critical value, µ∗b or µ∗f , of the far field chemical potential. In particular we identify regimes in which the geometric flow leads to growth or evanescence of each phase through curve shortening or regularized curve lengthening. Combining the curve shortening/regularized curve lengthening bifurcation with the pearling bifurcation results allows a characterization of the evolution of defect-free connected components of bilayer and filament morphologies. Our analysis predicts that codimension one and two structures do not generically coexist on the long, t = O(ε −1 ) time scale within the strong FCH gradient flow. Experimental results show transitions between distinct codimensional phases with relatively large margins of coexistence, however as remarked in the experimental literature, the transients associated to these experiments are long compared to experimental patience and transients required to achieve a single phase of codimensional morphology may require months to years, [5]. We compared the analytical bifurcation diagram to results from numerical simulations of the FCH gradient flow, to simulations of a self-consistent mean-field density functional model for amphiphilic polymers, and to three sets of experimental studies of amphiphilic polymers. We find that the self-consistent mean-field density functional model proposed in [2] predicts a bifurcation sequence of radial bilayers, pearled bilayers, and filaments with end-cap defects in strong qualitative agreement with the FCH bifurcation structure, in particular the structure of the pearled spherical bilayers computed by both models are in excellent agreement. The experimental bifurcation analysis of Jain and Bates, [4], was conducted at two polymer lengths, long polymers with NPB = 170 and shorter ones with NPB = 45. We find strong qualitative agreement between the FCH bifurcation structure and the experimental results for the longer chains, with the FCH results suggesting that the development of morphological complexity could be produced in a region in which bilayer pearling and filament curve shortening bifurcations lie in close proximity, see Section 4.4.2. The passage of the far-field chemical potential through these bifurcations engenders network morphologies with Y -junctions, end-caps, and stable pearled filaments. The analytical basis of the complexity suggests a mechanism for hysteresis:

the passage through these sequences of bifurcations will not be readily reversed by a non-adiabatic return of the state variables. The FCH model with the parameter choices presented herein does not qualitatively reproduce all of the experimental results. Both the solvent quality bifurcation experiments of Dicher and Eisenberg, [3] and the short polymer NPB = 45 experiments of Jain and Bates do not show evidence of pearling bifurcations. Complex morphology is not exhibited, and Dicher and Eisenberg show that for their experimental parameters the morphology is remarkably reversible: hysteresis is not observed. However, the Dicher and Eisenberg experimental bifurcation structure is well described by the dynamical favoritism arising from the morphological competition. To match the full range of experimental results the FCH free energy needs a tunable mechanism to robustly inhibit the pearling bifurcation. These can be achieved in two ways. The first is physically motivated: pearling is a modulation of bilayer width, and shorter polymers are stiffer and lead to bilayers with a less compressible width. The compressibility of the bilayer is tunable through the slope of the well W at the value u at which the bilayer density is greatest, generically the second W = 0 crossing. Tuning this slope to be large represents a stiffer diblock and may serve to inhibit the pearling mechanism. The second mechanism is mathematically motivated: the pearling bifurcation arises from a balance between the positive eigenvalue λb,0 of the linearization Lb,0 , see (2.7) about the bilayer profile, and the negative eigenspace of the Laplace–Beltrami operator, ∆s . The balance cannot occur if the operator Lb,0 is non-self adjoint and whose positive real part spectrum have non-zero imaginary parts. This arrangement can be tuned and detuned within a multicomponent model, providing precisely the desired mechanism for robust pearling inhibition. This mechanism was discussed in section 5 of [21] and forms the basis of the study of the singularly perturbed systems in [22]. There remains a considerable amount to address within the family of FCH energies. While the FCH is at some level a phenomenological model with parameter values that are not directly tunable from first principles, the FCH parameter values can be fit to experiments, or subscale molecular simulations, much like Flory–Huggins parameters. Nonetheless, fitting the full form of a complicated well W , especially for multicomponent models, may be challenging. Even within the simple model presented here, the role of micelles on morphological competition has not been addressed, and their stability, including their growth into dumbbells and end-caped filaments is a primary instability mode

20

A. Christlieb, N. Kraitzman and K. Promislow / Physica D 400 (2019) 132144

that leads adiabatically from a codimension three structure into a defective codimension two structures. As there are almost no analytical characterization of defect modes, the micelle to dumbbell instability is worthy of study on a purely mathematical basis. Acknowledgment The third author acknowledges support by the National Science Foundation through grants DMS-1409940 and DMS-1813203. All authors thank Zhengfu Xu for providing the numerical simulation which generated Fig. 10 (right). References [1] A. Blanazs, S.P. Armes, A.J. Ryan, Self-assembled block copolymer aggregates: from micelles to vesicles and their biological applications, Macromol. Rapid Commun. 30 (4–5) (2009) 267–277. [2] J. Fraaije, G. Sevink, Model for pattern formation in polymer surfactant nanodroplets, Macromolecules 36 (21) (2003) 7891–7893. [3] D.E. Discher, A. Eisenberg, Polymer vesicles, Science 297 (5583) (2002) 967–973. [4] S. Jain, F.S. Bates, On the origins of morphological complexity in block copolymer surfactants, Science 300 (5618) (2003) 460–464. [5] S. Jain, F.S. Bates, Consequences of nonergodicity in aqueous binary peo-pb micellar dispersions, Macromolecules 37 (4) (2004) 1511–1523. [6] N. Kraitzman, K. Promislow, Pearling bifurcations in the strong functionalized Cahn–Hilliard free energy, SIAM J. Math. Anal. 50 (3) (2018) 3395–3426. [7] K. Promislow, B. Wetton, PEM fuel cells: a mathematical overview, SIAM J. Appl. Math. 70 (2) (2009) 369–409. [8] N. Gavish, G. Hayrapetyan, K. Promislow, L. Yang, Curvature driven flow of bi-layer interfaces, Physica D 240 (7) (2011) 675–693. [9] S. Dai, K. Promislow, Geometric evolution of bilayers under the functionalized Cahn–Hilliard equation, in: Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, Vol. 469, The Royal Society, 2013, p. 20120505. [10] A. Doelman, G. Hayrapetyan, K. Promislow, B. Wetton, Meander and pearling of single-curvature bilayer interfaces in the functionalized Cahn–Hilliard equation, SIAM J. Math. Anal. 46 (6) (2014) 3640–3677.

[11] K. Promislow, Q. Wu, Existence of pearled patterns in the planar functionalized cahn–hilliard equation, J. Differential Equations 259 (7) (2015) 3298–3343. [12] D.A. Scherlis, J.-L. Fattebert, F. Gygi, M. Cococcioni, N. Marzari, A unified electrostatic and cavitation model for first-principles molecular dynamics in solution, J. Chem. Phys. 124 (7) (2006) 074103. [13] O. Andreussi, I. Dabo, N. Marzari, Revised self-consistent continuum solvation in electronic-structure calculations, J. Chem. Phys. 136 (6) (2012) 064102. [14] J. Zhu, N. Ferrer, R.C. Hayward, Tuning the assembly of amphiphilic block copolymers through instabilities of solvent/water interfaces in the presence of aqueous surfactants, Soft Matter 5 (12) (2009) 2471–2478. [15] J. Zhu, R.C. Hayward, Interfacial tension of evaporating emulsion droplets containing amphiphilic block copolymers: Effects of solvent and polymer composition, J. Colloid Interface Sci. 365 (1) (2012) 275–279. [16] S. Dai, K. Promislow, Competitive geometric evolution of amphiphilic interfaces, SIAM J. Math. Anal. 47 (1) (2015) 347–380. [17] G. Hayrapetyan, K. Promislow, Spectra of functionalized operators arising from hypersurfaces, Z. Angew. Math. Phys. 66 (3) (2015) 631–662. [18] I.R. Cooke, M. Deserno, Coupling between lipid shape and membrane curvature, Biophys. J. 91 (2) (2006) 487–495. [19] S. Dai, K. Promislow, Geometrically constrained minimizers in dispersions of strongly hydrophobic amphiphilic polymer, Preprint. [20] Y. Chen, K. Promislow, Regularized curve lengthening for the strong FCH flow, Preprint. [21] K. Promislow, Q. Wu, Existence, bifurcation, and geometric evolution of quasi-bilayers in the multicomponent functionalized Cahn–Hilliard equation, J. Math. Biol. 75 (2) (2017) 443–489. [22] A. Doelman, K. Promislow, F. Veerman, Freeways, tollroads, and pearling inhibition in multicomponent functionalized systems, Preprint. [23] R.L. Pego, Front migration in the nonlinear Cahn-Hilliard equation, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 422 (1863) (1989) 261–278. [24] N. Alikakos, P. Bates, X. Chen, Periodic traveling waves and locating oscillating patterns in multidimensional domains, Trans. Amer. Math. Soc. 351 (7) (1999) 2777–2805. [25] N.D. Alikakos, P.W. Bates, X. Chen, G. Fusco, Mullins-sekerka motion of small droplets on a fixed boundary, J. Geom. Anal. 10 (4) (2000) 575–596. [26] O.G. Mouritsen, Lipids, curvature, and nano-medicine, Eur. J. Lipid Sci. Technol. 113 (10) (2011) 1174–1187.