Variational Models of Network Formation and Ion Transport: Applications to Perfluorosulfonate Ionomer Membranes

We present the functionalized Cahn-Hilliard (FCH) energy, a continuum characterization of interfacial energy whose minimizers describe the network morphology of solvated functionalized polymer membranes. With a small set of parameters the FCH characterizes bilayer, pore-like, and micelle network structures. The gradient flows derived from the FCH describe the interactions between these structures, including the merging and pinch-off of endcaps and formation of junctions central to the generation of network morphologies. We couple the FCH gradient flow to a model of ionic transport which incorporates entropic effects to localize counter-ions, yielding a flow which dissipates a total free energy, and an expression for the excess electrochemical potential which combines electrostatic and entropic effects. We present applications to network bifurcation and membrane casting.


Introduction
A central goal of polymer chemistry is the design of materials with novel macroscopic properties by controlling the spontaneous generation of nanoscaled, phase separated networks [1,2]. A primary mechanism to generate such networks is through the "functionalization" of hydrophobic polymer chains by the addition of acid terminated side-chains, such as in the casting of perfluorosulfonic acid (PFSA) membranes, see Figure 1. The acid-tipped end-groups interact exothermically with solvents, driving the generation of charged polymer-solvent interface. By virtue of their tethered charges, these networks are selectively conductive, inhibiting the transport of co-ions while facilitating the transport of counter-ions. The broad class of ionomer membranes have important applications in efficient energy conversion devices such as polymer electrolyte membranes for fuel cells [3][4][5][6], dye sensitized solar cells [7], bulk-heterojunction solar cells [8,9], and lithium ion batteries [1,10,11]. Considerable computational and experimental effort has been invested in the design and optimization of membranes with strongly charge-selective transport characteristics. We present a continuum approach which affords a surprisingly compact description of the network morphology in these complex systems. Perhaps the most successful example of a reduction of a complex system is the hydrodynamic limit of an over-damped, highly collisional gas. By identifying the equilibrium modes of the collision operator, namely the Maxwellian velocity distributions, and projecting the residual convective terms onto the tangent plane of the Maxwellians, the associated Boltzmann equation can be rigorously reduced, in various limits, to the Euler or the Navier Stokes equation, see [12] or the review article [13]. The relaxation of ionomer-solvent mixtures is also heavily over-damped, but the underlying physics is considerably more complex than the simple conservation of momentum and energy embodied in collision operators. Nevertheless, we propose, without rigor, that the richly structured composition of solvent and ionic groups that forms within interfaces, if properly characterized, can play a role similar to the Maxwellians of a collision operator. That is, the polymer-solvent configurations form a rich family of local quasi-equilbrium, or adiabatic, structures which parameterize the configurational space occupied by the full system as it evolves under a long-term flow. Moreover, the residual that drives the slow flow, the convective terms in the analogy to the Boltzmann equation, can be characterized by the entropy of the interaction of solvent with the tethered ionic groups and the counter ions, which are protons for PFSAs.
We incorporate these effects into a general framework, the Functionalized Cahn-Hilliard (FCH) energy, for an interfacial energy landscape whose over-damped (gradient) flows recover much of the morphological complexity observed in the polymeric membranes. We do not propose a mechanism to determine the explicit form of the constitutive relations from the underlying structure of the chemical system, but rather view the parameters of the FCH energy as targets of upscaling from particle based simulations and experimental data, in the spirit of the viscosity and other macroscopic parameters of a hydrodynamic flow [14]. The FCH is a possible upscaling target for data on preferred states of ionic solvation which can be integrated into a continuum framework for the evolution of network morphology, a framework whose computational cost is many orders of magnitude lower than that of the underlying particle system.
Even at a course-grained level, molecular and dissipative particle dynamics simulations of solvated, hydrophilic-hydrophobic polymer systems cannot compute the morphological development to its equilibrium [15], even in state-of-the-art simulations resolved to milliseconds of real time [16]. Indeed experimental evidence suggests PFSAs can possess 10-20 h transients after changes in external environment [17]. Continuum models of network formation in Nafion [18][19][20], have successfully predicted geometric features, such as pore radius and overall solvent uptake as a function of ionic group density, liquid phase pressure, and effective elastic response of the polymer phase. However they do not predict the structure of the over-all solvent network, nor can they describe the geometric evolution of the network, which they take as a model assumption.
Semi-crystalline PFSAs are widely recognized as the benchmark materials in polymer electrolyte membrane fuel cell applications [21,22]. While these polymers inherently possess outstanding chemical and physical properties, the harsh conditions encountered in operating fuel cells and the current demands for long-term durability have revealed that these benchmark materials are susceptible to degradation, particularly in the cast form required to make thin membranes currently in use, as opposed to the thicker extruded version of Nafion. Consequently, the fuel cell community has shifted attention to the importance of the polymer morphology in optimizing macroscopic properties.
The casting of PFSAs involves their dissolution in solvent, and then evaporation of the solvent under various heat treatments. It is thought that ionic groups are not uniformly distributed along the backbone, and as solvent evaporates during the casting, regions of low ionic densities may lead the backbone groups to crystallize either into isolated clusters or into a percolating matrix with length scales on the order of 10 to 100 nm [21]. SAXS data suggests that in regions of high ionic density the ionic groups cluster forming a network of 1-5 nm diameter pores filled with solvent and lined with the ionic groups. Key fuel cell membrane characteristics such as water uptake, gas transport, proton conductivity, and mechanical stability are linked to the presence and organization of the ionic and crystalline domains in PFSAs. Future improvement in PFSA performance requires a fundamental understanding of the formation and function of these essential morphological features.
While the the continuum approach we present has a wide range of potential applications, including the development of reduced dimensional models of the morphology of amphiphilic diblock copolymers, vesicle fusion of cellular bilayers, and a broad range of issues in cellular morphogenesis, we restrict our attention here to PFSA membranes. This paper is organized as follows. In Section 2 we consider binary mixtures and discuss an interpretation of the model in the context of PFSAs. In Section 3 we present applications of the binary model to replication of SAXS data for well-hydrated Nafion, and discuss extensions of the energy to include ion transport within the solvent phase. In Section 4 we present extensions to multi-component models and a brief discussion of their application to the casting of PFSAs.

