A Molecular Thermodynamic Model of Coacervation in Solutions of Polycations and Oppositely Charged Micelles.

To guide the rational design of personal care formulations, we formulate a molecular thermodynamic model that predicts coacervation from cationic polymers and mixed micelles containing neutral and anionic surfactants and added salt. These coacervates, which form as a result of dilution of conditioning shampoos during use, deposit conditioning agents and other actives to the scalp or skin and also provide lubrication benefits. Our model accounts for mixing entropy, hydrophobic interactions of polycation with water, free energies of bindings of oppositely charged groups to micelles and polycations, and electrostatic interactions that capture connectivity of charged groups on the polycation chain and the micelle. The model outputs are the compositions of surfactants, polycation, salt, and water in the coacervate and in its coexisting dilute phase, along with the binding fractions and coacervate volume fraction. We study the effects of overall composition (of surfactant, polycation, and added salt), charge fractions on micelles and polycations, and binding free energies on the phase diagram of coacervates. Then, we perform coacervation experiments for three systems: sodium dodecyl sulfate (SDS)-JR30M, sodium methyl cocoyl taurate (Taurate)-JR30M, and sodium lauryl alaninate (Alaninate)-JR30M, where JR30M is a cationic derivative of hydroxyethylcellulose (cat-HEC), and rationalize their coacervation data using our model. For comparison with experiment, we also develop a parametrization scheme to obtain the requisite binding energies and Flory-Huggins χ parameter. We find that our model predictions agree reasonably well with the experimental data, and that the sulfate-free surfactants of Taurate and Alaninate display much larger 2-phase regions compared to SDS with JR30M.

Brownian dynamics simulations of telechelic polymer - latex suspensions under steady shear.

We carry out coarse-grained Brownian dynamics simulations of shearing flow of a colloidal suspension bridged by telechelic polymers with "sticky" end groups and vary sticker strength ε over a range from 3 to 12 in units of kBT, motivated by an interest in simulating the rheology of latex paints. The most extensive results are obtained for dumbbells, but the trends are confirmed for 3-bead trumbbells and chains of up to 11 beads. The numbers of colloids and of polymers are also varied over a wide range to confirm trends established for smaller, more computationally affordable, systems. The dynamics are the result of an interplay of the shear rate and three different times scales: the time τBridge for a sticker on a bridging chain to be released from a particle surface, which scales as exp(0.77ε), the time for the polymer chain to relax, τR, which scales as the square of polymer chain length, and the time τD for a colloid to diffuse a distance comparable to its own radius, R, which scales as R3. The scalings of the bridge-to-loop and loop-to-bridge times namely τBL ∝ exp (0.75ε) and τLB ∝ exp (0.71ε), are similar to those of τBridge, for ε values above around 5 kBT, because of the relatively short chains considered here (i.e., 60 Kuhn steps). However, τR becomes more dominant for longer chains, as shown by Travitz and Larson. The zero-shear viscosity η0 is estimated from the Green-Kubo relation, and found to scale as exp (0.69ε), similar to that of τBridge. A weak influence of η0 on τD is observed, with the influence expected to become stronger when τD becomes larger, as shown previously by Wang and Larson. At shear rates in the nonlinear regime, shear-thinning is found with exponents ≈ -0.10 to -0.60, and the first normal stress difference is positive, consistent with some of the experimental data of Chatterjee et al. on model latex paint formulations. The weakness of the shear thinning, relative to that of hydrophobically modified ethoxylated urethane (HEUR) solutions without colloids, is likely due to the observed insensitivity of the loop-to-bridge and bridge-to-loop transition times to the imposed shear rate. This preliminary study provides the first mesoscale simulations of these suspensions, useful for assessing and improving both more accurate multi-scale models and eventually constitutive equations for these complex suspensions.

Extracting free energies of counterion binding to polyelectrolytes by molecular dynamics simulations.

