We develop a general framework for transport of polyions, solvent and salt, with intended application to Layer-by-Layer (LbL) assembly of polyelectrolyte monolayers (PEMs). The formulation for the first time includes electrostatics, chemical potential gradients, and mechanical stress gradients as driving forces for mass transport. The general model allows all species to be mobile throughout the process and avoids the assumptions of stepwise instantaneous equilibrium and/or immobilized structures typical of previous approaches, while reducing to these models in appropriate limits. A simple constitutive equation is derived for a mixture of oppositely charged polyelectrolytes that accounts for network strand dilution and cross-chain ion pairing by appending reactive terms to the Smoluchowski probability diffusion equation for network strand end-to-end vectors. The resulting general framework encompasses the Poisson equation describing the electrostatic potential distribution, an osmotic pressure balance, a stress constitutive equation, and a generalized flux law of polymer transport. The computational domain is split into a PEM phase and an external solution phase with an appropriate boundary condition derived for the interface between the two. The mobile species (water and small salt ions) are taken to be in a state of dynamic equilibrium with their distributions enslaved to the perturbations in the two polyion compositions. The proposed model captures the swelling response of PEM films to external solutions. For the first time, we studied the effects of the temporal evolution of electrostatic and stress distribution on the rate of chain loss and absorption during rinsing and dipping of an idealized and arbitrarily selected and rigid brush layer into external solutions. The temporal evolution provides a kinetic basis for the ability of LbL films to grow under conditions that thermodynamics alone suggests would cause them to be washed away and to account for partial desorption during washing. The proposed transport framework constitutes a solid basis for eventual quantitative modeling of LbL assembly and transport in polyion networks more generally.

## I. INTRODUCTION

Oppositely charged polyelectrolyte (OCP) assemblies have attracted ever more attention in soft matter applications due to their ease of fabrication and versatility, especially in growing polyelectrolyte multilayers (PEMs). PEMs are typically manufactured via the popular Layer-by-Layer (LbL) assembly technique where two OCPs are alternately deposited onto a suitably primed substrate by either dipping the substrate in each polyelectrolyte (PE) solution or by spraying the solution onto the substrate with possible rinsing steps with pure buffer intervening between two consecutive polymer deposition steps. For simplicity, we refer to the polymer deposition step as the “dipping step” regardless of the actual method used to do so, while the term “rinsing step” is reserved for exposure of the PEM to polymer-free buffer solution. The primary constituents of PEMs include water-soluble positively and negatively charged PEs referred to as “polycation” and “polyanion,” respectively. Additionally, counterions, coions, water, and occasionally some desirable payload for which the PEM serves as a vehicle or reservoir make up the rest of the OCP assembly. Promising biomedical applications of OCP assemblies have emerged over the past two decades in areas that include DNA condensation for gene delivery,^{1,2} non-viral gene therapy,^{3–6} delivery of various biotherapeutics,^{7–12} cell encapsulation^{13–15} and culture,^{16} tissue engineering,^{17} and, for PEMs in particular, in microencapsulation,^{18} nanocomposite assembly,^{19} selective patterning,^{20} enzyme-active coating,^{21} drug delivery,^{22,23} and sensor fabrication.^{24}

The formation rate and performance of PEMs largely depend on the transport of components in and out of the film during the LbL process, which in turn is controlled by the physiochemical conditions such as pH, ionic strength, PE concentration, mixing ratio, etc.^{25} Despite the burst of applications, only a handful of theoretical studies have been published. The time and length scales involved in a typical LbL process, namely, 10-1000 s, and 1-1000 nm, respectively, put it well beyond the scope of ordinary atomistic molecular simulations. Even faithful coarse-graining attempts have yet to contend with extreme sensitivity to PE chemo-specificity of short-ranged electrostatic interactions such as charge regulation, cross-chain ion pairing, etc.^{26} There is, however, a great need for a continuum dynamic model of the LbL process to be used as a robust predictive tool in design and optimization of the PEMs.

Previous continuum models of LbL fall under two categories: (1) sequential *diffusion* and (2) *sequential equilibrium* models; each with their strengths and shortcomings, the most notable of which is briefly reviewed here. Following the pioneering work of Lavalle *et al.*^{27} which confirmed that vertical chain diffusion across a portion, or all, of the LbL film can occur depending on dipping conditions, Hoda and Larson^{28} proposed a sequential diffusion transport model that introduced transient chain diffusion to predict the exponential growth of an LbL assembled PEM and a transition to linear growth as they became thicker, while other phenomenological models^{29,30} predict exponential growth without relying on chain diffusion. Sequential equilibrium models treat the LbL process as a sequence of equilibrated, yet irreversible, adsorption steps.^{31–34} They are in essence generalizations of an equilibrium self-consistent field theoretic (SCFT) treatment of PE adsorption on charged substrates^{35,36} to a multi-step process.

The bound charge distribution inside the film is not explicitly accounted for in the determination of electrostatic potential in the sequential diffusion model but it is included in the SCFT treatment, which, however, overlooks the dynamics altogether from the outset. In the sequential diffusion model of Hoda and Larson,^{28} only one of the PEs, chosen to be the polycation, was assumed to have enough mobility for diffusion deep into the PEM. In both sequential diffusion and sequential equilibrium models, the electrostatic potential distribution has been taken to be stationary throughout any given deposition step. Neither type of model considers PEM swelling/deswelling and instead imposes the accumulated composition profile from preceding steps as a frozen constraint. Neither model allows for real-time tracking of film thickness controlled by the net effects of in and out diffusion of components nor PEM re-dissolution of weakly bound chains. Such transient film thickness predictions are needed for comparison with experimental measurements such as QCM-D (quartz crystal microbalance with dissipation monitoring) tracking of LbL assembly.^{37} Instead, in both types of models film re-dissolution can only be modeled phenomenologically at the conclusion of each step.

Stress evolution is an important feature of OCPs that is entirely missing in previous models that is the key, in our view, to improving their realism by providing, in particular, an explicit accounting of swelling/deswelling and partial, time-dependent, film redissolution. However, a constitutive equation for OCPs incorporating modulus variations due to solvent imbibition and ion-pair formation and break-up does not exist to date. Here, we first develop a general framework for transport of mixtures of multiple charged polymers and small ions that include chemical potentials, electric fields, and mechanical stresses, utilizing the Rayleighian formalism of Doi and Onuki.^{38} The two previous modeling strategies, sequential equilibrium and sequential diffusion, then become special cases of our proposed model. Next, we lay out an appropriate, though oversimplified, single-mode constitutive equation for stress in a mixture of two OCPs that could be generalized to a more realistic multi-mode sticky diffusion model in the future. We then focus on a few illustrative examples after the general model is simplified for a 1-dimensional Cartesian system appropriate for typical LbL experiments to highlight the crucial features of transport in idealized dipping and rinsing steps, deferring to future work the rigorous development of a robust numerical strategy for solution of the combined set of nonlinear and highly coupled equations.

## II. MODEL

### A. General formulation

The general conservation equation in a non-reactive system of water (*W*), simple ions (“+” and “−”) in addition to two PE chain types of opposite charges (*A*: polyanion, *C*: polycation) assuming no species volume change on mixing is given by the following equation:

where $\varphi i$ is the volume fraction of species *i*. Each species velocity vector *v*_{i} is determined from a flux law representing a balance of the friction and diffusion driving forces exerted on individual species. For a multicomponent mixture of simple fluids, these flux laws take the form of Stefan-Maxwell equations in which the diffusion driving forces are related to species velocities through a mobility matrix involving Onsager coefficients. However, it is not obvious how to extend the Stefan-Maxwell equations to incorporate elastic stresses and electrostatic fields subject to the incompressibility constraint. Thus, rather than using the Stefan-Maxwell equations, we derive the relevant flux laws in a self-consistent manner using the Rayleighian *R* of the system, Eq. (2), following the multi-fluid model of Doi and Onuki.^{38} The variational derivative of *R* with respect to each component velocity yields consistent dynamic equations that govern the system. Consequently, the dynamics of any non-inertial system is described through minimization of *R* with respect to species velocities comprising the system

The overdots above imply the time derivatives of the terms on the right side, where *S* and *F* denote the system total dissipation and free energy, respectively. Acceleration could also be included in *R* by adding a net momentum flux and momentum rate of change with time. However, inertial effects can be safely left out of the equations for sluggish motions of PEs over the time scales of interest to us.

The incompressibility constraint further implies that the volume fractions of species add up to unity, as given in Eq. (3), and the volume-averaged velocity $v\u2261\u2211\varphi ivi$ is then solenoidal

Introducing a Lagrange multiplier *p*, the incompressibility constraint is enforced in the Rayleighian, Eq. (4), where ** r** denotes the position vector. In addition to minimizing with respect to the species velocities, the Rayleighian now also needs to be minimized with respect to

*p*,

The rate of mutual frictional dissipation between two species is proportional to the square of their relative velocity. Another source of energy dissipation originates from deformation gradients of purely viscous small molecules, particularly of the solvent, characterized by a Newtonian viscosity $\eta $. We model polymers as typical bead-spring chains, in which stretching of the spring due to velocity gradients leads to elastic stress arising from the conformational entropy which thus should be incorporated in the free energy functional. Here we take these springs, as usual, to be instantaneously elastic so that polymer dissipation occurs only through frictional motion of their beads past the surrounding molecules. However, deformation gradients only dissipate energy through instantaneously viscous components, characterized by viscosity $\eta i$, which include the salt ions and solvent, but do not include the PEs. With this exclusion, $S\u0307$ from the two dissipative processes is given by the following equation:

Here *B*_{ij} are the elements of the inverse of the Onsager coefficient matrix. These are proportional to the mutual friction coefficient between *i* and *j*, denoted by $\zeta ij$, and in the case of no volume change on mixing, these are inversely proportional to the mutual diffusivity according to the following equation:

Here $v$ is a reference volume taken later to be the volume of a water molecule. We should note that in writing Eq. (6), we use the multicomponent diffusivities *D*_{ij} in mixtures, defined by $Bij=CkBTxixjDij$ which utilizes total molecular concentration *C* and species mole fractions.^{39} Note that the values of *B*_{ij} or *D*_{ij} might also depend on molecular weight and solvent quality and in general depend on temperature. If the reference volume is taken to be the same for each species (except for the salt ions, which we here take to be point charges), and there is no excess volume of mixing, we obtain Eq. (6).