Network Formation in Binary Mixtures
A widely cited work of Hsu and Gierke [23], describes a cluster-network morphology of hydrated Nafion which they ascribe to a balance of elastic energy, from deformation of stiff backbone material, against hydrophilic surface interactions involving solvent molecules and the ionic end-groups. We embed these ideas within a novel mathematical framework which naturally incorporates a selection mechanism between competing morphologies. We first consider a binary mixture, viewing well-hydrated PFSAs as a hydrophobic polymer network and a solvent phase, with a fixed charge density of ionic groups at the interface. We describe the morphology in terms of the dimensionless volume fraction, u ∈ [0, 1], of the solvent phase relative to an irreducible water volume fraction. The fundamental assumption is that the network structures observed in polymer electrolyte membranes can be well-described by solutions of a Cahn-Hilliard (also known as Ginzburgh-Landau) critical-point equation where u is a function of the material point x within the spatial domain Ω ⊂ R 3 . The function W b is a double-well potential for the binary mixture that ascribes energy densities to the various mixture values. It generically has unequal-minima at the pure-state volume fractions u = 0, 1. The parameter > 0 is a dimensionless ratio of length scales, and ∆ denotes the Laplacian operator. The rich family of solutions of Equation (1) are the critical points (maxima, minima, or saddles) of the Cahn Hilliard (or Ginzburgh-Landau) interfacial free energy The mass constrained minimizers of the Cahn-Hilliard energy have a distinguished history [24]. It is well understood that, independent of the values of the self-energies of the two pure states (the values of W b at its two minima) the global minimizers of E CH are comprised of spatial domains of approximately pure states, u ≈ 0 or 1, arranged so that the interface between the domains has minimal surface area. Moreover u varies monotonically from 0 to 1 in the direction transverse to the interface, whose width scales with . These structures are called single-layer interfaces. In particular for the minimizer, u, as → 0 + , the morphology of the spatial domains where u ≈ 0 and u ≈ 1, does not change, the width of the interface between them merely shrinks to zero.
However, for the applications we consider it is not the minimizers of the Cahn-Hilliard but the vast families of saddle points that are physically relevant. Indeed, it is considerably less well-known that the critical points of the Cahn-Hilliard energy encompass large classes of network structures: domains for which the characteristic width of the minority phase domain scales with , giving rise to long, thin, percolating structures: bilayers, pores, pearled-pores, and micelle clusters [25]. The network structures are fundamentally distinct from the more familiar surface area minimizers, as → 0 + , the network morphologies grow thinner, and longer, and do not approach a fixed limit. It is these network structures that play the role of "Maxwellians" within our energetic framework.
Our approach rests upon a reformulation of the Cahn-Hilliard energy, generating a new free energy which possesses the same families of critical points as the Cahn-Hilliard but which permits a facile mechanism to select the critical points with desired attributes as its minimizers. Most significantly we show that the selection mechanisms have physically meaningful interpretations; they provide the parameters that serve as macroscopic characterizations for molecular-scale properties. The starting point for the reformulation of the Cahn-Hilliard energy is its variational derivative at a configuration u, whose zeros are the solutions of the critical point Equation (1). The square of the variational derivative, when integrated over the spatial domain Ω, gives an energy which possesses the same critical points as the Cahn-Hilliard energy, but for which all critical points are self-evidently global minimizers with minimum energy zero. To effect the balance described by Hsu and Gierke, we subtract small terms which describe the interfacial area and relative entropy of the solvent-counter ion mixture. Since this new energy incorporates the effect of the functionalization of the underlying hydrophobic polymer backbones, and to a lesser extend because of its origins in the functional calculus [26], we refer to it as the Functionalized Cahn-Hilliard free energy An interpretation of the parameters is provided in Section 2.2, however for bilayer and pore networks the first term of the FCH represents elastic energy, as measured by the square of an interface's curvature while the negative terms parameterized by η h and η m represent the hydrophilic interactions associated with solvation of ionic groups and mixture entropy of the counter ions and solvent. It has been shown that the FCH energy is bounded below and has global minimizers over natural spaces of functions [27].
For positive values of η m and η h the competitors for the minima of the FCH energy are critical points of the Cahn-Hilliard energy whose surface area is maximal in an appropriate measure.
It is natural to compare the FCH energy to a sharp-interface model for bilayer structures, whose energy depends only upon the curvatures κ 1 and κ 2 of the two-dimensional interface Γ embedded in R 3 . The quintessential sharp-interface energy is the Canham-Helfrich energy [28,29], which considers the most general symmetric, quadratic polynomial in the two curvatures. In terms of the mean, H = 1 2 (κ 1 + κ 2 ), and Gaussian curvatures, K = κ 1 κ 2 , the energy takes the form The parameters a 1 and a 4 , denote the energy density attributed to the curvatures, with a 2 specifies the intrinsic, or zero-energy, value of the mean curvature H. The parameter a 3 denotes the energy per unit of interfacial area. Since the integral of the Gaussian curvature, K, over a closed interface depends only upon the interface's Euler characteristic, the Gaussian term is usually omitted, see [30] for further discussion. In Section 2.3 we show that the sharp interface limit of the FCH energy at a bilayer reduces to a Canham-Helfrich energy with a 2 = a 4 = 0 and a 3 = a 3 (η m , η h ). Indeed we choose the Cahn-Hilliard energy as the starting point for our energy precisely because it is a finite-width regularization of surface area whose variational derivative reduces to mean curvature. The functionalization of any energy with these two properties would yield an equivalent formulation. However we emphasize that the FCH energy encompasses a wide-range of interfaces beyond bilayers, including pore and micelle networks, as well as Y-junctions and bilayer-pore hybrids. It also permits, but penalizes defect structures, such as end-caps. Indeed, depending upon parameter values, the FCH energy prefers to either merge the end-cap into a network structure to form a lower-energy pore network, to pinch them off into isolated micelles, or to evolve into pearled pore structures, see  (6) with periodic boundary conditions. Each simulation starts from an identical initial data, which was chosen to be randomly 0 or 1 at each grid point with a 20% solvent volume fraction. The difference between the simulations lies in the value of the mixture parameter, η m . The binary well W b has tilt τ = −0.4, = 0.03, E b = 1, and η h = 5 , yield (left to right) bilayer-pore mixture, a highly branched pore network, a pore network, disconnected pores, a pearled pore and micelle mixture, and isolated micelles.  Initial Data t = 640 t = 10, 000