We use all-atom molecular dynamics simulations to extract ΔGeff, the free energy of binding of potassium ions K+ to the partially charged polyelectrolyte poly(acrylic acid), or PAA, in dilute regimes. Upon increasing the charge fraction of PAA, the chains adopt more extended conformations, and simultaneously, potassium ions bind more strongly (i.e., with more negative ΔGeff) to the highly charged chains to relieve electrostatic repulsions between charged monomers along the chains. We compare the simulation results with the predictions of a model that describes potassium binding to PAA chains as a reversible reaction whose binding free energy (ΔGeff) is adjusted from its intrinsic value (ΔG) by electrostatic correlations, captured by a random phase approximation. The bare or intrinsic binding free energy ΔG, which is an input in the model, depends on the binding species and is obtained from the radial distribution function of K+ around the charged monomer of a singly charged, short PAA chain in dilute solutions. We find that the model yields semi-quantitative predictions for ΔGeff and the degree of potassium binding to PAA chains, α, as a function of PAA charge fraction without using fitting parameters.

Surfactant desorption and scission free energies for cylindrical and spherical micelles from umbrella-sampling molecular dynamics simulations.

HYPOTHESIS: The free energies associated with adsorption/desorption of individual surfactants from micelles and the fusion/scission of long micelles can be used to estimate the rate constants for micellar kinetics as functions of surfactant and salt concentration. EXPERIMENTS: We compute the escape free energies △Gesc of surfactant from micelles and the scission free energies △Gsciss of long micelles from coarse-grained molecular dynamics simulations coupled with umbrella sampling, for micelles of both sodium dodecylsulfate (SDS) in sodium chloride (NaCl) and cetyltrimethylammonium chloride (CTAC) in sodium salicylate (NaSal). FINDINGS: For spherical micelles, △Gesc values have maxima at certain aggregation numbers, and at salt-to-surfactant molar concentration ratios R near unity, consistent with experiments. For cylindrical micelles, SDS/NaCl shows a minimum, and CTAC/NaSal a maximum in △Gesc, both at R ~ 0.7, while △Gsciss of CTAC micelles also peaks at around R ~ 0.7 and that of SDS micelles increases monotonically with R. We explain the non-monotonic dependence of escape and scission free energies on R by a combination of electrostatic screening and the decrease of micelle radius with increasing R. Transitions from predominantly spherical to cylindrical micelles, and between adsorption/desorption and fusion/scission kinetics with changing salt concentration can be inferred from the free energies for CTAC/NaSal.

Hysteretic Swelling/Deswelling of Polyelectrolyte Brushes and Bilayer Films in Response to Changes in pH and Salt Concentration.

Abstract We use a quartz crystal microbalance with dissipation (QCM-D) to investigate the swelling/de-swelling and hysteresis in brushes of weakly ionizable polyanion poly(acrylic acid) (PAA) brushes and bilayers containing a PAA brush and a poly(ethylene imine) (PEI) overlayer [...].

Free Energy Cost of Interdigitation of Lamellar Bilayers of Fatty Alcohols with Cationic Surfactants from Molecular Dynamics Simulations.

Cationic surfactant mixed with fatty alcohol as cosurfactant in excess water can form stable emulsions, known as "lamellar gel networks," that contain extended and interconnected networks of swollen bilayers, including ones with in-plane liquidlike disorder (Lα phase) and solidlike order (Lß phase). To study their structure and thermodynamics, molecular dynamics (MD) simulations with lateral pressure and temperature scans along reversible pathways were used to drive reversible phase changes, including formation at negative lateral pressure of the LßI phase with interdigitated tails of opposing leaflets. Thermodynamic integration, with extrapolations to infinitely slow scans, yielded a free energy difference between the interdigitated LßI and non-interdigitated Lß phases of 2.4 ± 0.5 kJ/mol, which is consistent with the spontaneous formation of the Lß phase under atmospheric pressure in simulation. Thermodynamic cycles involving temperature and lateral pressure for which the free energy difference is identically zero were constructed as negative controls to verify the method. Using lateral pressure, including negative lateral pressure, helps avoid kinetic bottlenecks that occur when temperature alone is used as the control variable. The method, using negative lateral pressure, should be widely applicable to other bilayers to identify molecular properties that control interdigitation and other bilayer properties.