The free energy of systems involving PEs, in particular of OCPs, is an emerging area of research in molecular thermodynamics, and a definitive closed form free energy functional does not exist as yet^{.}^{40} Here we represent the free energy functional as a combination of a mean-field contribution that encompasses electrostatic field and non-local effects, an elastic free energy stored in polymer chains, as well as a non-mean-field contribution from correlated electrostatic fluctuations

The mean-field free energy, $fMF=fFH+fGr$ in turn, encompasses the conventional Flory-Huggins (FH) mixing energy, electrostatic, and non-local (square gradient, Gr) free energy contribution expressions, which will be given later in Appendix B. *f*_{corr} denotes the correlation free energy arising from chain connectivity and electrostatic field fluctuations. Taking the time derivative of *F*, we arrive at the following form:

For the sake of brevity in notation, we leave off the subscripts conventionally used to indicate which variables are held fixed in the variational and partial derivatives in Eq. (8); instead, we note here that all variables, except for the one with respect to which the derivative is taken, are held fixed. The time derivative of the elastic free energy contributed by the two polyions A and C is now expressed as

Defining electrochemical potentials $\psi i\u2261\delta F\delta \varphi i$ and substituting for the time derivatives of volume fractions in Eq. (8) using the species conservation Eq. (1), the time derivative of free energy is now given by the following equation:

The electrostatic potential, mixing free energy, and non-local effects contribute to the electrochemical potentials $\psi i$ in Eq. (10). The variational derivative of *R* with respect to velocities yields the force balance on each component in the absence of inertia

The Onsager reciprocal principle stipulates that *B*_{ij} = *B*_{ji}. The last term in Eq. (11) is the stress driving force and is given by Eq. (12) for the solvent and salt molecules, while those of the polymeric species need to be determined according to an appropriate constitutive equation, which we derive in Appendix A,

In Eq. (12), $\eta i$ for each salt ion and for water should be equal to each other, and taken to be the viscosity of the salt-in-water solution that constitutes the non-polymeric species in our model. (This assumption is of no consequence here since we will end up neglecting the viscous stresses from these three species.) The final dynamic equations involving the two PEs are similar to those derived by Doi and Onuki in their two-fluid model,^{38} with the exception that they did not consider electrostatics, considered a mean-field network stress and restricted their discussion to binary systems, i.e., a polymer and a solvent, or two polymers. Extension of the two-fluid model to a multi-component system has been proposed for a system of solvent/non-solvent/neutral polymer where the polymer viscoelasticity was ignored.^{41} The current development is thus a generalization of the two-fluid model to charged multicomponent systems involving more than one polymeric species that also allows for stress relaxation. For an M-component system, there are M + 1 unknowns, namely, M velocities and the Lagrange multiplier *p*. The incompressibility condition, Eq. (3), along with Eq. (11) applied to each component, affords the requisite M + 1 equations, once an appropriate constitutive equation for stress and a free energy density expression are chosen. The mean-field electrostatic field $\Phi $ is obtained by minimizing *R* or equivalently *F* with respect to $\Phi $, which, as we demonstrate in Sec. II C, gives the Poisson equation if electrostatic correlations are ignored

### B. Constitutive equation for elastic stress

The present development is not meant to provide a rigorous constitutive equation which would require deriving self-consistently both the electrochemical potential $\psi i$ and the stress $\sigma i$ from a microscopic theory, in addition to inclusion of multi-mode polymer relaxation dynamics. Here we will focus on developing a single-mode constitutive equation for elastic stress relaxation that accounts for changes in modulus due to swelling by solvent and time-dependent changes in PE concentration. We consider a network of strands of opposite charge, temporarily bound to each other by ion pairing, which can break and reform, ultimately giving rise to relaxation of polymer configurations and therefore of stress. (Note that this “ion pairing” that we assume to be responsible for elastic stress is not included in the expression for free energy used here since we wish to avoid overcomplicating the theory at this point. However, a complete theory would need to include ion-pairing effects in both the electrochemical potential and the constitutive equation for mechanical stress, in a self-consistent way.) In the absence of composition gradients, this would give rise, in the simplest case, to an ordinary upper-convected Maxwell (UCM) equation for the stress. However, there are two polymeric species (polycation and polyanion) that are inter-diffusing, which leads to gradients of the density of ion pairs, and therefore variations in the modulus. In addition, there is swelling of the network by solvent, leading to dilution of the modulus. The dependence of modulus on network junction concentration produced by affine network swelling differs from that produced by diffusion of PE, and hence care must be taken in generalizing the UCM model to account, even qualitatively, for these effects. For this reason, we need to re-derive the constitutive equation from a starting point that is general enough to account for both deformation and spatial variations in network density. In this section, we present the following constitutive equation that we have derived in Appendix A for a mixture of OCPs, accounting for solvent imbibition and complexation:

where *i* can be A or C and the subscript “(1)” represents the upper convected time derivative defined by Eq. (A3) in Appendix A. The above constitutive equation at the first order level resembles the conventional UCM model and in fact reduces to it if the deformations are volume preserving, and the network reaction due to ion-pairing is absent. While the above may seem to be a rather trivial modification from the ordinary UCM equation, it contains two significant differences, both of them related to the modulus *G*_{i} = *n*_{i}*k*_{B}*T*, which is a function of time, through its dependence on the network density *n*_{i}. Both diffusive motion of the two PEs and swelling make the modulus time dependent, captured, for example, by setting *G*_{i} proportional to $\varphi A\varphi C$ assuming instant local equilibrium of ion pairing in the weak ion pairing limit; see Appendix A. However, modulus change due to swelling is accompanied by deformation, and this is captured by the upper convected derivative, modified by the volume-change term $\u2207\u22c5vi\sigma i$. Modulus change due to strand diffusion changes *G*_{i} but does not involve the term $\u2207\u22c5vi\sigma i$. Correct development of these terms, and accounting properly for these two ways in which network modulus changes, required starting from the general Smoluchowski probability equation, along with the network creation/destruction terms, as laid out in Appendix A.

### C. Model reduction

Here we lay out our simplifying assumptions and apply them to the general framework in Sec. II A and summarize the fundamental equations that make up our model. We adopt a one dimensional geometry in the direction normal to the film, say z and assume that the film is translationally invariant in the x and y axes. Consequently, only the zz-component of the tensor variable needs to be considered hereafter. All bold-faced vector and tensor variables in Sec. II A are thus replaced with their ordinary z and zz components, respectively (e.g., $vi\u2261vi,z$, $\sigma i\u2261\sigma i,zz$, and so forth) throughout the remainder of this work. Let us introduce the exchange electrochemical potentials as $\psi iw\u2261\psi i\u2212\psi W$. In Appendix B, we show that the pressure gradient in the force balance on each PE type Eq. (1) can be eliminated in which case it now reads for the polycation and polyanion as follows, respectively:

Here we have defined the net diffusion driving force $di$ for each PE type, *i* = *A*, *C*. In deriving Eq. (15), we have assumed that water and simple salts move through the system much faster than the mutual diffusion time scale of the two PE types past one another such that these mobile species are in a state of dynamic equilibrium at all times only to move in response to perturbations of either PE composition due to diffusion. The rapid equilibration assumption thus decouples the fast processes from the dynamic equations, which facilitates the numerical solutions by allowing much larger time steps to be taken. However, it is possible to derive explicit expressions only for relative polymer velocities using Eq. (15). This is a direct consequence of requiring rapid equilibration of water and salt ions, which leave the relative velocity of the polymer chains as the only dynamically relevant term. The rapid equilibrium assumption leads to another equation stemming from the overall force balance, as derived in detail in Appendix B, which for a special case of constant dielectric constant simplifies to the following algebraic form:

Here, $\Pi \u0303$ denotes the zz-component of the non-dimensional osmotic pressure tensor defined in Appendix B by Eq. (B15), while $\Phi \u0303\u2261eo\Phi kBT$ denotes the non-dimensional electrostatic potential, where *k*_{B} and *T* represent the Boltzmann constant and temperature, respectively, and *e*_{o} denotes the elementary charge. In Eq. (16), there appears the non-dimensional zz-component of the Maxwell electrostatic stress given by $\sigma \u0303E=l38\pi \u2113B\u2207z\Phi \u03032$, where *l* = 0.31 nm denotes the cube root of the volume of a water molecule and we have defined the Bjerrum length, $\u2113B\u2261eo24\pi \epsilon kBT$ which quantifies the relative importance of electrostatic and thermal energies. Similarly, $\sigma \u0303C$ and $\sigma \u0303A$ are non-dimensional zz components of polycation and polyanion elastic stress according to $\sigma \u0303i\u2261\sigma il3kBT$. The last equation is easily recognizable in the context of swelling of hydrogels with buffer, with the exception that the electrostatic stress is included, although it is negligible except at boundaries. However, we need to retain $\sigma \u0303E$ in Eq. (16) owing to the presence of thin interfaces between the film and the solution and at the substrate.

We further take the simple salt ions to be point charges occupying no volume for simplicity even though this assumption can be readily relaxed at additional numerical cost. Minimization of *R* with respect to electrostatic potential is carried out in Appendix B leading to the Poisson equation below:

The charge valence of a PE, even a strong one, is generally a function of composition; i.e., there is charge regulation. However, here, for simplicity, charge valencies *z*_{i} are taken to be insensitive to composition and electrostatic potential given that correlation free energies are not considered here. As a result, we take all polymer repeat units hereafter to be monovalent and fully ionized, i.e., $zC=\u2212zA=1$. Note that we have used the number concentration of the salt ions, $C\xb1$, rather than their volume fractions, as they do not occupy any volume.

Recalling the adoption of a one-dimensional geometry, the working constitutive equation for our model, Eq. (14), reduces to the following form after some rearrangement:

In general, the modulus *G*_{i} = *n*_{i}*k*_{B}*T* is proportional to the *i*th PE strand density, *n*_{i}, which is assumed to be equal to the density of ion pairs. Here we are assuming that for every strand of the polycation, there is a strand of the polyanion to which it is bound. There can be differing concentrations of the two polymers, which would lead to longer strands of the more concentrated polyion, but equal numbers of strands of each since each junction is assumed to bind together a polycation to a polyanion chain. The modulus can be affected by the instantaneous equilibrium between ion pair formation and break-up reactions, which produces changes in *G*_{i} additional to that produced by network dilation. As discussed in Appendix A, the rapid equilibrium assumption implies an instantaneous response of the modulus to changes in composition so that the modulus is slaved to the concentration of polycations and polyanions at any given position. This requires that ion pairs form and break up faster than the transport or deformations of polyions. The equilibrium ion-pair density, and thus the strand density, in a mixture of two OCPs has been shown to be a highly nonlinear function of individual PE composition as well as of the salinity, charge density, etc.^{26,42} However, at low enough salinity and in the weak ion-pairing limit where the fraction of PEs binding with oppositely charged repeat units is much smaller than unity, $ni\u221dKip\varphi A\varphi C$, where *K*_{ip} is the thermodynamic equilibrium constant of the ion pairing reaction;^{26} see also Appendix A for details. Therefore, we can set $Gi=Gip\varphi C\varphi A$ in this limit, which allows for the variation of modulus due to the complexation reaction, as well as due to swelling and de-swelling in our model. Since both PEs share the ion-pairs, *G*_{A} = *G*_{C}, so that the total mixture stress and modulus are simply twice those of individual PEs. We note that if entanglements, rather than ion pairs, between PEs dominate the modulus, then each *G*_{i} would be proportional to *ϕ*_{i}(*ϕ*_{A} + *ϕ*_{C}) such that *G*_{A} + *G*_{C} ∝ (*ϕ*_{A} + *ϕ*_{C})^{2}.

## III. RESULTS AND DISCUSSION

The proposed dynamic model in the equilibrium limit using appropriate parameters reduces to the sequential equilibrium model of LbL reviewed in the Introduction.^{12} In the very first step of LbL assembly, PEs from the solution adsorb onto the bare oppositely charged substrate. In the supplementary material, we demonstrate in detail that our model predicts equilibrium adsorption of PEs onto a flat solid surface, which has been extensively studied using SCFT and was the foundation of the sequential adsorption model of LbL assembly.^{35,36} If the coefficient of the square gradient term in the SCFT formulation, *κ*_{ii}, is set to $l212\varphi i$ in our model, which is valid for binary polymer-solvent systems and is also derived from a random phase approximation^{43} for $\varphi A\u226a1$ in the large wave vector limit, we recover the self-consistent-field theory (SCFT) result of Wang^{33} in the ground-state dominance limit for deposition of the first (base) layer at equilibrium. Now if the base layer is frozen so that its composition profile is not allowed to change when OCP invades, the overall force balance Eq. (16) becomes redundant. If we then solve equilibrium adsorption equations similar to those solved for the base layer, i.e., the Poisson equation and equality of the electrochemical potential of the invading PE throughout, setting the stress *σ*_{A}, *σ*_{C} = 0, and include the bound charges arising from the previously deposited base layer as a frozen charge distribution in the Poisson equation, we recover the second step of the sequential equilibrium adsorption model of LbL assembly.^{32,33} Subsequent adsorption steps under the same assumptions are similarly equivalent to the sequential equilibrium adsorption model.

### A. Equilibrium swelling

Water and simple salt ions always adopt the equilibrium distributions in our model as discussed in Sec. II C. PEMs are highly responsive to external stimuli such as changes in the bulk pH or salinity due to the redistribution of water and salt ions. In fact, every deposition step is typically accompanied by a change in bulk salinity and occasionally the pH, causing an instantaneous film shrinkage/swelling that precedes the slower diffusion of the PEs in and out of the film.

Equilibrium swelling of macroscopic ionizable gels where non-local and interfacial effects are absent is well understood.^{44} However, the interfacial and non-local effects cannot generally be ignored in thin LbL films. Nonetheless, the LbL models to date have not considered the film swelling/deswelling process. The present framework captures the swelling/deswelling of thin films with interfacial effects, although we only consider the equilibrium response in this section and focus on the case of a bilayer deposited in the first step of a typical LbL experiment adjacent to a substrate responding to a sudden change of the external solution with which the multilayer is in contact. It is imperative that the chain segments of the base layer form practically irreversible bonds to the substrate, at least relative to the time scale of the subsequent chain depositions. If not, film growth is thwarted, as equilibration with buffer during the rinsing step would wash the film entirely off the substrate. Accordingly, here, the first layer is taken to behave as an elastic permanently adsorbed PE “brush” layer (although the polymers can be adsorbed at positions other than their ends), the modulus of which is set by the areal number density of such bonds of PE to the substrate. For simplicity, we consider an initially dry bilayer of arbitrary and non-uniform composition that responds instantaneously to pure buffer.

To account for swelling, the equilibrium conditions in our model need to be imposed by solving the overall force balance Eq. (16) along with the Poisson equation, Eq. (17). We chose an initial coordinate frame *z*_{i} corresponding to the bilayer prior to swelling, and another final coordinate corresponding to equilibrium swollen state, *z*_{f}, with the transformation between the two, $\lambda zi$, denoting the local stretch ratio defined below, where “*i*” and “*f*” denote the initial (reference) state and final (equilibrium) state, respectively,

An elastic response implies that both chain types are stretched to the same extent and there is no relative motion between them, as would be expected if the ion pairing interactions remain intact during swelling. Mutual diffusion of one polymer with respect to the other commences immediately after the initial equilibration with buffer. The Poisson equation in the final state is now recast into the final reference coordinate system, using the transformation below:

We also need an expression for the moduli of each polyion type, which are proportional to the number density of their elastic strands times *k*_{B}*T*, in turn proportional to the density of network junction points formed through ion pairing. In the first bilayer, polyanions are anchored to both the substrate and the polycation chains via ion pairs. While the density of ion-pairs decreases upon swelling due to dilution of strands, the number of anchor points on the substrate is assumed here to remain unaffected by solvent ingress. We can thus take the dimensionless polyanion modulus $G\u0303A$ to be the sum of two contributions.

where $Gsub$ and $Gip$ have in the above been rendered dimensionless in the same way as $GA$ itself has been. The former contribution, $Gsub$, is proportional to the areal density of anchor points on the substrate while the latter, $G\u0303ip$, is a prefactor that is proportional to the equilibrium constant of an ion pairing reaction. Polycation chains on the other hand share the ion-pairs with the polyanion such that their elastic strands and thus its non-dimensional modulus is taken as $G\u0303C=G\u0303ip\varphi A\u2009f\varphi C\u2009f$. The total mixture stress to be used in Eq. (16), $\sigma A+\sigma C$, for the case of instantaneous swelling is obtained by setting the effective relaxation time of either PE in the constitutive Eq. (18) to infinity, i.e., in the elastic limit, and substituting the respective moduli with the proposed expressions for $G\u0303A$ and $G\u0303C$. The result derived in Appendix C is given below in the non-dimensional form

The first term on the right side of Eq. (22) is the contribution of fixed anchor points of the adsorbed polyanion layer to the substrate. This contribution, without which the film will be washed off eventually, would more rigorously be replaced by applying a boundary condition on polyanion composition and solving Eq. (16), including the gradient free energy terms, subject to the requirement that the total mass of the layer should be conserved upon swelling. Since the purpose of the present simulations is to illustrate the capability of the model in capturing the swelling/deswelling and the numerical challenges involved in the more rigorous approach, we capture the anchoring of the polyanion onto the solid substrate and chain connectivity during swelling by using a constant, uniform, modulus *G*_{sub}. Here, we only consider two contributions to the mean-field free energy density and thus to the osmotic pressure, namely, the Flory-Huggins (FH), given in Eq. (23a) and the square-gradient (non-local, Gr) contribution penalizing concentration inhomogeneities, given in Eq. (23b). In the latter, each term is weighted by respective stiffness coefficients *κ*_{ij} = *κ*_{ji}. Expressions for FH and square gradient free energy densities are given by Eqs. (B7) and (B9) in Appendix B, from which Eqs. (23a) and (23b) are derived. The resulting mean-field osmotic pressure $\Pi \u0303=\Pi \u0303FH\u2212\Pi \u0303Gr$ is the sum of a FH contribution, Eq. (23a), and one originating from non-locality arising from $f\u0303Gr$,

FH theory includes the van der Waals (VdW) interactions between species *i* and *j* characterized by χ_{ij}. We considered an initially dry and arbitrarily selected bilayer composed of polyanions and polycations alone. The volume fraction of the latter was picked as a decaying exponential given by $\varphi Ai=0.9e\u22122zi$ and a cutoff at the bilayer height of *z*_{i} = 2 nm and $\varphi Ci=1\u2212\varphi Ai$. In Fig. 1(A), the equilibrium structure, i.e., the final state, denoted by ( *f*), as well as the initial dry state, denoted by (*i*), of a bilayer 2 nm thick designated by the dashed curves is depicted using the parameters listed in Table I. Note that we chose a negatively valued $\chi CA$ to partly mimic the attractive short-ranged interactions between OCPs without having to account for nonlinearities associated with the determination of extent of ion-pairing,^{26} as discussed in Appendix A. Also, we consider the effect of chain stiffness coefficients on the swelling of the same bilayer examined in Fig. 1 with otherwise identical parameters in Fig. I-3 in the supplementary material. The film reaches an equilibrium thickness of nearly 6.4 nm, as inferred from Fig. 1(A). The final state is determined by a subtle balance of different forces in the bilayer, namely, osmotic pressure, electrostatic, and elastic stresses, the distribution of each of which is plotted in Fig. 1(B). Note that the osmotic pressure arising from non-local terms as well as elastic stresses are zero in the solution and their plot is cut off at around the equilibrium thickness ∼6.4 nm in Fig. 1(A). Any non-homogenies in FH osmotic pressure in the solution are created by Boltzmann distribution of simple ions which in turn is balanced by the electrostatic stress, as seen upon close inspection of Fig. 1(A) for z > 6.4 nm. The electrostatic potential distribution in the final state is plotted in Fig. 1(C). The overall charge in the film and solution balances the areal surface charge density set by the $ES/4\pi \u2113B$ ∼ 0.112 nm^{−2}, where $\u2207z\Phi \u0303(z=0)=\u2212ES$ is the boundary condition on the electrostatic potential. A cutoff value of $z=L\u221e$ was chosen such that the $\Phi z=L\u221e<10\u221210$ corresponding to an approximately zero electrostatic field at $z=L\u221e$, which is taken as the second electrostatic boundary condition. This approach is consistently used in the simulations discussed hereafter. The computational domain along the *z* direction is discretized with a non-uniform grid with a fine mesh close to the substrate ensuring that the results in Fig. 1 are insensitive to further grid refinement. We used an error relaxation method based on a damped Newton-Raphson method for solving the nonlinear discretized equations. In the vicinity of the solid substrate (*z* ∼ 0) and of the interface between the film and solution (*z* ∼ 6 nm), a “Debye layer” and a “double layer” storing electrical charge are established, respectively, as can be inferred from Fig. 1(D). An internal double layer overlapping with the substrate Debye layer in Fig. 1(D) emerges around *z* ∼ 1 nm where a change of the sign of net charge arising from both PEs occurs, which is due to the structure of the initially dry bilayer that is selected here. We believe that this is the first prediction of swelling by solvent and small ions of a dried adsorbed film containing two PEs. With proper determination of parameters, this approach could be used to predict swelling of a PE multilayer, as a first step in modeling LbL growth upon switching to the next LbL step.