Gradient Flows of the FCH Energy
We are interested not only in equilibrium configurations, but also in the influence of material history on morphology. Since the slow evolution of morphology is heavily over-damped in comparison to the time scales required to establish interfacial structures, it is appropriate to model the time evolution as a mass-preserving gradient flow of the FCH energy. A significant advantage to a phase-field approach, such as the FCH energy, over a sharp interface energy, such as the Canham-Hilfrich energy, is that its gradient flows automatically incorporate merging and pinch-off events, as well as higher co-dimension and defect structures. Moreover, as we detail below, the FCH energy also couples naturally to other physical processes, such as mass and momentum balances [31]. The mass-preserving gradient flow on the FCH energy takes the form, In the simulations considered here, we supplement Equation (6) with periodic boundary conditions on a cube Ω ⊂ R 3 . The resulting flow decreases the FCH energy along orbits, while preserving the total mass fraction of each component [32]. Other boundary conditions are also natural to consider, in particular the normal component of the flux of water phase can be prescribed at a membrane-air interface. However modeling the air-interface boundary of PFSA membranes involves non-trivial physical considerations, including capillary pressure jumps, wettability of liquid phases in contact with the membrane, and non-uniform swelling of the membrane. Within this article we consider only the bulk phase of PFSA membranes, in the absence of external fluxes, which is consistent with periodic boundary conditions. To emphasize the impact of the parameter η m on morphology, Figure 2 presents the end-state of six simulations of Equation (6) from identical initial data (random 0 or 1) with solvent volume fraction 20%. By varying only the value of the parameter, η m , quasi-equilibria are obtained which are a mixture of bilayer and pore solutions, pore networks with differing levels of branching, pearled pore and micelle mixtures, and isolated micelles. To show merging and network formation, Figure 3 presents the time evolution of the gradient flow from random initial data for parameters in the transition from pore to micelles. Micelles initially form by agglomerating the random initial distribution of solvent, however after the initial transient the micelles either merge or become unstable and grow into dumb-bell like structures, essentially a short pore with two end-caps. The pore-stems elongate and where end-caps encounter another pore stem, they merge into a Y junction, while the intersection of two end-caps results in a merging into a longer pore. Eventually a pore-network structure forms with several remaining end-caps. Depending upon parameter values, the pores may absorb the end-caps, grow longer and meander, or as depicted in Figure 3, relax to a local equilbria with a mixture of end-cap structures and micelles. In Figure 4, we present a time evolution starting from a low-surface area initial data, a solid sphere of solvent corresponding to a macroscopic mixture of water and polymer. For η h > 0 the large solid sphere is unstable and grows end-cap structures which elongate and evolve into a high-surface area pore network with a mixture of Y junctions and end-caps.