Overcharging of polyelectrolyte complexes: an entropic phenomenon.

Overcharging in complex coacervation, in which a polyelectrolyte complex coacervate (PEC) initially containing equal moles of the cationic and anionic monomers absorbs a large excess of one type of polyelectrolyte species, is predicted using a recently developed thermodynamic model describing complexation through a combination of reversible ion binding on the chains and long-range electrostatic correlations. We show that overcharging is favored roughly equally by the translational entropy of released counterions and the binding entropy of polyelectrolytes in the polyelectrolyte complex, thus helping resolve competing explanations for overcharging in the literature. We find that the extent of overcharging is non-monotonic in the concentration of added salt and increases with both strength of ion-pairing between polyions and chain hydrophobicity. The predicted extent of overcharging of the PEC is directly compared with that of multilayers made of poly(diallyldimethylammonium), PDADMA, and poly(styrene-sulfonate), PSS, overcompensated by the polycation in two different salts: KBr and NaCl. Accounting for the specificity of salt ion interactions with the polyelectrolytes, we find good qualitative agreement between theory and experiment.

Forward Flux Sampling of Polymer Desorption Paths from a Solid Surface into Dilute Solution.

We compute desorption rates for isolated polymers adsorbed to a solid wall with a rare event sampling technique called multilevel splitting, also known as forward flux sampling. We interpret computed rates with theories based on the conjecture that the product tdesDRg2 of the desorption time tdes and diffusivity D divided by squared radius of gyration Rg scales with exp(h/Rg) where h is the equilibrium ratio of adsorbed surface concentration of polymer Γ to bulk concentration of polymer c. As the polymer-wall interaction energy is increased, the slope of lntdesDRg2 vs. NVMFkBT nearly approaches unity, as expected for strongly-adsorbing chains, where N is the degree of polymerization and VMF is the height-averaged monomer-wall interaction energy for a strongly adsorbed chain. However, we also find that this scaling law is only accurate when adsorption strength per monomer exceeds a threshold value on the order of 0.3-0.5 kBT for a freely jointed chain without or with excluded volume effects. Below the critical value, we observe that tdesDRg2 becomes nearly constant with N, so that tdes∝Nα, with α≈2. This suggests a crossover from "strong" detachment-controlled to a "weak" diffusion-controlled desorption rate as VMF/kBT drops below some threshold. These results may partially explain experimental data, that in some cases show "strong" exponential dependence of desorption time on chain length, while in others a "weak" power-law dependence is found. However, in the "strong" adsorption case, our results suggest much longer desorption times than those measured, while the reverse is true in the weak adsorption limit. We discuss possible reasons for these discrepancies.

Correction: Time-dependent shear rate inhomogeneities and shear bands in a thixotropic yield-stress fluid under transient shear.

Correction for 'Time-dependent shear rate inhomogeneities and shear bands in a thixotropic yield-stress fluid under transient shear' by Yufei Wei et al., Soft Matter, 2019, 15, 7956-7967, DOI: 10.1039/C9SM00902G.

Salt- and pH-induced swelling of a poly(acrylic acid) brush via quartz crystal microbalance w/dissipation (QCM-D).

We infer the swelling/de-swelling behavior of weakly ionizable poly(acrylic acid) (PAA) brushes of 2-39 kDa molar mass in the presence of KCl concentrations from 0.1-1000 mM, pH = 3, 7, and 9, and grafting densities σ = 0.12-2.15 chains per nm2 using a Quartz Crystal Microbalance with Dissipation (QCM-D), confirming and extending the work of Wu et al. to multiple chain lengths. At pH 7 and 9 (above the pKa ∼ 5), the brush initially swells at low KCl ionic strength (<10 mM) in the "osmotic brush" regime, and de-swells at higher salt concentrations, in the "salted brush" regime, and is relatively unaffected at pH 3, below the pKa, as expected. At pH 7, at low and moderate grafting densities, our results in the high-salt "salted brush" regime (Cs > 10 mM salt) agree with the predicted scaling H ∼ Nσ+1/3Cs-1/3 of brush height H, while in the low-salt "osmotic brush" regime (Cs < 10 mM salt), we find H ∼ Nσ+1/3Cs+0.28-0.38, whose dependence on Cs agrees with scaling theory for this regime, but the dependence on σ strongly disagrees with it. The predicted linearity in the degree of polymerization N is confirmed. The new results partially confirm scaling theory and clarify where improved theories and additional data are needed.