$NA=NC$ . | $\kappa AC$ . | $\kappa AA=\kappa CC$ . | $\chi CW=\chi AW$ . | $\chi CA$ . | $ES$ . | $\u2113B$ . | T
. | $Gsub$ . | $Gipo$ . |
---|---|---|---|---|---|---|---|---|---|

1000 | 0 | 0.1 nm^{2} | 0.5 | −0.5 | +1 nm^{−1} | 0.71 nm | 300 K | $104\u2009Pa$ | $108\u2009Pa$ |

$NA=NC$ . | $\kappa AC$ . | $\kappa AA=\kappa CC$ . | $\chi CW=\chi AW$ . | $\chi CA$ . | $ES$ . | $\u2113B$ . | T
. | $Gsub$ . | $Gipo$ . |
---|---|---|---|---|---|---|---|---|---|

1000 | 0 | 0.1 nm^{2} | 0.5 | −0.5 | +1 nm^{−1} | 0.71 nm | 300 K | $104\u2009Pa$ | $108\u2009Pa$ |

### B. Polymer transport and relaxation dynamics

Dynamics enter the model through two mechanisms; namely, (1) the mutual diffusion between PEs governed by Eq. (15) and (2) stress relaxation set by the effective relaxation time of each PE type, $\tau eff,i$. Substituting $BAC$, given by Eq. (6), into either Eq. (15a) or Eq. (15b) and recognizing $Ji\u2261\varphi ivi$ as the z-component of the species volumetric flux, we arrive at the pseudo-binary flux law, Eq. (24), governing the PE diffusion process

Here, we have introduced the non-dimensional diffusion driving force for polycation migration $d\u0303C\u2261dCl3kBT$. Recall that $\psi \u0303CW$ is the nondimensional electrochemical exchange potential between polycation and water. Note that an equivalent and dependent diffusion equation could be written based on the polyanion diffusion driving force, using Eq. (15b) given below:

Recalling that $BCA=BAC$ and that Eq. (15) implies that $d\u0303C+d\u0303A=0$, it is readily verified that Eqs. (24) and (25) are linearly dependent. Since small molecules and ions exert negligible friction on either PE, the total force exerted on each PE type is exactly counterbalanced by the force on the other. It should be noted that water and salt counterions are free to move only in response to the slow diffusion of PEs. The salt counterion distribution is set by the Boltzmann distribution and the electrostatic potential distribution, as demonstrated in Appendix B. The overall force balance Eq. (16), incompressibility, and Poisson equations set the water distribution.

The flux law, either Eq. (24) or Eq. (25), coupled with Poisson Eq. (17), overall force balance, Eq. (16), the constitutive equation for mechanical stress, Eq. (18), and the PE conservation equation (1) provide the required framework for predicting two sequentially and cyclically occurring steps in LbL assembly, namely, rinsing with buffer and dipping in PE solution.

Due to complete lack or minimal concentration of either PE type in the solution phase during the rinsing and dipping steps, respectively, the full solution of the developed dynamic equation throughout the domain requires inordinately small time steps. It is thus sensible to treat each LbL step as a moving boundary problem, wherein PE conservation, the elastic stress constitutive equation, as well as the flux law for PE mutual diffusion are solved within the PEM film and are terminated by introducing a jump condition at the PEM-solution “interface,” while the Poisson equation and the overall force balance, Eq. (16), are still to be solved over the entire domain $0<z<\u221e$. While this method speeds the simulation, jump conditions in polymer composition across the interface pose numerical challenges, especially if the non-local terms in free energy are included.

In Secs. III C and III D, test cases are presented that allow us to reduce the mathematical complexity of our dynamic framework by considering only a bilayer and additionally freezing one of the PEs at the substrate and introducing a jump boundary condition on the mobile PE composition during the rinsing and dipping steps. This simplification allows us to focus on the role of time varying electrostatic and elastic stress fields in the rate of chain absorption and desorption. Consider an arbitrarily chosen polyanion brush layer of uniform composition 0.5 (50% hydration level by volume) with a constant thickness *L*. We take this polyanion layer to be firmly tethered to the substrate and its compositional profile to be frozen in place at all times such that the corresponding volume fraction profile for the polyanion is imposed as below:

### C. Dipping step

All derivatives of concentrations and stresses are discontinuous at the sharp PEM-solution interface. The condition that the polycation volumetric flux *J*_{C} is continuous and is non-infinite throughout the simulation at the PEM-solution interface can be used to derive a viscoelastic jump condition at the interface, Eq. (27), which is derived in Appendix D,

where we have used the far-field boundary condition on the electrostatic potential, i.e., $\Phi \u0303\u221e=0$ and split the non-dimensional polycation exchange electrochemical potential into a bare non-dimensional polycation exchange chemical potential, $\mu \u0303CW$, and a purely electrostatic potential, $\Phi \u0303$ as discussed in Appendix D. The structure of the solution phase, i.e., $L<z<L\u221e$, is resolved according to Eq. (28) which is coupled with the Poisson equation to be solved throughout the simulation domain, i.e., in the film and solution

Inside the film, i.e., $0<z<L$, polycation conservation should also be solved with Eq. (28) governing the instantaneous equilibrium polycation distribution on the solution side subject to the matching condition at the interface location $z=L$, Eq. (27). Even with the polyanion frozen, the rigorous simulation of viscoelastic diffusion of the polycation into the polyanion brush layer proved to be numerically challenging using the finite difference method. The primary reason lies with the intricate coupling between the electrostatic potential, the concentration of the polycation, and the mechanical stress borne by the incoming polycation chains as well as the fact that the polycation flux is set by the gradients in all three said parameters. The explicit finite difference schemes we employed were rarely stable using manageable time increments in simulating the dipping step. The implicit approach, however, could in principle address the instability of explicit methods but given the number of equations involved in our dynamic model and the number of grid points, the implicit finite difference methods would require the solution of a nonlinear set of equations involving a large number of variables. In order to tackle the stiffness problem, we here consider the infinitely fast relaxation limit of polycation stress; that is, the ion-pairs are taken to be so short-lived and/or so scarce (as might be the case in the high salt conditions) compared with characteristic diffusion time, that the polycation chains can move in a virtually stress-free state such that $\sigma C\u223c0$ in Eq. (24). (A slow mutual diffusivity *D*_{AC} is only consistent with a very fast effective relaxation time if *D*_{AC} is controlled not by the life-time of ion pairs but by the bare reptation, that is, snake-like diffusion, of PEs in the PEM.) While we are here neglecting stress relaxation during dipping, we will, however, consider stress relaxation for the rinsing step in Sec. III D 2. The overall force balance, Eq. (16) can now be used to derive the local elastic stress in the polyanion $\sigma \u0303A$, as a response function, since the concentration of the polyanion is fixed. Since the velocity of the polyanion is zero, the flux of the polycation can now be derived explicitly from Eq. (24) and is given by the following equation:

The polycation volumetric flux *J*_{C} has thus been reduced to a simple pseudo-binary (Fickian) flux expression with a generalized driving force. Equation (29) is similar to the flux expression derived by Hoda and Larson^{28} except that they used the regular (not the exchange) chemical potential of the polycation $\mu C$. A notable difference between the model of Hoda and Larson^{28} and the present work is the absence in our present work of any imposed electrostatic boundary condition at the film-solution interface, i.e., *z* = *L*, which implies that in our new work the Poisson equation needs to be solved dynamically in conjunction with conservation Eq. (1) and the simplified polycation flux equation, Eq. (29), throughout the simulation and thus extending to the solution. (Hoda and Larson^{28} did not consider the evolution of the electrostatic potential during a given LbL step.) Taking the solution side to be well mixed, the jump condition (27) on the polycation composition at the solution-PEM interface simplifies to that given by the following equation:

Here, with the assumption of fast relaxation of polycation stress, we consider the effect of degree of polymerization of the polycation, *N*_{C} and salinity of overlaying solution $C\xb1\u221e$ on transient chain absorption using $DAC$ = 0.1 (nm^{2}/s) for diffusion of the polycation into the rigid polyanion adsorbed layer described by the composition profile given by the Heaviside step function, Eq. (26), with $\varphi Ao=0.5$ and *L* = 10 nm. As depicted in Fig. 2, this idealized film is initially devoid of polycations. The far-field concentrations of salt and polycation in the solution are either both 1000 mM or both 100 mM with the former corresponding to $\varphi C\u221e$ ∼ 0.0002, which shows as nearly zero in the solution side (i.e., > 10 nm) in Fig. 2. The polycation concentration profiles and the corresponding height integrals *m*_{C} (*t*) of the polycation volume fraction profiles, for the two solution salt concentrations, are plotted in Fig. 2. These height integrals are reported in units of thickness (in nm) of the pure polycation. The electrostatic potential distribution profiles corresponding to panels in Fig. 2 are supplied in the supplementary material. Away from interfaces of the film with the substrate and the solution, the film remains electroneutral throughout the dipping period, while a double layer develops at the film free surface and a Debye layer develops in the vicinity of the charged substrate. A sample total charge distribution plot is provided in the supplementary material.