Interpretation of Energetic Terms and Parameters
The competition for minimization of the FCH, and the flow along families of quasi-minimizers, balances four competing effects. First the phases should separate into interfacial structures, the critical points of the Cahn-Hilliard energy, parameterized primarily by the well-tilt parameter τ b . Curvature of the interfaces, parameterized by the effective bulk modulus E b , should be minimized due to polymer backbone stretching. Surface area, as parameterized by η h , should be maximized to separate charged groups and enhance access to solvent. Protonic groups within the solvent phase maximize their entropy relative to available hydrogen bond densities, as characterized by η s and W s .
The dominant term in the FCH energy is the square of the variational derivative of the Cahn-Hilliard energy, all other terms are small relative to the coefficient E b , which has units of energy per volume. A mixture distribution u which minimizes the FCH energy must have a small residual with respect to the Cahn-Hilliard critical point Equation (1). The principle measurement of smallness is the length, ξ int , over which variation in interfacial structure is observed. In the context of PFSAs, ξ int corresponds to the length of the functional side-chains plus the solvation shells of the charged groups, typically a fraction of a nanometer. The parameter is the ratio of ξ int to the width of the fixed physical domain, Ω. In applications with acid group densities on the order of 1 molar, the Debye length is typically many multiples of ξ int . For computational and organizational simplicity we fix the spatial domain, Ω, and take ≈ 1 − 5 × 10 −2 1, to fit a large number of interfacial structures within the spatial domain. The parameters η h and η m are also small, and to study structures arising from balances between terms it is natural to scale them with . In applications to functionalized polymers in good solvent, η h is generically positive while η m which corresponds to a relative entropy, and plays a key role in morphology selection, can be either positive or negative.
The most important element of the binary double-well potential W b is the difference in self-energy of the pure states: We consider binary wells of the form The choice Equation (8) assures that the minima remain at 0 and 1 independent of the well-tilt parameter τ . Moreover for small τ , this is the optimal well tilt as measured by the bifurcation rate of double-layers from a single-layer interface [32]. We emphasize that the value u = 0 corresponds not to "pure polymer" but to a polymer blend with an irreducible amount of solvent, for many PFSA membranes this is widely viewed as two-three water molecules per acid group.
The square of the variational derivative. The critical point equation is singularly-perturbed, possessing boundary layers (fast or inner structures) separating equilibria characterized by the pure states u = 0, 1 (outer or slow structures). Unlike the minimizers of the Cahn-Hilliard energy, whose unique inner structure is a single-layer interface, the potential minimizers of the FCH free energy can possess three dominant boundary-layer structures characterized by their distinct numbers of fast and slow directions. We refer to the number of fast directions as the co-dimension, since the slow directions are naturally viewed as tangent to the center-line or center-plane of the interface, while the fast directions (or direction) are measured along the normals to the interface. Co-dimension one solutions correspond to membranes or bi-layer structures in R 3 , with a center-plane given by a two-dimensional surface Γ 1 with a single normal (fast) direction. The in-plane coordinates are denoted s = (s 1 , s 2 ) and through-plane distance, scaled by , is denoted z. Local to each fixed interface Γ 1 there is a standard change of variables from cartesian to (z, s) coordinates, in which the Laplacian transforms to where ∆ s is the Laplace-Beltrami operator, which represents surface diffusion on Γ 1 , and κ(z, s) = H(s) + zκ 1 (s) + O( 2 ) has leading order term H(s), the mean curvature of Γ 1 at s. At leading order, co-dimension one solutions of Equation (1) are described by an (arbitrary) closed two-dimensional interface Γ 1 and the solution φ 1 (z) of fast component of the critical point equation, If the potential W b has unequal depth wells, τ > 0, then this equation possesses a unique homoclinic solution which approaches 0 as z → ±∞, corresponding to a polymer backbround φ 1 = 0 on either side of a thin solvent bilayer, see Figure 5 (right). In three space dimensions the coordinate z can be viewed of a function of Cartesian position x (and Γ 1 ), that is z = z(x; Γ 1 ), see Figure 6, so that φ 1 induces a solution u 1 (x; Γ 1 ) := φ 1 (z(x; Γ 1 )). For this solution the center of the bilayer, z = 0, corresponds to the interface Γ 1 , see Figure 6 (left). If Γ 1 has non-zero mean curvature, then the dominant residual in the square of the variational derivative is Since φ 1 decays to zero at O( ) distances from the interface, the mean-curvature squared term is localized on the interface. In this sense, for co-dimension one bi-layer structures, the square of the variational reduces to the square of the mean curvature of the center-plane-which has the interpretation as the deformational or elastic energy of the thin bi-layer structure. α =3 Figure 6. (Left) The co-dimension one coordinate system has a fast normal z and two slow directions s = (s 1 , s 2 ) lying in the tangent plane of Γ 1 ; (Right) The co-dimension two coordinate system has a two dimensional fast sub-system parameterized by cylindrical coordinates (R, θ) where R is the distance from the pore center and θ is the polar angle, and a one dimensional slow system parameterized by distance s along Γ 2 . s s 1 Co-dimension two solutions are pore or spaghetti-like: thin, solid cylinders comprised of one phase surrounded by the complementary phase. A distinction should be drawn between a co-dimension two solution and a Lysosome structure formed by wrapping a bi-layer (co-dimension one) into a cylindrical (hollow) tube. Co-dimension two structures are described by a one-dimensional curve Γ 2 and a local coordinate system comprised of two fast variables, naturally described in -scaled (fast) polar coordinates (R, θ), and a single tangential (slow) coordinate s which records distance along Γ, see Figure 6 (right).
In these variables the cartesian Laplacian takes the form where K 1 and K 2 depend upon the curvature and torsion κ = (κ 1 (s), κ 2 (s)) of Γ 2 and on the angle θ, while D 2 s = ∂ 2 s + O( R| κ|) is at leading order the second derivative with respect to arc-length along Γ 2 . We consider only inner solutions which are independent of θ, building leading-order, co-dimension two, critical points from the cylindrically symmetric solution of which satisfies ∂ R φ 2 (0) = 0 and the leading-order far-field condition φ 2 (∞) = 0. Again, given a one-dimensional curve Γ 2 , these induce cartesian solutions u 2 (x) := φ 2 (R(x; Γ 2 )). For this function u 2 = u 2 (x; Γ 2 ), the residual of the squared variational derivative term manifests itself as the square of the curvature vector of Γ 2 .
In R 3 , the co-dimension three variables reduce to the spherical coordinate system about center point x 0 . Spherically symmetric critical points satisfy subject to ∂ R φ 3 (0) = 0 and the leading-order far-field equation φ 3 (R = ∞) = 0. Local to x 0 we have the relation R = R(x) from which the critical point φ 3 induces the cartesian solution The well-tilt parameter: τ . The radius of solvent filled pores within the charged polymer has been described at a continuum level [20,33], as a balance between the solvent pressure against the elastic pressure exerted upon the solvent by the polymer phase. The solvent pressure, P , depends strongly upon the density of ionic groups on the surface of the polymer which determines the corresponding density of counter ions immersed within the solvent. The solvent pressure also depends upon capillary effects at the external boundary of the material. The elastic pressure of the cross-linked polymer phase arises from its compression by the swelling solvent phase. Literature values [20,23] for the effective elastic modulus of the polymer phase range from 100-200 J/cm 3 . The osmotic pressure within the liquid phase takes the form P osm = RT c p where R denotes the gas constant, and c p is the proton density. Assuming an interfacial charge of 0.5 C/m 2 [20], and a 5 nm radius yields an osmotic pressure of 80-100 J/cm 3 . Local variations in ionic group density as well as changes in liquid pressure due to external forces or changes in surface capillary pressure would impact this balance. Within the FCH framework, the well-tilt parameter, τ , encodes this equilibrium balance, determining the radius of bilayer, pore, and micelle structures, in units of , as depicted in Figure 5. In a general framework, τ will exhibit a functional dependence upon the total solvent pressure, P = P osm + P l , which incorporates the osmotic pressure of the ions through the local ion density, n, and the bulk elastic moduli, E b of the polymer, τ = τ (P, E b ).
The hydrophilic energy: η h . The hydrophilic energy term η h 2 2 |∇u| 2 is localized on the polymer-solvent interface, where the gradient of u is non-zero, and represents the free energy released by the formation of the polymer-solvent surface area. Other effects held constant, an increase in surface area facilitates the access of the charged groups to the solvent, enabling them to build up a solvation shell. Moreover, increases in surface area also permit the charged groups to spread out, particularly charged groups which are proximal but are tethered to distinct strands of the polymer backbone. Literature values of the Gibbs free energy of water sorption run from 25-50 kJ/mol [20,34]. At ionic group densities typical of Nafion, this value is 15-30 times higher than the energy density associated to the backbone stiffness. The magnitude of the hydrophilic energy parameter depends upon the local density of acid groups, η h = η h (n), as well as the nature of the solvent. Larger values of η h favor lower co-dimensional structures which produce more surface area per volume of solvent, see Figure 5.
The mixture entropy: η m . The critical point Equation (1) determines the profile of solvent volume fraction at the solvent-polymer interfaces, particularly the purity of the solvent phase at the structure's center (mid-point of a bi-layer, center-line of a pore, or center of a micelle), see Figure 5. This value is sensitive to the well-tilt parameter τ , with the higher co-dimension structures having the greater purity of solvent. For PFSA's, in which the counter ion is a proton and the solvent is water, the state of the solvent molecules, as reflected by the associated dielectric constant, is thought to vary significantly from a low value for water molecules solvating the acid groups, to a high value for bulk-like water, even in the proximity of the protonic groups [35], see Figure 7. Indeed protons evidence a strong preference to delocalize within a bulk-like water regime. The mixture entropy parameter, η m , encodes the difference in free energies, per unit volume, for the various states of the solvent-proton mixture, particularly the strength of the entropic preference for the protons to reside in bulk-like water over the low-dielectric water. For protonic counter ions, the solvation entropy term, W s , should be low for bulk-like water and high for water of solvation, that is W s (0) > W s (1). For simplicity we set W s = W b , since W s (0) = 0 the value of η m controls the sign and magnitude of the mixing entropy. Specifically, for η m < 0 the mixture term gives energetic preference to bulk-type solvent and hence favors higher co-dimensional structures, while for η m > 0 the mixture term gives relative preference to co-dimension 1 bilayers which have less bulk water than pores or micelles. The entropic parameter has a functional dependence upon the density of counter-ions, however due to electroneutrality over length scales large compared to , it is reasonable to take η m = η m (n). Figure 7. Statistical mechanical models of the dielectric ε and hydronium concentration, c H + in a pore at three water contents [36,37]. The dielectric ranges from 2 to 80 while the proton density decays to zero at pore walls. Reprinted (adapted) with permission from Chem. Rev. 2004, 104,