Time-dependent shear rate inhomogeneities and shear bands in a thixotropic yield-stress fluid under transient shear.

We study the rheological responses and shear-rate inhomogeneities and shear banding behaviors of a thixotropic fumed silica suspension in shear startup tests and flow reversal tests. We find that this suspension under transient shear exhibits not only viscoelasticity, yielding, kinematic hardening, and thixotropy, but also time-dependent shear inhomogeneities including bands when the apparent shear rate is below a critical value between 0.1 and 0.25 s-1. Through multiple shear startup tests and flow reversal tests, we find that thixotropy promotes flow heterogeneity while kinematic hardening suppresses it. We propose a simple thixo-plastic constitutive equation that can qualitatively predict the important features of the rheological response and banding dynamics in shear startup tests and flow reversal tests.

A metastable nematic precursor accelerates polyethylene oligomer crystallization as determined by atomistic simulations and self-consistent field theory.

Using PYS, TraPPE, OPLS-L, and Flexible-Williams (FW) force field models, atomistic simulations at temperatures ranging from 450 K to 600 K are performed to predict the melt density ρ, the persistence length Np, the nematic coupling constant α, and crystallization dynamics for pentacontane (C50). The coupling constant α arises from packing entropy of rodlike Kuhn segments and increases with increasing ρ and Np. Together with a self-consistent field theory, Np and α are then used to predict the isotropic-to-nematic (IN) transition temperature for polyethylene (PE) oligomers as a function of chain length. The nematic phase is found to be metastable since the IN transition temperature lies below the crystal melting temperatures for C50 in simulations using different force fields. Finally, isothermal simulations of crystallization for PE C50 oligomers and C1000 polymers show that crystal nucleation may be much accelerated by quenching below the IN transition temperature, where chains in the isotropic state first rapidly form nematic ordered domains, within which crystalline order then grows. We also find that the PYS, TraPPE, and FW models overpredict the melting temperature for C50 by around 50 K, while the most flexible OPLS-L model gives a melting temperature within around 10 K of the experimental value. Although giving a more accurate melting temperature, the slow crystallization kinetics of the OPLS-L model may limit its application in direct simulations of PE crystallization.

Computational self-assembly of colloidal crystals from Platonic polyhedral sphere clusters.

We explore a rich phase space of crystals self-assembled from colloidal "polyhedral sphere clusters (PSCs)," each of which consists of equal-sized "halo" spheres placed at the vertices of a polyhedron such that they just touch along each edge. Such clusters, created experimentally by fusing spheres, can facilitate assembly of useful colloidal crystal symmetries not attainable by unclustered spheres. While not crucial for their self-assembly, the center of the PSC can contain a "core" particle that can be used as a scaffold to build the PSC. Using Brownian dynamics simulations, we show the self-assembly of eight distinct crystalline phases from PSCs that correspond to the five Platonic polyhedra, and that are made of spheres with purely repulsive interactions. Strong crystalline order is seen in the centers of mass of the PSCs, or equivalently the core particles. The halo particles also may organize into crystal structures, usually with weaker crystalline order than the core particles. Notably, however, in crystals assembled from the octahedral and icosahedral PSCs, the halo particles are also well ordered, nesting within the crystals formed by the cores. Interestingly, despite the rounded nature of the PSCs, in some cases we obtain structures similar to those of the corresponding faceted polyhedra interacting only via excluded volume. Only the tetrahedral PSCs fail to self-assemble into a crystal, but we demonstrate that a pre-assembled crystal - whose halo particles sit on a close-packed face-centered cubic lattice, and whose core particles form a diamond structure - is stable at high density and melts into a hexagonal phase at lower density.