An implicit finite difference scheme and error relaxation were used to discretize and solve the coupled set of equations without the constitutive equation, Eq. (18), and overall force balance, Eq. (16). Given a polycation concentration at each grid point across the film at a given time step *p*, the discretized nonlinear equations were solved using a damped Newton-Raphson method for the polycation volume fractions in the film and electrostatic potential both in the film and in the solution at the next time step. A constant time step $\Delta t\u223c10\u22125\u221210\u22122$ (s) for a given simulation and non-uniform grid spacing of Δz ∼ 0.02-0.2 (nm) with finer meshing near the interfaces were selected for each simulation in this section and for the ones in Secs. III D 1 and III D 2, where the time step was chosen so that results in Fig. 2 are insensitive to further temporal and spatial grid refinements.

The driving force for diffusion of the polycation into a rigidly absorbed polyanion layer is not influenced by *N*_{C} for values significantly greatly more than 1, see Fig. I-6 in the supplementary material, as the chain translational entropy becomes vanishingly small. The main contribution to the driving force in this case stems from electrostatic attraction, which is taken to be insensitive to connectivity effects due to neglect of electrostatic correlation and VdW interactions in Fig. 2 and the negligible role of translational entropy. However, the external solution salinity does alter both the equilibrium polycation uptake as well as the diffusion driving force of the polycation and the uptake rate into the polyanion layer. The prediction in Fig. 2 that a PE deposition decreases at increased salinity has been experimentally observed at least for high enough salinity.^{45} In principle, the mutual diffusivity between OCPs increases with salinity. However, it was nevertheless held constant here to isolate the effects of chain length and salinity on the diffusion driving force alone, in the absence of their effect on the diffusion coefficient.

Just inside the film near the interface with solution, located at z = 10 nm, a discontinuity in the form of a sharp jump in the composition of the polycation persists at all times for simulations depicted in Fig. 2. This results from the imposed discontinuity in concentration of the underlying polyanion adsorbed layer as well as the instantaneous accumulation of polycations from the solution to the surface of the film, where an attractive negative electrostatic potential emanates from the polyanion adsorbed layer surface; see Fig. I-4 in the supplementary material. The electrostatic potential grows slowly less negative over time as the brush layer takes up positively charged polycations. The electrostatic potential eventually reaches equilibrium, which corresponds to the equilibrium polycation volume fraction profile given by the dotted lines for each salinity in Fig. 2. In other words, the polyanion absorbs as much the polycation from the solution as possible so as to neutralize its negatively charged repeat units with those of the polycation and thus minimize the electric potential inside the film, the value of which away from interfaces is commonly known as the Donnan potential. Polycation absorption is, however, opposed by expulsion of water (and salt ions, if they had been assigned a finite volume) from the film such that the balance between the electrostatic force drawing the polycation in and the excluded volume force exerted by water pushing the polycation out determines the equilibrium polycation distribution. For the selected choice of $\varphi Ao=0.5$ in Eq. (26), charge overcompensation over the entirety of the film is not achieved at equilibrium, as this would be physically impossible since the polycation volume fraction cannot exceed 0.5. However, the polycations from the solution extensively accumulate at the interface with the polyanion brush layer; see the composition profiles in Fig. 2. Even though a complete charge overcompensation of the brush layer considered here is not achievable by design, recall that our model reduces to the sequential adsorption model, which has been shown to capture charge overcompensation in favorable conditions, and would be achieved in our model if the brush layer had been obtained more rigorously from SCFT calculation. This scenario, however, would not offer any more insights into the dynamics.

In fact, the solution next to the film-solution interface contains roughly 85% and 60% by volume polycations at the beginning of the simulated dipping step at 100 and 1000 mM salt, respectively. This drops to 60% and 45% within a minute after the beginning of the simulated dipping step and stays more or less at the same level until the end of the dipping near equilibrium, as can be inferred from the dotted equilibrium curves in Fig. 2. (The progression of composition profiles toward the equilibrium dotted curve cannot be easily seen in Fig. 2 because of overlapping points.) This high concentration of the polycation in the immediate vicinity of the film surface is reminiscent of charge overcompensation which is believed to be a crucial factor for stable PEM formation during LbL experiments.^{46} The solution of the Poisson equation coupled with the diffusion equation in our model allows for our predictions to converge to the true equilibrium structure after long times during the dipping step. The framework developed here thus improves on the dynamic model of Hoda and Larson^{28} by allowing for the transient evolution of the electrostatic potential distribution consistent with polymer composition evolution both in the film and in the neighboring solution.

### D. Rinsing step

#### 1. Stress-free diffusion

In Sec. III C, we allowed the frozen polyanion layer to become equilibrated with polycation solution. In this section, we consider rinsing of the same polyanion layer after it has been fully equilibrated with each polycation solution as described in Sec. III C. Rinsing simulations are carried out here for films prepared under each of the two salinities corresponding to those of the dipping steps in Sec. III C, where the rinsing is carried out with a pure buffer solution of the same salinity used in the prior dipping step but devoid of either PE type. While the salinity of the buffer solution for the rinsing step need not be the same as that of the dipping steps, we make them the same here to compare the diffusion driving forces present during the dipping and rinsing steps under similar conditions. Thus, the equilibrium polycation volume fraction distribution inside the film, i.e., the dotted lines in Fig. 2, just before the beginning of the rinsing with pure buffer in each case is used as the initial condition here.

Since the buffer solution during the rinsing step is devoid of the polymer, the exchange chemical potential in the solution far-field $\mu \u0303CW\u221e$ tends toward negative infinity due to the diverging contribution of translational entropy of the chains in this limit. Therefore, Eq. (27) in the absence of the mechanical stress term requires that $\mu \u0303CWz=L,t$ should be negative infinity, which in turn implies that $\varphi Cz=L,t=0$. Unlike the dipping-step simulation, since there is now no polymer in the solution, there is no longer a need for determination of the structure of the solution dynamically during a rinsing step. The mutual diffusivity of entangled neutral polymers is inversely proportional to the molecular weight squared and we likewise expect a strong dependence of PE mutual diffusivity on molecular weight, especially given that we consider the weak ion-pairing limit in this section. However, here we take *D*_{AC} = 0.1 (nm^{2}/s) independent of molecular weight to assess the effect of molecular weight on the polycation diffusion driving force alone and on the simplified diffusion outflux given by Eq. (29).

The volume fraction profiles and the corresponding instantaneous polycation content *m*_{C} (*t*) versus time defined earlier are depicted in Fig. 3 corresponding to each simulated dipping step in Fig. 2 for two solution salt concentrations. The electrostatic field opposes the polycation out-diffusion at the interface while promoting the in-diffusion during dipping step. The main driving force for out-diffusion is provided by the entropy of water whose concentration is lower in the film due to the presence of polycations. In the simulated dipping steps in Fig. 2, the higher salt concentration leads to a lower equilibrium polycation content, expelling less water compared with the simulation at lower salt, in turn leading to a smaller driving force for water coming back into the film during the ensuing rinsing step considered in Fig. 3. This is also reflected in the faster rate of chain loss at lower salt inferred from Fig. 3, given the fact that we hold *D*_{AC} fixed and neglect all VdW and correlation effects in Fig. 3.

As with the dipping step, the chain length does not affect the polycation desorption if *N*_{C} ≫ 1 during the rinsing, over the time scales simulated here (again assuming a constant diffusivity). Interestingly, as shown in the supplementary material, the total charge distribution is qualitatively similar to that which develops during the dipping step with the film remaining electroneutral away from interfaces, while a negatively charged Debye layer in the vicinity of the positively charged substrate and an electrical double layer straddling the free surface develops and remains fairly static over time. Even though the equilibrium is not achieved until all residual polycations have escaped, the out-diffusion rate slows down substantially since *N*_{C} ≫ 1 in Fig. 3 as the attractive electrostatic potential grows progressively stronger as the polycation exits the film to the point that the polycation desorption appears to grind to a near halt after 10 min. This of course provides a purely kinetic explanation as to why multilayer films are not entirely washed away during rinsing steps. Moreover, the polycation desorption proceeds to ∼50% and 40% of the equilibrium state, i.e., *m*_{C} (*t* → ∞) = 0, for $C\xb1\u221e$ = 100 and 1000 mM, respectively, while both simulated dipping experiments in Fig. 2 achieve near equilibrium within 2 and 4 min for $C\xb1\u221e$ = 100 and 1000 mM, respectively, in Fig. 2. The faster completion of in-diffusion during the dipping step relative to the out-diffusion of polycations in the rinsing step, despite identical molecular weight and mutual diffusivity, is a direct consequence of the dynamically varying electrostatic potential created by the underlying polyanion adsorbed layer, which accelerates in-diffusion but impedes out-diffusion during rinsing and dipping steps. Our model thus provides the first direct prediction of this fundamental asymmetry between in-versus out-diffusion of PE, including the influence of polymer chains length and salt concentration. Since this dynamical asymmetry in the transport rate is the basis for the LbL assembly method, an ability to predict both its existence and dependence on solution parameters represents a fundamental advance in understanding, prediction, and design, of the LbL process.

#### 2. Effects of viscoelasticity on diffusion

Without any stress build-up during the rinsing step in Sec. III D 1, the polycation concentration at the sharp interface drops to zero instantly and the chains are washed away by the buffer solution. However, if the PE relaxation time is comparable to the characteristic diffusion time, then any release of the polycation into the solution at the surface is accompanied by a counteracting stress build-up that retards polycation out-diffusion. This physical picture follows from the interfacial jump condition given by Eq. (27), once we regularize the singularity arising from the zero polycation concentration in the solution. This singularity results in a divergence in the exchange electrochemical potential due to the contribution of polycation translational entropy. To obtain the proper boundary condition, we carry out a regularization in Appendix D and present the final result in the following equation:

For an infinitely fast stress relaxation, as was the case in Sec. III D 1, $\varphi CL=0$ follows from the above relation. Otherwise, the stored stress in the polycation chain maintains a positive PEM surface concentration despite being in contact with a polymer-free buffer solution. We now consider the same frozen polyanion brush layer given by Eq. (26) and used in Secs. III C and III D 1, and consider the effect of stress relaxation time on the dynamics of polycation release. We start the simulation with an initially uniform polycation volume fraction of 0.455. This value was achieved in the top two panels in Fig. 2 by exposing the frozen bare polyanion layer to a polycation solution in a dipping step until attainment of equilibrium. (We also used this initial concentration in our simulation of rinsing in Fig. 3.) We here consider both a fully charged polycation, *z*_{C} = +1, and a neutral chain, *z*_{C} = 0, to isolate the role of coupled electrostatics and stress relaxation. (The neutral chain concentration is artificial since it could not have arisen from a prior dipping step, in the absence of the electrostatic driving force for its entry into the film.) We use an explicit finite difference method to solve the coupled set of conservation, Poisson, and constitutive equations along with the jump condition at the PEM interface given by Eq. (31). Stable numerical solution of the coupled set of this viscoelastic diffusion problem, especially in the case of charged polymer, requires progressively smaller time steps for smaller relaxation times and in some cases also selective grid size refinement close to the film-solution interface. For instance, a stable numerical solution for reduced polycation relaxation times $\tau \u0303eff,C\u2261\tau \u0303eff\u22616\tau effDACNCl2$ of 0.1 and 0.001 required Δt = 10^{−2} and 10^{−5} s, respectively, for identical spatial grid resolution.

Figure 4 depicts the normalized polycation release profile defined as the ratio of the remaining polycation to the initial polycation content of the PE bilayer for various $\tau \u0303eff$. We intentionally held the mutual diffusivity fixed irrespective of the relaxation time and electrostatic charge valency so as to probe the impact of relaxation time alone on polycation transport, without the accompanying effect of reduced diffusivity. (Reduced diffusivity is likely to accompany the reduced relaxation time since both are controlled by the lifetime of ion pairs. Also, the neutral chain would likely to have a higher diffusivity since it would have no ion pairs, but we ignore these effects here.) The first observation from Fig. 4 is that the charged polymer escapes somewhat more slowly than does the neutral counterpart. The slower outflux of the charged polymer relative to its neutral counterpart is due to the progressively steeper uphill electrostatic potential forming during the course of rinsing encountered in the immediate vicinity of the PEM interface when the outgoing polymer is charged. A much greater slow-down would be seen if we allowed ion pairing to slow down the diffusivity of the charged case, relative to the neutral polymer.

The stress relaxation time is found to have a non-monotonic effect on the release of the polymer during the rinsing of both a fully charged polycation, see Fig. 4(A), and an analogous neutral chain, see Fig. 4(B). Below a critical relaxation time, the chain release is sped up by the finite relaxation time relative to that for the stress-free chain diffusion limit characterized by $\tau \u0303eff=0$. Above this critical value of the reduced relaxation time, namely, 0.01 and 0.1 for charge and neutral chains in Fig. 4, respectively, the stress build-up restrains the mobile chain segments at the interface that are bleeding out and swept away into the solution. Indeed, for the elastic limit of the constitutive equation, which is nearly reached once $\tau \u0303eff=10$ in Fig. 4, both with and without electrostatic driving force, nearly 70% of the polycation remains locked inside the film, held in by the elastic stress, and only able to escape when the built-up stress relaxes sufficiently. If the relaxation time were infinite, this polycation would be retained indefinitely. The entropic diffusion driving force dragging the polymer into the solution is thus eventually counterbalanced by the built-up stress (along with the electric field in the case of charged polymer), anchoring it to the underlying polyanion brush layer in the fully elastic limit. It is worth noting that even in the purely elastic limit, some chains escape. In this limit, the layer of the polycation is uniformly strained uniaxially to the extent necessary to generate an elastic stress sufficient to offset the osmotic stress driving the polycation out of the film. This uniaxial strain produces a dilation of the polycation layer so that it is thicker than the frozen polyanion layer, and the “overhang” of the polycation that extends beyond the polyanion layer is assumed to be in the free solution, and hence is dissolved away. In reality, if the polycation has an infinite relaxation time because it is bound by ion pairs to the polyanion, the overhanging polycation chains would remain tethered to polyanions by ion pairs and, though the polycation chains are stretched, they would not be free to diffuse into the far-field bulk solution. The retention of the overhanging polycations is a nonlocal effect that could be modeled using a SCF representation of the chains so that polymer composition would not jump to zero at a sharp interface as it does in the present case study, which leaves out the non-local terms in electrochemical potentials.

To understand the release profiles better, we singled out the two reduced relaxation times at which the release was the fastest for each of the two cases in Fig. 4 and examined the corresponding chain concentration and stress evolution distribution profiles inside the PEM in Fig. 5. The electrostatic potential distribution for the *z*_{C} = +1 is provided in Fig. I-5 in the supplementary material. The polymer concentration at the PEM surface is seen to remain greater than zero during the rinsing step, gradually tending toward zero as the chain stress relaxes away. While the polymer concentration profile at the substrate-PEM interface (*z* = 0) is flat in Fig. 5(B) due to a combination of uniform polyanion composition and the zero flux boundary condition there, the same zero flux boundary condition at the substrate along with a positive surface charged used in the simulations leads to an uptick with increasing z in the charged polymer concentration profile, see Fig. 5(A).

It can be inferred from Figs. 5(C) and 5(D) that the spatial gradient of the stress remains positive throughout the film peaking at the PEM free surface for both charged and neutral chains. This positive stress gradient creates a net negative diffusion driving force according to Eq. (24) and thus a net positive outflux that results in the boost in the polymer release rate observed for relaxation times below the critical values in each of the charged and neutral chain cases. The polymer mass is essentially “towed” to the free surface by the gradient in elastic stress, which is largest at the free surface and thus pulls the polymer there. For relaxation times longer than these critical values, however, the chain relaxation process at the PEM free surface becomes the rate limiting step such that the composition profile of the polymer becomes enslaved and nearly equilibrated to the concentration at the free surface and only drops when and if the polymer stress relaxes. In other words, the polymer chains diffuse to the surface, assisted by the stress gradient, only to be bottle-necked there by the slow relaxation of stress. Since the stress is zero in solution, the stress that pulls the polymer to the free surface also holds it there where the stress is maximum, until the stress can relax. This bottleneck leads to a decline in the polymer release rate. Even though the greater unrelaxed stress at the surface for simulations with long relaxation times creates a large diffusion driving force and outflux at the early stages of the rinsing step, as could be inferred upon close inspection of Fig. 4, the rapid outflux quickly fades, because of the build-up in stress and the concomitant drop in chain outflux. These results illustrate the subtle and rich phenomena present even in the oversimplified model of PE transport presented here. The results presented here should provoke careful experimental examination of the interesting effects of the interplay of electrostatics, mechanical stress, and diffusion in PEMs.

## IV. CONCLUSIONS

The next step toward a complete modeling of the LbL process involves allowing for the film-solution interface to be tracked dynamically by treating each step as a moving-boundary problem. In which case, two jump conditions, one on each PE type, during each of the rinsing and dipping steps, and the overall force balance equation will also need to be solved coupled to the other dynamic equations, thus allowing for swelling/deswelling of the PEM in real-time, even as other transport processes are occurring simultaneously. Our formulation is complete enough to allow this to be done, although the numerical challenges are formidable. Contributions to the free energy from short-ranged electrostatic interactions, including ion pairing of PE chains, counterion condensation of salt ions, electrostatic fluctuation free energy, as well as nonlocal free energy terms, could also be included in the free energy density functional, supplementing the mean-field expressions used here. In principle, once strategies to model the rinsing and dipping steps are developed given an initial PE composition profile and stress, the entire LbL experiment could be sequentially modeled, with the final state resulting from a previous step providing the initial condition for the subsequent one. However, a more reasonable strategy in the short and medium term would be to make detailed measurements of transport in single steps of the LbL process or polycation transport through pre-made PE films,^{47} to compare these with predictions of models such as ours, to confirm the validity of the basic approach, and to allow estimates of thermodynamic and transport coefficients to be made. We are undertaking such an effort at present and expect to begin comparing predictions with experimental data in the near future. A detailed numerical description dealing with moving and diffuse PEM-solution interphase and dynamic swelling/deswelling is left to a future contribution.

## SUPPLEMENTARY MATERIAL

See supplementary material for the effect of molecular weight and chain stiffness coefficient on the equilibrium polymer composition profile and electrostatic potential distribution for equilibrium adsorption of the polyelectrolyte onto an oppositely charged substrate as predicted by our model. Calculations in Fig. 1 have been performed for higher stiffness coefficients. Evolution of electrostatic potential distribution during rinsing and dipping steps corresponding to Figs. 2 and 5 is supplied. Negligibility of chain length for long chains on the dynamics of polymer sorption in the dipping step has been demonstrated. Finally, typical plots of total charge density for rinsing and dipping steps are presented.

## ACKNOWLEDGMENTS

The research presented in this paper was supported by the National Science Foundation under Grant No. 1707640. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of NSF.

### APPENDIX A: DERIVATION OF CONSTITUTIVE EQUATION

To develop a suitable starting point, we will combine the theory of spatially inhomogeneous dumbbells of Bhave *et al.*^{48} with the network theory of Yamamoto *et al.*,^{49–51} to produce a one-mode theory of stress relaxation in inhomogeneous networks of strands, where the dumbbells of Bhave *et al.*^{48} will be re-interpreted as mean-field network strands in our theory. A more complete theory would involve multiple relaxation modes and employ the “sticky diffusion” concept of Leibler *et al.*^{52} to link the breakage time of an ion pair to the distribution of relaxation times of the molecules and of the network.

In the kinetic theory of a dilute inhomogeneous solution of Hookean dumbbells, Bhave *et al.*^{48} developed the Smoluchowski equation, Eq. (A1), governing the evolution of the grand configuration density function $\Omega r,Q,t$ which depends on spatial coordinate vector ** r**, end-to-end dumbbell vector

**, and time t. $\Omega r,Q,t$ is the product of a local configurational distribution function $\omega r,Q,t$, which is normalized to unity when integrated over**

*Q***at each position**

*Q***, and a local number density of dumbbells $nr,t$. The subscripted nabla operator $\u2207Q$ in Eq. (A1) denotes the gradient in conformation space, while $\u2207$ denotes the conventional spatial gradient in**

*r***,**

*r*Bhave *et al.*^{48} derived the following constitutive equation by taking the second moment of Eq. (A1) in configurational space ** Q**. Their result is reproduced below up to the second order given by Eq. (A2), where the subscript “(1)” denotes the upper-convected derivative defined in Eq. (A3) for an arbitrary tensor

**and the superscript “t” denotes the transpose of the tensor**