4637-4678. Copyright (2012) American Chemical Society
The effective bulk elastic energy, 2 E b . As discussed above, the physical values of the bulk elastic energy are typically 25-50 times smaller than the hydrophilic interactions associated to η h and η m . However, because of the form of the second variation, the coefficient E b obtains a natural scaling by 2 , see Equation (20). In our simulations we have taken values for η h and η m on the order of , while we choose E b = 1. For ≈ 0.03, these choices yield effective values of bulk elastic energy, 2 E b , and hydrophilic interaction, η h with a ratio of 2 E b /η h ≈ ≈ 0.03, in harmony with literature values.

Sharp-Interface Reductions of the FCH Energy
An intuition for the nature of the FCH energy can be enhance by examining its sharp-interface reductions, the → 0 + limit, for fixed co-dimension one, two, or three structures without end-caps or junctions. A sharp-interface energy of an interface depends only upon the interface's size and curvatures in its slow coordinates. A co-dimension one quasi-equilibria of the FCH is parameterized by its two-dimensional center-line, Γ 1 , which we take to be without boundary (end-caps) or self-intersection (junctions). Recalling the Cartesian form of the co-dimension one critical point u 1 (x) = φ 1 (z(x); Γ 1 , τ ), we evaluate F at u 1 . In a neighborhood of Γ 1 we change to the (z, s) variables in the integral defining F. Keeping only leading order contributions (in ) we obtain the sharp-interface (leading-order) reduction F (u 1 (x; Γ 1 , τ )) = σ 1 (τ ) where σ 1 (τ ) = 1 2 ∞ −∞ (φ 1 (z; τ )) 2 dz. This is a precise mathematical translation of Hsu and Gierke's, elastic and hydrophilic balance: the square of the mean curvature, H, of the interface balances against the hydrophilic interactions characterized by the solvation and mixture parameters η h and η m . So long as the integrand is negative, the free energy can be decreased by increasing the surface area, |Γ 1 | of the bilayer structures.
A co-dimension two quasi-equilibria is determined by its one-dimensional center-line, Γ 2 , which generates a local (R(x), θ(x), s(x)) coordinate system and induces the critical point u 2 (x) = φ 2 (R(x); Γ 2 , τ ). Evaluating F at u 2 and changing variables in the integral yields the sharp-interface free energy where σ 2 = 1 2 ∞ 0 (φ 2 (R; τ )) 2 RdR and κ(s) is the vector curvature of Γ 2 at position s. Finally, co-dimension 3 quasi-equlibria are determined by the locations Γ 3 = { x 1 , . . . , x N } of the N micelles. Their contribution to the free energy takes the reduced form where Naturally enough the co-dimension three micelles have no curvature contribution to their energy as they have no spatially extended directions.
A key step in the reduction to the sharp interface is the evaluation of the volume integrals of the mixture entropy term. For the particular choice W s = W b of the mixture entropy potential, the volume integrals have simple closed forms In the reduced energy formulation, these integrals lead to the coefficient of α−1 which multiplies the mixture entropy, η m . This co-dimension distinction highlights the role the mixture entropy plays in selecting energetically favored structures. A second, apparent, distinction between co-dimensions is the α pre-factor on the reduced energies, however this coefficient is balanced by the fact that for a fixed solvent volume, V , the area (or length or number) of co-dimension α structure is where M 1 (τ ) = ∞ −∞ φ 1 (z; τ )dz, and M α (τ ) = ∞ 0 φ α (R; τ )R (α−1) dR, for α = 2, 3, are the solvent volume per unit co-dimension α structure.
For many parameter values, the FCH energy supports coexistence of co-dimensional structures. A composite morphology, u, comprised of non-intersecting bilayers, pores, and micelles with respective solvent volume fractions V α leads to an over-all free energy The free energy associated to each co-dimensional α structure has a pre-factor, σ α (τ )/M α (τ ), that depends upon the binary-well tilt, τ , multiplying an integral that depends upon the curvature and the hydrophilic and mixture entropy parameters. It is important to remark that the curvature terms arise only from curvature of the centerline or plane. For example the radial curvature intrinsic in a pore structure does not explicitly contribute to the free energy, except through integrals of φ 2 . We emphasize that these reduced expressions do not account for end-caps and junctions, which are essential features of the morphology depicted in Figures 2-4. Moreover, it appears that interfacial curvature only increases the free energy so that minimizing morphologies would seek to be flat. However interfaces often grow while constrained by connections to other parts of the network which may be too "inertial" to be pushed, so that bending or folding is the only mechanism available to extend the interfacial area. This is particularly true in PFSAs where the interface much grow within a confining elastic matrix formed by the crystalline phase of the backbone.