Slip-Spring and Kink Dynamics Models for Fast Extensional Flow of Entangled Polymeric Fluids.

We combine a slip-spring model with an 'entangled kink dynamics' (EKD) model for strong uniaxial extensional flows (with Rouse Weissenberg number W i R ≫ 1 ) of long ( M w > 1   Mkg / mol for polystyrene) entangled polymers in solutions and melts. The slip-spring model captures the dynamics up to the formation of a 'kinked' or folded state, while the kink dynamics simulation tracks the dynamics from that point forward to complete extension. We show that a single-chain slip-spring model using affine motion of the slip-spring anchor points produces unrealistically high tension near the center of the chain once the Hencky strain exceeds around unity or so, exceeding the maximum tension that a chain entangled with a second chain is able to support. This unrealistic tension is alleviated by pairing the slip links on one chain with those on a second chain, and allowing some of the large tension on one of the two to be transferred to the second chain, producing non-affine motion of each. This explicit pairing of entanglements mimics the entanglement pairing also used in the EKD model, and allows the slip spring simulations to be carried out to strains high enough for the EKD model to become valid. We show that results nearly equivalent to those from paired chains are obtained in a single-chain slip-spring simulation by simply specifying that the tension in a slip spring cannot exceed the theoretical maximum value of ζ ' ϵ ˙ L 2 / 8 where ζ ' , ϵ ˙ and L are the friction per unit length, strain rate and contour length of the chain, respectively. The effects of constraint release (CR) and regeneration of entanglements is also studied and found to have little effect on the chain statistics up to the formation of the kinked state. The resulting hybrid model provides a fast, simple, simulation method to study the response of high molecular weight ( M w > 1   Mkg / mol ) polymers in fast flows ( W i R ≫ 1 ), where conventional simulation techniques are less applicable due to computational cost.

From well-entangled to partially-entangled wormlike micelles.

We combine mechanical rheometry, diffusing wave spectroscopy (DWS), and small angle neutron scattering (SANS) with a simulation model, the "pointer algorithm", to obtain characteristic lengths and time constants for wormlike micelle (WLM) solutions over a range of salt concentrations encompassing the transition from unentangled to entangled solutions. The solutions contain sodium lauryl ethylene glycol sulfate (SLE1S), cocamidopropyl betaine (CAPB), and NaCl. The pointer algorithm is extended to include relaxation of unentangled micelles, allowing micelle parameters to be extracted from the rheology of partially entangled solutions. DWS provides the data at high frequency needed to determine micelle persistence length accurately. From pointer algorithm fits to rheology, we observe a salt-induced rapid change in micellar length as the solution enters the well-entangled regime and a weaker growth with surfactant concentration consistent with mean-field theory. At a lower surfactant concentration, micelle length and persistence length from SANS are roughly consistent with values from rheology once the lower surfactant concentration used in SANS is accounted for. This is, to our knowledge, the first time that quantitative comparisons of structural features including micelle length are made between rheology and SANS. Finally, scaling laws for micelle diffusion and recombination times indicate that micelle kinetics are reaction controlled leading to mean-field recombination with surrounding micelles over the entire range of concentration of interest except at very low and very high surfactant concentrations where either short micelles or branched micelle clusters are dominant.

Search for the Source of an Apparent Interfacial Resistance To Mass Transfer of CnEm Surfactants To the Water/Oil Interface.