*B*The dumbbell translational diffusivity in the solvent *D* is given by $kBT2\zeta $, where $\zeta $ is the drag coefficient of a dumbbell bead. The dumbbell spring constant is denoted by *H*, and $\tau =\zeta 4H$ is the relaxation time in this single-mode description of a dumbbell. The strain rate tensor ** D** is given by $2D=\u2207v+\u2207vt$. The number density of dumbbells

*n*whose centers of mass are located at position

**is obtained by integration of the probability distribution function Ψ over**

*r***space**

*Q*The velocity $vr,t$ in Eq. (A1) in the work of Bhave *et al.*^{48} is the mass-averaged velocity. The fluid, which is considered to be a simple solvent, is taken to be incompressible. Here, the relevant velocity that determines the deformations of each PE is the respective velocity $vi$ (*i* = *A*, *C*) for each of which $\u2207\u22c5vi\u22600$ since solvent imbibition causes the dilution of strands and a decrease in network modulus. Hereafter, we replace $v$ with $vi$ in all equations including the upper-convected derivative given by Eq. (A3). As a result, the diffusion equation (A1) can be developed for each PE individually leading to two constitutive equations, one for each PE type. Here, for simplicity, we assume that both PEs share the same modulus since the modulus derives from the ion pairing between them. Next, we recognize a few distinctions between dilute polymer solutions, for which the last three equations were derived, and a mixture of OCPs that form a network. First, the network strands in the latter are composed of two types of PEs that form a network through entanglements and ion-pairs. Second, LbL assembled PEMs are generally not dilute, and therefore the terms $ln\Omega $ in Eq. (A1) need to be replaced by a generalized thermodynamic potential function to account for non-idealities and the electrostatic field. However, for simplicity, here we leave the $ln\Omega $ in Eq. (A2) intact.

Finally, elastic strands can be created and destroyed through the formation and break-up of ion pairs between oppositely charged segments. Therefore, we need to include a reaction term in our modified probability density distribution equation to account for strand formation and disappearance. To do so, we borrow from Yamamoto’s network model,^{49–51,53} which is written below for each PE, with no spatial variations, as

The final two terms of this equation are, respectively, the rate of creation of new network strands of *i* = *A*, *C*, which have configurations of the equilibrium (no-stress) ensemble, $\omega o,i$, while the last term is the destruction rate, set by the time required for junction breakage $\tau o$ which is shared by both PE types. The above equation assumes that strands deform affinely until a junction breaks, after which the strand takes on an equilibrium configuration. Detailed balance, in which strands can relax, but are not destroyed or created, implies that $k=1/\tau o$. The Yamamoto model is most applicable to simple telechelic polymer networks with network-forming stickers confined to the two ends of the polymer molecule, where breakage of a single junction allows the connected strands to quickly relax completely, on a time scale of $\tau o$. For a PE network, however, with multiple stickers per molecule, loss of a single junction permits only partial relaxation. To relax the stress completely, even within a single strand, the entire molecule must relax, on a time scale much longer than the sticker lifetime $\tau o$. The chain as a whole relaxes according to Rouse or reptation diffusion dynamics, with diffusion rate set by the time constant $\tau o$, a process known as “sticky diffusion.” Theories for sticky diffusion are able to relate $\tau o$ to a much longer time scale $\tau $ for relaxation of the entire chain, and hence of the stress. To generate a one-mode model for relaxation of stress governed by breakage of ion pairs for each PE type *i*, in a spatially inhomogeneous network, we append a modified strand creation term and a destruction term to the diffusion equation of Bhave *et al.*^{48} given by Eq. (A1) for each PE type comprising the network

Again, the normalized probability distribution function $\omega i$ reverts at equilibrium to $\omega o,i$ and both of distribution functions integrate to unity in $Q$ space. In Eq. (A6), *k*_{i} is the rate of network formation per unit volume, which in turn depends on the local composition of the network. Note that, unlike the Yamamoto model, we do not allow breakage of a network junction to completely relax the strand since it is strongly constrained by the remaining junctions along the chain. Thus, we introduce the coefficient $\alpha i$ in Eq. (A6) as a microscopic parameter that controls the portion of the orientation that a freed strand relaxes before forming a new ion pair. The parameter $\alpha i$ is clearly a function of chain microstructure, chain length, and cooperativity of complexation, which could be determined by, e.g., molecular simulations. We have also replaced the dumbbell diffusion coefficient in the equation of Bhave *et al.*^{48} with *D*, the polymer network diffusivity. Integrating Eq. (A6) over $Q$ space yields the strand conservation equation up to second order gradient terms for each PE type

In the limit in which the exchange rate between newly formed ion pairs and the ones breaking up is much faster than the transport of the strands (which would be the case if there are many ion pairs per chain), the reaction term in Eq. (A7) dominates and we can write

The *i*th PE stress $\sigma i$ is related to its second moment of the probability distribution function $QQi$ according to Eq. (A9). At equilibrium, $\Omega i=\Omega o,i$ is a Gaussian function and the polymer stress vanishes. Accordingly, $QQi$ becomes isotropic at equilibrium and equal to $nikBTH\delta $,

Taking the second moment of the configuration probability distribution Eq. (A6), we arrive at an expression for the evolution of the second moment of each PE given below:

Substituting for $niQQi$ and *k*_{i} in Eq. (A10) using Eqs. (A8) and (A9) we arrive at the following constitutive equation with the third and higher order gradient terms dropped:

The second term on the left side of Eq. (A11) can be substituted out using the strand conservation Eq. (A7). We notice that the microscopic details embodied in the parameter $\alpha $ can be absorbed into an effective single-mode relaxation time $\tau eff,i\u22614H\zeta +\alpha \tau o\u22121$. Since the lowest order network stress is zero, at least the first order terms in the constitutive equation need to be retained, and given below after some rearrangement

Due to solvent imbibition and consequent swelling of the network, Eq. (A12) has an additional term $\u2207\u22c5vi\sigma i$ not present in Eq. (A2) which accounts for the decrease in modulus *G*_{i} = *n*_{i}*k*_{B}*T* as the network takes up more solvent. Formally, *n*_{i} is given by the product of the total density of *i*th repeat units and the fraction of ion-paired repeat units, $\beta i$ which depends on local composition, salinity, pH, etc.^{26} In principle, a kinetic (reaction) law accounting for formation and breakage of ion pairs needs to be coupled to the constitutive equations and thus the rest of the model. The ion pairing kinetic law would involve the equilibrium ion pairing constant that embodies chemospecificity of OCPs as well as the lifetime of a single ion-pair. However, we already assumed that the exchange rate of individual ion pairs is much faster than the characteristic transport time in deriving Eq. (A8). At the very beginning of rinsing or dipping steps, an instantaneous swelling such as the case presented in Sec. III A occurs causing rapid deformation without there being any significant relative transport of the OCPs. Following this, any subsequent deformation can only come from the relative transport of OCPs and the concomitant rearrangement of solvent and salts such that after the initial instantaneous swelling the deformations and transport occur at a similar rate. Due to the much faster exchange rate ion pairs, we can assume that ion pair density reaches a state of local equilibrium instantly. Even in this limit, $\beta i$ is a nonlinear function of composition, salinity, pH, etc. However, at low enough salinity for strongly charged OCPs and in the weak ion-pairing limit, where $\beta i\u226a1$, it can be shown that $\beta C\u221dKip\varphi A$ and $\beta A\u221dKip\varphi C$, where *K*_{ip} is the equilibrium ion pairing reaction constant. Finally, we can set $G=Gip\varphi A\varphi C$ if network junctions are predominantly made up of ion pairs.

For the purpose of transport in the LbL assembled films, the first order terms in the constitutive equation suffice, particularly given that the thermodynamic non-idealities and electrostatic potential, were left out of the probability distribution Eq. (A6) to begin with. The second order terms [i.e., the last two terms in Eq. (A12)] would be strictly accurate anyway only in the dilute regime for neutral polymer solutions. Dropping the second order terms in Eq. (A12), we arrive at the desired constitutive equation appropriate to ion-pair forming PE mixtures subject to swelling/de-swelling.

### APPENDIX B: MODEL REDUCTION

Approximating salt ions as point charges and the rapid equilibration assumption decouple them from dynamic equations simply because they pose no friction when moving in the system. In fact, as we shall see, the salt ion concentrations become enslaved to the electrostatic potential distribution. The force balance, Eq. (11), when applied to *i* = *W*, ±, namely, on water and simple ions denoted collectively by “±” in our shorthand notation, reduces to the simple forms below where the friction forces other than those involving polymer-polymer friction, as well as the viscous stress terms, have been dropped

Note that the Lagrange multiplier *p* does not apply to point charges anymore, as they occupy no volume. Retaining the friction term between the two polymer types, *i* = *A*, *C*, the force balance Eqs. (15a) and (15b), respectively, simplifies to the following:

As can be seen, the velocity $vW$ of the solvent decouples from these two dynamic equations so that, of the conservation equations in Eq. (1), we only need to consider those for the two polymers. The Lagrange multiplier *p* in Eqs. (B3) and (B4) is a response function that can be determined by solving for the composition profiles. Eliminating *p* and introducing exchange electrochemical potentials as $\psi iw\u2261\psi i\u2212\psi W$, Eqs. (B3) and (B4) are rewritten below:

Here we have defined the net diffusion driving force *d*_{i} for each PE type, *i* = *A*, *C*. Given a free energy density to evaluate exchange electrochemical potentials and a constitutive equation for *σ*_{i}, the two polymer conservation equations, Eq. (1) for *i* = *A*, *C* coupled with Eqs. (B3) and (B4) seemingly provide four independent equations to solve for the four field variables of interest, namely, $\varphi C,\varphi A,vC,vA$. However, it is only possible to derive explicit expressions for relative polymer velocities using Eqs. (B3) and (B4). This is a direct consequence of requiring rapid equilibration of water and salt ions, which leave the relative velocity of the polymer chains as the only dynamically relevant term. The next equation is obtained by summing the last two equations, yielding the overall force balance in the system