Applications of the Binary FCH Energy
We present applications of the binary mixture energy to Small Angle X-ray Scattering (SAXS) data for well-hydrated Nafion, and an extension of the energy to incorporate of ionic conduction in a framework that preserves the gradient flow structure.

A Study of Nafion: Comparison to SAXS Data
SAXS data provides a standard characterization of morphological structure in Nafion based upon the strength of scattering intensity as a function of wave-length. For Nafion SAXS data typically displays an "ionomer peak" associated to the network structure of the solvent phase at the 2-4 nm length scale. On a longer length-scale, the inhomogeneous distribution of the ionic groups, and the pretreatment of Nafion, often by extrusion, have their impact. Regions of Fluorocarbon backbone with low ionic density may crystallize, forming a stiff scaffolding which provides the elastic structure within which the water and the higher ionic content, uncrystallized polymer phase-separate [39,40]. The crystalline network occupies the 10-20 nm length-scale, often refereed to as the "matrix knee" [6]. However in cast forms of Nafion that have not been extensively annealed, the matrix knee may not be present [41].
Unfortunately SAXS is an indirect measurement and there is considerable ambiguity with respect to its interpretation. There has been controversy regarding the morphology of the water network as numerical and experimental investigation of the microstructure of Nafion and related PFSA membanes have lead to the proposition of bi-layer morphologies [42,43], cylindrical pores [44], an inverted pore morphology with solvent groups surrounding cylinders of crystallized backbone [6,45], spherical clusters [39], as well as more complex morphologies suggested by atomistic simulations [3,35,46,47]. Regarding the crystalline domains, proposed models include lamellar models [48], the rod-like model [45], a fringed-micelle model [40], and bi-modal networks [49]. The work of Schmidt-Rohr and Chen [44], was the first to use simulations of SAXS curves from three-dimensional data to systematically derive a model for Nafion structure. Their model features long parallel water channels in cylindrical inverted micelles, and long crystallites which are parallel to the water channels. Their simulated SAXS data yields an impressive agreement with experimental SAXS data. However, as the authors note, "the great length of the crystallites, more than 50nm, is surprising and needs to be critically checked in future work".
We demonstrate a fit to SAXS data, comparable to that obtained in [44], from a more amorphous, thermodynamically motivated, network structure. To compare the quasi-minimizers of the FCH to data on well-hydrated Nafion, we simulate two networks, a water-pore network driven by the hydrophilic interaction with the high-ionic group phase of the polymer network, and a crystalline network, which generates longer length-scale correlations. The water-pore network (blue) is generated using a 20% water volume fraction, and = η h = 0.03 and τ = −0.4, see Figure 8. The crystalline network (green) is generated from a second simulation with well tilt, τ = −0.1 and larger values of . The numerical SAXS data is computed using the procedure developed by Schmidt-Rohr [44]. In the absence of the crystalline network, we observe a SAXS curve with a single peak, but no matrix knee, while the inclusion of the crystalline network induces a matrix knee, see Figure 8C. A limitation of our approach is that the two networks are computed independently, and do not interact. However the crystalline network is substantially more sparse, and the interactions are limited to a small fraction of the spatial domain.  [45] in a log(I) versus log(q) plot for Nafion at 20 vol% of H 2 O, and simulated SAXS curve from a water pore network produced by the FCH energy (black dotted line) and SAXS data from a superimposed water pore network and a crystalline network, both obtained separately from the binary FCH energy (solid red line).