Experiments have shown that relaxation of oil/water interfacial tension by adsorption of alkyl ethoxylate surfactants from water onto an oil droplet is delayed relative to diffusion-controlled adsorption. We examine possible causes of this delay, and we show that several are implausible. We find that redissolution of the surfactant in the oil droplet cannot explain the apparent interfacial resistance at short times because the interface will preferentially fill before any such redissolution occurs. We also perform umbrella sampling with molecular dynamics simulation and do not find any evidence of a free-energy barrier or low-diffusivity zone near the interface. Nor do we find evidence from the simulation that premicellar aggregation slows diffusion enough to cause the observed resistance to interfacial adsorption. We are therefore unable to pinpoint the cause of the resistance, but we suggest that "dead time" associated with the experimental method could be responsible-specifically a local depletion of surfactant by the ejected droplet when creating the fresh interface between the oil and water.

Tension-Induced Nematic Phase Separation in Bidisperse Homopolymer Melts.

We use an analytical mean-field theory and all-atom molecular dynamics (MD) simulations to predict that external tension, together with the nematic coupling interactions, can drive phase separation of long chains from short ones in bidisperse homopolymer melts. The nematic coupling parameter α for polyethylene (PE) oligomers under applied tension is extracted from the MD simulations and used in the mean-field free energy to predict the phase boundary for bidisperse melts in which the longer chains are stretched by uniaxial tension. The predicted phase diagram is validated by direct MD simulations. We also show that extensional flow, and possibly even shear flow, may lead to nematic phase separation in molten PE oligomers, because the flow can impose a stronger tension on the longer chains than the short ones.

A transport model and constitutive equation for oppositely charged polyelectrolyte mixtures with application to layer-by-layer assembly.

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.

Role of electrostatic correlations in polyelectrolyte charge association.

Reversible ion binding equilibria in polyelectrolyte solutions are strongly affected by interactions between dissociated ionic species. We examine how the structural correlations between ionic groups on polyelectrolytes impact the counterion binding. Treating the electrostatic correlation free energy using the classical Debye-Hückel expression leads to complete counterion dissociation in the concentrated regime. This unphysical behavior is shown to stem from improper regularization of the self-energy of dissociated ions and polyions and is mitigated by smearing point-like charges across a finite width. The influence of the self-energy on counterion binding is elaborated on by generalizing the Debye-Hückel free energy to polyelectrolytes with variable fractal dimension and stiffness. In the dilute regime, a greater propensity for binding is found for chains with more compact architectures, which in turn reduces the harsh self-repulsions of tightly packed arrangements of charge. In the concentrated regime, the effects of electrostatic correlations weaken due to screening and the extent of binding is governed by a balance of short-ranged interactions and the translational entropy of ions.

Inertial migration of a rigid sphere in plane Poiseuille flow as a test of dissipative particle dynamics simulations.

After reviewing and organizing the literature on the problem of inertial cross-stream migration of rigid spheres in various geometries including tubes and channels, we use Dissipative Particle Dynamics (DPD) simulations to study the simplest case of migration of a single neutrally or non-neutrally buoyant sphere with diameter 20% of the gap in plane Poiseuille flow and assess the potential and limitations of DPD simulations for this and similar problems. We find that the neutrally buoyant sphere lags by up to 6% behind the surrounding fluid and is focused at a position around 50% of the distance between the channel center and the wall. With Re increasing from around 100 to 500, the sphere migrates closer to the channel center. With flow driven by gravity, a much denser non-neutrally buoyant sphere leads the surrounding fluid and is focused at a position closer to the wall, around 60% the distance from the channel center to the wall, in qualitative agreement with previous work. The lower values of the Schmidt number Sc in DPD simulations relative to real fluids, due to the relatively large diffusivity of DPD beads, are shown to not significantly affect the consistency of our DPD results with literature results although they make results noisy at low Re (i.e., ≲ 50). However, the increase in Ma and Wi with increasing Re leads to compressible flow effects and in some cases viscoelastic effects at high Re depending on the DPD parameters chosen. Even for optimally chosen parameters, we require Re ≲ 500 to avoid strong compressibility effects. Thus, the relative simplicity of the DPD method for complex fluid flows is offset by the need to control the effects of unphysically high values of other parameters, such as Ma and Wi, which seriously limits the range of conditions under which DPD simulations give valid results in fluid transport problems.