We will return to Eq. (B6) later and recast it in a more familiar form. We next turn our attention to the exchange electrochemical potentials by choosing explicit free energy densities in *F*. We point out that the particular choice of the free energy model here does not alter the general dynamic framework. In general, a definitive theory of correlation free energy of OCPs is still an emerging area and, despite recent progress,^{40} rigorous closed-form expressions requisite for a transport model are still lacking. Therefore, here we focus only on the mean-field contributions to the free energy in order to reduce the complexities arising from charge regulation, complexation, etc. The first contribution to the mean-field free energy is from the Flory-Huggins (FH) theory which includes the van der Waals (VdW) interactions between species *i* and *j* characterized by $\chi ij$ and the translational entropy of each component weighted inversely by *N*_{i} the ratio of the molecule size and the volume of a lattice site *l*^{3} taken to be volume of a water molecule ∼30^{3} Å^{3}. All parameters with an overlying tilde are non-dimensional, hereafter

Here, *V* represents the system volume. The mean-field electrostatic free energy of a molecule of component *i* with a charge valence of *z*_{i} in an electrostatic potential $\Phi $ is simply $zi\Phi $. Additionally, the energy stored in an electrostatic field is proportional to square of the gradient of the electrostatic potential. Therefore the normalized mean-field electrostatic free energy $f\u0303ele$ is given as

The charge valence of PEs, even strong ones, is generally a function of composition; i.e., there is charge regulation. However, here, for simplicity, charge valencies are taken to be insensitive to composition and electrostatic potential given that correlation free energies are not considered here. As a result, we take all polymer repeat units hereafter to be monovalent and fully ionized, i.e., $zC=\u2212zA=1$. The electrostatic potential in Eq. (B8) has been made non-dimensional according to $\Phi \u0303\u2261eo\Phi kBT$, where *e*_{o} denotes the elementary charge. Note that we have used the number concentration of the salt ions, rather than their volume fractions, as they do not occupy any volume. We have also defined the Bjerrum length, $\u2113B\u2261eo24\pi \epsilon kBT$ which quantifies the relative importance of electrostatic and thermal energies. Each lattice site size is taken to be volume of a water molecule ∼30^{3} Å^{3}.

The non-local (gradient) free energy *F*_{Gr} penalizes concentration non-homogeneities and leads to smoothly varying composition profiles across interfaces. The non-local free energy is proportional to the sum of the square gradients weighted by their respective stiffness coefficients *κ*_{ij} = *κ*_{ji},

The square gradient terms for non-polymeric components are not considered here. The normalized mean-field free energy is thus given by $f\u0303MF=f\u0303FH+f\u0303elec+f\u0303Gr$. Minimization of *R* with respect to electrostatic potential in Eq. (13) can now be carried out leading to the Poisson equation below:

The electrochemical potential $\psi i$ can be evaluated by taking the variational derivative of the free energy functional with respect to component volume fractions

The dimensionless chemical potential $\mu \u0303i\u2261\mu il3kBT$ encompasses all contributions to free energy except the electrostatic potential and is evaluated as follows:

The electrochemical potentials in the overall force balance Eq. (B6) can be nondimensionalized using Eq. (B11) resulting in Eq. (B13), where $\sigma \u0303\u2261\sigma l3kBT$,

Enforcing the rapid equilibration of simple salt ions by setting the variation of the free energy functional with respect to their concentration to zero, resulting in Eq. (B2), yields the Boltzmann concentration distribution as demonstrated below, given that the volumes of simple salt ions are taken to be negligible and the only contribution of simple salts to free energy is through their translational entropy and electrostatic energy

In the last equation, we have introduced the bulk salt concentrations, i.e., infinitely far from the film, $C\xb1\u221e$, where we take the reference electrostatic potential $\Phi \u221e$ to be zero. We can recast Eq. (B13) by introducing the osmotic pressure tensor defined as below:

where in general the above is a tensor equation if the free energies have spatial gradients due to nonlocal terms. Here, we consider the scalar form only, as we neglect nonlocal terms. Due to the linear dependence of mean-field electrostatic free energy $fele$ on composition in Eq. (B8), $\Phi $ does not contribute to the osmotic pressure and thus the variational derivatives can be replaced with exchange chemical potentials for the two PEs, $\mu AW\u2261\mu A\u2212\mu W,\mu CW\u2261\mu C\u2212\mu W$ and the regular chemical potentials for simple point-like ions, $\mu \xb1$. Additionally, $f=fMF$ is the only contribution to the free energy density in the absence of correlation terms. Due to the one-dimensional nature of our model, we only consider the zz-component of the osmotic pressure tensor in the non-dimensional form $\Pi \u0303=\Pi zzl3kBT$. In this case, taking the divergence of the osmotic pressure tensor simplifies to taking the gradient of $\Pi \u0303$, which allows us to convert Eq. (B15) to the following expression in the non-dimensional form where we have used the non-dimensional chemical potentials defined earlier

We have used the Gibbs-Duhem relation given by Eq. (B17) to obtain Eq. (B16) by substituting for the volume fraction of the water using the incompressibility constraint

We can now substitute for the chemical potential gradients in Eq. (B17) with their electrochemical counterparts using Eq. (B11) and the Boltzmann distribution of simple ions, Eq. (B14), where $\u2207z\mu \u0303\xb1=\u2213\u2207z\Phi \u0303$,

The right-hand-side of the above can now be eliminated using the overall force balance Eq. (B13),

We can further simplify the result by substituting out the terms in the parentheses in Eq. (B19) using the Poisson equation (B10), leading to the following:

A further simplification is possible if $\epsilon $ is taken to be constant, implying that the Bjerrum length $\u2113B$ is also constant throughout, in which case Eq. (B20) can be recast in a form that is readily integrated to yield the following:

Here, there appears the non-dimensional zz-component of the Maxwell electrostatic stress given by $\sigma \u0303E=l38\pi \u2113B\u2207z\Phi \u03032$. The last equation follows from the overall force balance, Eq. (B13) and is one that is easily recognizable in the context of swelling of hydrogels with buffer, with the exception that the electrostatic stresses at the interfaces are insignificant at a macroscopic level. However, we need to retain $\sigma \u0303E$ in Eq. (B21) owing to the presence of sharp interfaces between the film and the solution and at the substrate. Throughout the system, the difference between the local osmotic pressure and the bulk $\Pi \u0303\u221e$ is given by the sum of the Maxwell electrostatic and mechanical elastic stresses. Our one-dimensional model including Eq. (B21) could be easily extended to multi-dimensional problems in which case the zz-components on the left-hand-side can be replaced by their full tensorial form with the right-hand-side of Eq. (B21) multiplied by the identity tensor.

### APPENDIX C: ELASTIC LIMIT OF CONSTITUTIVE EQUATION

The response of the constitutive Eq. (18) in the elastic limit is obtained by setting the effective relaxation time $\tau eff,i\u2192\u221e$. We seek here the sum of the two PE stresses that enter the osmotic pressure balance given by Eq. (16). Since the ion pairs do not break up in the elastic limit, there is no relative movement of PEs inside the PEM. As such, both PE chain types move in concert so that they are stretched to the same extent so that $vC=vA$. Thus, we can solve for the elastic limit of individual PE stress $\sigma i(i:A,C)$ by dropping the term in Eq. (18) involving the effective relaxation rate of each PE type, resulting in the following equation:

Replacing the gradient of the velocity in Eq. (C1) by the local network stretch ratio $\lambda $ shared by both PE types gives

Now introducing $\epsilon \u2261ln\lambda $, Eq. (C1) can be rewritten in the following form where the time derivatives in (C1) cancel out. As expected, time does not affect the elastic stress, which is solely a function of $\epsilon $ given by the following equation:

As discussed in Sec. III A, the polyanion modulus permanently anchored to a substrate and forming ion-pairs with polycation chains is taken to be of the form $GA=Gsub+Gip\varphi A\varphi C$, while $GC=Gip\varphi A\varphi C$. A mixture of PEs is taken to be at rest in the initial state “*i*,” where $\epsilon =0$. Expressing the volume fraction in terms of local stretch ratio can be done according to Eq. (19) in the main text. Therefore, the individual moduli can now be expressed in terms of $\epsilon $ as below:

Plugging Eqs. (C4) and (C5) into Eq. (C3) and solving the resultant ordinary differential equation (ODE) for the elastic limit of the stress of respective polyions in terms of $\epsilon $, we arrive at the following general solutions:

Applying the initial condition $\sigma i\epsilon =0=0$, the integration constants *C* in Eq. (C6) and *C*′ in Eq. (C7) are obtained. The final result is obtained by adding up the last two equations and after some rearrangement is given by the following equation:

Replacing $\epsilon \u2261ln\lambda $ and multiplying through by $l3kBT$, we arrive at Eq. (22). Note that the result is markedly different from one of an incompressible UCM model $G(\lambda 2\u22121)$ in which the modulus *G* remains constant.

### APPENDIX D: JUMP CONDITION AT FILM-SOLUTION BOUNDARY

#### 1. Dipping step

We start by writing the formal definition of one-dimensional differentiations in the general polycation flux expression given by Eq. (29) as below:

Here, the dummy variable *h* is taken to be positive so as to evaluate the right-handed derivatives in Eq. (29). Multiplying both sides of Eq. (D1) by *h*, taking the limit, and recognizing that the stress outside the PEM film is zero and that $JC$ is finite, we arrive at the following:

Rapid equilibration and well-mixing on the solution side imply that the diffusion driving force vanishes everywhere on the solution side. Setting $d\u0303C$ Eq. (24) to zero in the solution then yields the following:

#### 2. Rinsing step

We split into two parts the exchange electrochemical potential of the polycation $\psi Cw\u2261\psi C\u2212\psi W$, recalling that the electrochemical potential is defined as $\psi i\u2261\delta F\delta \varphi i$. The splitting is done by separating the translational entropy, $ln(\varphi C)NC$ from all contributions to the non-dimensional exchange electrochemical potential of the polycation $\psi \u0303CW\u2261\psi CWl3kBT$ denoting all other contributions by $\psi \u0303CWex$, leading to

Plugging the above equation into the definition of the diffusion driving force [Eq. (24)], we get

Now a jump condition resembling Eq. (27) can be obtained by first plugging Eq. (D5) into the simplified polycation flux Eq. (29). We then use the formal definition of the right-handed differentiations in the result at the film-solution interface so as to approach it from the solution side by taking the limit of the dummy variable *h* → 0^{+} as below:

Since the buffer solution during the rinsing step is devoid of the polycation, we can write in the infinite dilution thermodynamic limit