Extensions to Incorporation of Ionic Conduction
It is unclear to what extent imposed electric fields, which drive protonic current flow in Nafion, have an impact on the morphology. However it is clear that the protonic flux drags water and can generate macro-scale deformations of the material [50]. The development of detailed macroscopic models of electro-osmotic water drag, and its impact on the distribution of water within the membrane, have been a long-running focus of modeling efforts [51]. We present a framework which incorporates the interaction of long-range electrostatic forces with a dynamic charged interface, producing a flow which dissipates a total energy. To do so, the model must balance capacitive effects associated with the size and composition of the solvated region and the variation of the dielectric constant with the structure of the interface, against the longer range electostatic forces accounting for ionic motion. The minimal amount of self-consistency that one can ask of a model of this complexity is that the evolution decrease the free energy.
The structure of charge layers and the impact of solvent and dielectric on this structure has been actively discussed in the literature for more than 50 years. In the 1970's and 1980's Kornyshev introduced a family of nonlocal models for the dielectric which suggested its decomposition into high-frequency electronic transitions leading to fluctuations in the dipole moment, fast oscillations due to infrared intramolecular vibration of dipoles, and Debye fluctuations of the dipole orientations that hinder rotation of dipoles [52][53][54]. A final resolution is not yet at hand within the literature [55][56][57]. We do not propose a form for the functional dependency, rather we propose a framework in which the functional form should reside. In spirit we view the electrostatic interactions splitting into two distinct parts, a near-field that induces the complex structure of the solvated layers and couples to the interfacial morphology-this has been accounted for within the FCH, and a far-field potential, the potential of mean-force, which drives the motion of the counter-ions within the solvent phase. Moreover we derive an expression for the entropic effects within the potential of mean force. In the context of protons within a hydrated pore, the entropic force serves to confine the protons within the high-dielectric regions of the solvent, see Figure 7, where the long-range force induces their net flux.
To develop an energetically self-consistent model we appeal to a general framework [5], which couples elecrostatic energy to configurational and entropic energies while engendering an energy dissipating flow. We incorporate the long-range electrostatic field, φ, within an action functional, A and define the induced free-energy to be the action evaluated at the electric field which satisfyies the Poisson-Boltzmann equation. We show that a gradient flow generated by this action engenders novel coupling between the phase field and the charged groups while dissipating the induced free energy. The action is written in terms of the density of positive charges, p (hydronium for Nafion), the density of negatively charged ionic groups, which for simplicity we take as a prescribed function of the phase function u, n = n(u), and the dielectric constant ε = ε(u), and the electric potential φ. Guided by Kornyshev, we take the dielectric, ε to be slaved to the morphology of the mixture. The dielectric is small within the interface were waters solvate the ionic groups, but is large in regions of bulk-like water. Scaling distance x by the interfacial length scale ξ int , we propose an action which balances entropic and electrostatic contributions, Here q is the unit of electrostatic charge. A novel element of the action is the that Gibbs-Boltzmann entropy for the protons is taken to be relative to the hydrogen bond density, h = h(u), which is not spatially uniform but correlates with the dielectric ε(u). The entropy drives the protons to diffuse not to a spatially uniform state, but towards a state that gives them uniform occupancy of potential hydrogen bonds within the water network. The far greater propensity for bulk-like waters to form hydrogen bounds is incorporated by taking h much larger for u ≈ 1, i.e., in the core of the water-network. The key to generating an energy dissipating flow is requiring φ be a critical point of A, yielding Poisson's equation and subsequently evaluating the action at its critical point φ, which yields the energy Within the framework [5] the gradient flow generated by the action couples the phase field, u, to the electric field φ, and the protonic density p, where the potential of mean force (PMF) is given by The system Equations (22), (24) and (25) incorporates effects which are difficult to derive from force-balance laws. The coupling of φ and p to the evolution of the phase field u represents the influence of capacitive effects on the morphology. Conversely the proton density p satisfies a Nernst-Planck (NP) equation in terms of potential of mean force which depends upon entropic effects arising from solvation, finite volume, and electrostatic interactions of solvent with the tethered ionic groups. This is not unrelated to the classical models of excess chemical potential, based largely upon finite volume effects, such as proposed by Bikerman, see [56] for a detailed discussion. Most important however is that under homogeneous boundary conditions, the flow generated by A dissipates E Ionic [58], This framework is independent of the nature of the dependency of n, ε, and h upon u, and leads to a rich class of models with novel coupling of the electric field to the solvent motion. Extensions of this model to address water drag within an ionomer membrane could be achieved by incorporating convective flux. This requires replacing the ∂ t terms in Equations (24) and (25) with the material derivative ∂ t + v · ∇ and coupling to a momentum balance for the velocity field v. Models of this form which dissipate a total energy have been constructed, see [5,31]. However, validation of such a model would require an understanding of the effective coupling between water molecules and protonic groups, as well as the dependence of the fluid viscosity on the fluid state, i.e., bulk-like water or waters of solvation, which are points of active discussion in the molecular literature. The numerical resolution of convective models in a full 1+3D setting presents additional challenges.

Multi-Component Extensions
Binary models of ionomer membranes have natural limitations. While they incorporate competition among ionic groups for scarce solvent phase, it is less natural to incorporate local variation in densities of ionic groups as may arise from swelling of the pore network as a membrane imbibes liquid from an external source. In effect, the binary model assumes a constant charge density on the pore walls. Indeed, this is motivation for considering the binary phase model with periodic boundary conditions which fix the net ratio of water to ionic groups. In this section we outline a multi-component extension of the FCH energy which accommodates a diverse range of chemical interactions between components via linear constitutive dependencies of the functionalization parameters τ , η h , and η m on the local mixture composition. In particular, applications of the multi-component model to casting and annealing processes are outlined.
For a generic N + 1 phase model, the multi-well potential W m confines the volume fractions U i of the i'th phase to lie in the admissible set A = {0 ≤ U i ≤ 1 and U 1 + · · · + U N +1 = 1}. The second constraint represents incompressibility, and permits the N + 1'st phase to be slaved to the other phases. The independent variable is thus U = (U 1 , . . . , U N ), and the multi-component Cahn-Hilliard energy takes the form where D = diag{D 1 , . . . , D N } is a diagonal matrix of relative mobilities and W m is the multiphase potential well. The interfacial structures which determine the morphology are, at leading order, the solutions of the multi-component Cahn-Hilliard critical point equation The multi-component Functionalized Cahn-Hilliard takes the form where the parameters E, τ, η m and η h depend upon the phase field vector U .
As an example, consider a three component mixture of solvent with a hydrophilic polymer, labeled type 1, and a hydrophobic polymer, labeled type 2. The three component well W m depending upon U = (U 1 , U 2 ), is depicted in Figure 9, with three local minima at points 1, 2, and 3, corresponding to equiibria states of phases U 1 , U 2 and U 3 respectively. The phase U 3 denotes the solvent, while U 1 and U 2 denote the respective polymer species. The solutions of the critical point Equation (29) which connect the equilibria are depicted as the curves A, B, and C. By adjusting the value of η m (U ) for U along the curve A we can tune the energetic preference for the co-dimension of this connection, choosing bilayers, pores, or micelles without affecting the other interfaces. To reflect the hydrophobicity of polymer 2, we energetically disfavor the C interface by tuning η h to be negative along this curve while it is positive along A. For a mixture with polymer U 1 in the majority phase, the minimizers would consist of a network of solvent within the U 1 phase. Complementary to the solvent network will be a region of polymer 2 within the polymer 1 domain. In this transparent manner, we develop a model for a solvent which forms a "wetting" network within the hydrophilic polymer, while the hydrophobic polymer avoids contact with the solvent. In addition, the morphology of the solvent and the hydrophobic polymer networks can be independently adjusted by tuning the average values of η m and η h along the orbits A and B respectively. As an extension of the well-tilt parameter, we can adjust the characteristic radius of the pore or micelle component of the networks by modifying the relative values of W m at the local minima points, that is the differences W m (1, 0)−W m (0, 0) and W m (1, 0)−W m (0, 1). Moreover, only the integral of the parameter values over the interfacial critical orbits impacts the morphology; it suffices to consider a linear variation with the composition, for example η m (U ) = η m 0 + η m 1 · U for a constant η m 0 and a vector η m 1 ∈ R 2 . Moreover this model can form the basis for a mass-preserving transport mechanism for the ionic phase within PSFA membranes. Considering two polymer phases with distinct ionic group densities, n 1 > n 2 , the parameters of the FCH can be functions of the local ionic density n = U 1 n 1 + U 2 n 2 . The resulting flow drives the ionic groups via a local force balance which couples the ionic group density back to the evolving network morphology.

Casting and Annealing of PFSAs
The chief among these is that the ionic groups are not uniformly distributed, and the high-ionic group regions may experience rearrangement as hydration levels change. This effect is particularly significant during the casting and extrusion processes used to form the PFSAs. Indeed, Nafion membranes have typically been formed by a high-temperature melt extrusion process, which serves to form a network of crystalline structures which provide the membrane its structural rigidity. Recently the industry has shifted to a relatively low-temperature solution casting to produce the thinner membranes required in modern PEMFC applications, which has unfortunately resulted in a compromise in the crystalline content of the resulting membranes, which may disintegrate rapidly when exposed to boiling methanol, in contrast to the old extruded membranes which swell but maintain their mechanical integrity.
It is well understood that membrane casting and annealing procedures affect the morphological development of the crystalline phase within PFSA films and membranes [59,60]. It is desirable to develop processing and thermal annealing procedures for controlling morphological organization in order to systematically vary the degree of crystallinity and order within the ionic domains, however the time-scales for casting greatly exceed those accessible to particle-based simulations. Indeed, the FCH energy represents one of the few viable routes to computation of casting and annealing of PFSA membranes.
Within the context of the casting process the inhomogeneity of the ionic groups which functionalize the polymer backbones can be modeled by considering two polymer phases, a high ionic-group density polymer with volume fraction U 1 , a low ionic-group density polymer with volume fraction U 2 , and a volume fraction of casting solvent, U 3 . The U 2 polymer is substantially less hydrophilic, perhaps even hydrophobic, and its solvent-polymer interface, labeled C in Figure 9 is energetically disfavored by taking η m negative along the line C. The evolution of the mixture is a gradient flow on the FCH energy, As the casting process is time sensitive, it is important to track the mobility M of the phases, which depends upon the phase composition, U , and the temperature, T . Indeed, high solvent volume fraction, U 3 , and high temperature increase the mixture mobility, while the mobility of the crystalline phase, point 2 on Figure 9, is zero. The rate of evaporation of the solvent from the exposed boundary of the domain, Ω, is a key control parameter in a casting process. As solvent evaporates, the gradient flow draws the high ionic-group polymer to the solvent-polymer interface and embeds the low ionic-group density polymer within the polymer matrix, where at sufficient purity it crystallizes. However temperature plays a significant role. At sufficiently high temperatures, the crystalline phase will melt, the barrier in W m between the 1 and 2 phases along the critical path B disappears, see Figure 9 (right). The goal of annealing is to raise the temperature sufficiently so that the increased mobility permits a phase separation of the two polymer species within the time-frame provided by evaporation of the solvent. However the temperature cannot be too high that the crystalline domains melt. Once the low-ionic content phase separates from the high-ionic content phase, lowering the temperature locks in the crystalline phase. Importantly, the energetics of the multi-species phase separation, controlled by the parameters η m and η h , are also typically temperature dependent, and the development of a percolating network of crystallites, verses a disconnected domain which will not induce structural robustness to the membrane, will depend sensitively upon this dependence. As well, a percolating solvent phase is desired to enhance the solvent extraction, reducing residual solvent in the final cast membrane. It is exactly these constitutive and thermal dependencies of the parameters that must be obtained by upscaling of particle-based simulations, or from experimental data. In particular it is intriguing to consider the design of experiments to quantify the constitutive dependencies of the functionalization parameters.

Conclusions
The FCH modeling framework extends the scope of the continuum variational approach to incorporate the influence of the solvation entropy of ions on the network morphology of phase separated materials. Its energetic framework captures the dynamics of network formation, including the merging and pinch-off events required to form or degrade a network structure. It affords competition between co-dimension one (bilayer), two (pore), and three (pearled-pore) networks, and provides for a clear selection mechanism based upon the entropic competition for solvation between co-and counter ions and an effective interfacial stretching driven by electrostatic considerations. Most importantly, the model is computable within the framework of existing CPUs, in a fully three-dimensional setting, to physically meaningful length and time-scales.
We have shown that the FCH energy naturally couples to standard models for electrostatic energy as well as mass and momentum balance, and is readily extensible to multi-component mixtures. These models are amenable to studies of membrane casting as well as hysteresis arising from time dependent, inhomogeneous boundary conditions. A promising avenue for future research involves the coupling of the FCH models to particle based simulations as a means to identify the constitutive relations between local composition and the functionalization parameters, see [14,61].