Interfacial properties of polymeric complex coacervates from simulation and theory Free

Oppositely charged polyelectrolytes can undergo an associative liquid-liquid phase separation in salt solution in a process known as complex coacervation.^1–4 This phase behavior has been widely exploited in both industry and research due to its sensitivity to the ionic environment; for example, the food and personal care industries have long used coacervates as encapsulants and viscosity modifiers,^5,6 and research has considered coacervation as a route to controlled drug release vehicles,^7–9 as tissue engineering scaffolds,^10,11 and in synthetic underwater adhesives.^12,13 More recently, coacervation-driven block copolymer assembly has been used for multimodal probes and reactive sensors.¹⁴ This diverse range of applications stems from the ability to readily tune the strength of coacervation by changing the chemical and physical properties of the coacervate-forming polyelectrolytes and their environment.

Building on longstanding interest in coacervate materials,^2,3 there has been a resurgent effort to characterize their equilibrium and dynamic behavior on a fundamental level. These include extensive experimental efforts; for example, the phase behavior of coacervation has been explored for a number of different polyelectrolyte systems, as a function of salt concentration,^15–19 salt valency and identity,²⁰ temperature,¹⁶ chirality,^21,22 pH,^16,20 molecular weight,^15,17 monomer sequence,²³ and polymer architecture.²⁴ These largely rely on turbidity as a proxy for phase separation since while phase separation is macroscopic, initial droplets are resistant to coarsening due to a low surface tension.^25,26 Further investigations have characterized a number of physical properties in these same systems, such as chain conformation,^27,28 bulk coacervation thermodynamics,^29,30 bulk rheology,^31,32 and surface tension.^25,33,34 Coacervates are promising because they are affected by a large number of chemical and physical parameters; however, this also highlights the need for good theoretical models to aid the exploration of such a large parameter space.

Bulk coacervate studies have gone hand in hand with efforts to understand oppositely charged block copolyelectrolytes that assemble via coacervation. Building on longstanding efforts to develop coacervate-core micelles for drug delivery,^7–9 both diblock and triblock coacervates self-assemble into a variety of (for example) BCC and hexagonally packed cylinder ordered morphologies.^35–37 Initial phase behavior has been mapped for a few simple systems;³⁵ however, there has yet to be a full physical understanding of coacervate-driven assembly akin to the developed theory of χ-driven block copolymer solution assembly.^38–41 This is in part due to the rapidly evolving state of our understanding of bulk coacervates,⁴ but also due to a still-nascent understanding of the interface between coacervate and supernatant. Indeed, existing theory on standard χ-driven assembly relies on classic work on polymer-polymer interfaces;^40–43 however, there have so far been only a few attempts to understand this interface in coacervates.^33,44

Most recent attempts to understand coacervate interfaces are built around the traditional theoretical understanding of complex coacervation, known as Voorn-Overbeek theory.^15,33,45,46 The free energy of the system is given by

where α is an electrostatic interaction parameter, σ_i is the charge density of species i (polycation, polyanion, cation, anion, or water), and ϕ_i is the volume fraction of species i. The first term on the right side of the equation accounts for the mixing entropy for each species. Electrostatic interactions are captured in the second term using the results of Debye-Hückel theory,⁴⁷ and the third term accounts for short-range interactions via the Flory χ-parameter.⁴⁸ This free energy has been used to determine phase diagrams,⁴⁵ which can match some experimental phase diagram measurements.¹⁵ Likewise, experimental measurements of interfacial tension exhibit scaling behaviors predicted from Voorn-Overbeek-based calculations by Qin et al.³³ Despite this success, it is widely understood that Voorn-Overbeek theory does not adequately capture the physics of coacervation.^3,4,49–52 This has long been appreciated because this theory (1) uses Debye-Hückel theory that is only accurate in the low-salt limit⁵³ and (2) neglects polymer connectivity,^4,49–51 instead of treating salt charges and polymer charges equivalently. Recent work by the authors and others has demonstrated using simulation and theory that these effects are important^{50,51,54–56} and that matching to Voorn-Overbeek is due to a fortuitous cancellation of errors due to the neglect of both the charge connectivity in the polyelectrolytes and the finite size of the charges.^18,49 There have been a number of efforts to systematically incorporate these effects using more sophisticated field theories, such as the random phase approximation (RPA) and field theoretic simulations.^{44,51,54–60} However, these techniques still have limitations; recent particle-based simulation results have shown that local charge correlations play a primary role in coacervation.¹⁸ This limitation can be overcome by either (1) using fully fluctuating fields using, e.g., complex Langevin methods^44,59,60 or (2) using theory or simulations in order to resolve local charge correlations and thus inform field theories.⁶¹ The first approach has resulted in investigations using fully fluctuating field theories,^44,59 and in this article, we will pursue the second approach.

Recent studies by the authors have developed a new theoretical approach to coacervation that is inspired by the concept of counterion release, which is a widely used concept in the pair complexes of dilute polyelectrolytes.^62–64 Here, counterions condense onto polyelectrolytes to neutralize highly charge-dense polyelectrolytes at the cost of their translational entropy.^65–67 When complexing polyelectrolytes condense onto one another, the counterions are released increasing their translational entropy, driving complex formation.^62,63 This concept has been included into a number of coacervate models,^68–70 for example, via “doping” of solid complexes^19,30 or in combination with Voorn-Overbeek arguments.⁵⁰ Recently, the authors have developed a transfer matrix (TM) theory that formalizes these ideas by mapping complex coacervation onto a one-dimensional adsorption model.^71,72 The substrate is a test polyelectrolyte chain, which can adsorb the various other oppositely charged species. The statistics of this adsorption can be parameterized by Monte Carlo (MC) simulation to yield a free energy expression^71,72

Here, a transfer matrix M was developed that captures the Boltzmann factors of the potential adsorption states along a test polyelectrolyte of length N. The vector $ψ→1$ sets the initial Boltzmann factors for the first monomer, and $ψ→0$ is a vector of ones. Additional parameters include κ, which sets the magnitude of a phenomenological third-order term that captures excluded volume, and Λ, which reflects the difference in excluded volume between the polymer and monomer. This theory can nearly quantitatively match simulation results with reasonable choices for fitting parameters⁷¹ and shows consistency between the local adsorption statistics (i.e., the local charge correlations) as well as the phase behavior.⁷¹ This model can be straightforwardly adapted to reflect changes in molecular features such as charge spacing, polyelectrolyte rigidity and architecture, solvent dielectric constant, and salt valency.^71,72

In this article, we start by describing a range of theoretical and computational methods (Sec. II), demonstrating how we use both Monte Carlo (MC) simulation and our recent transfer matrix (TM) approach to inform self-consistent field theory (MC-SCFT⁶¹ and TM-SCFT). In Sec. III A, we show that both TM-SCFT and MC-SCFT can quantitatively match both coacervate bulk thermodynamics and the structure of the interface seen in particle-based molecular dynamics (MD) simulations, demonstrating the efficacy of field theory techniques developed by the authors.⁶¹ We subsequently demonstrate in Sec. III B that the two SCFT-based methods, informed by theory or simulation, exhibit the same interfacial thermodynamics; this includes the surface tension, interfacial width, and the surface excess of salt ions. In Secs. III C–III D, we modify both MC-SCFT and TM-SCFT to include a neutral polymer and to understand interfaces similar to those found in block copolyelectrolytes. We show that increasing neutral polymer density drastically increases the immiscible region of the phase diagram due to the interplay between the excluded volumes of the polymeric species. Indeed, this can lead to large regions of small, but positive, surface tension; this understanding is crucial to understand the thermodynamics of coacervate-driven assembly.

We use a combination of theoretical and computational methods to study complex coacervate interfaces, in part to explore their limitations and capabilities. Specifically, the MD simulations require the fewest assumptions because all chains and salt ions are explicitly represented; however, this comes at a significant computational cost. SCFT is far more efficient, but comes with the trade-off that all chains are treated at a mean-field level. However, we use methods (both TM and MC) that provide information about the local charge correlations that are crucial to understanding coacervate thermodynamics.^49,73 TM and MC provide different routes to this correlation information, with MC accounting for charge correlations directly from simulation^18,24,61 while TM reflects a theoretical model developed by the authors.^71,72 We outline all three methods here, leaving details for  Appendixes A and  B.

There are approximations present in all of our methods. In particular, all methods consider the same restricted primitive model (RPM) representation of all the charges in the system.⁷⁴ The RPM is a standard approximation that removes atomistic detail, leaving charged species and their corresponding hydration shells as charged beads in an implicit solvent medium.^73,74 The charged beads represent both the small-molecule ions as individual beads and polyelectrolytes as connected chains of beads. By neglecting atomistic detail and solvent degrees of freedom, we are unable to describe specific ion and Hofmeister effects,⁷⁵ as well as any water structure effects that can become important at high polymer concentrations.^76–78 These effects are largely thought to prevent quantitative, but not qualitative matching to experiment, and indeed most comparisons between theory and simulation do not rely on these details.⁴ This model also neglects dielectric effects,^79–85 which also may play a role in the thermodynamics of coacervation. We similarly do not expect this to play anything more than a quantitative role.

We consider a system composed of polyelectrolytes, neutral polymers, salt, and water and use SCFT methods standard in the literature.^86–90 To simplify our SCFT calculation, we assume both local electroneutrality and charge stoichiometry between both polymers and ions; therefore, we only need one polyelectrolyte species representing both polycations and polyanions and one salt species representing both cations and anions. This is a standard assumption common to most coacervate models and is used to simplify experiments;^3,4,45 this can in principle be relaxed.⁷¹ For this paper, the degree of polymerization is N = 100 and is the same for both the polyelectrolyte and the neutral polymer. Species-specific quantities are denoted with subscripts; A for polyelectrolytes, B for neutral polymers, S for salt, and W for water. We thus define a set of volume fractions ϕ_i(x), where i = A, B, S, W that is a function of a spatial position x.

The homogeneous free energy F of a coacervate system can be written as^51,61

Here, we introduce the per-volume free energy f. The first four terms on the right side account for the translational entropies of all the species. The fifth term is an excess free energy that accounts for all remaining contributions to F and is generally a function of the set of volume fractions ϕ_j. $fEXC{ϕj}$ is the excess free energy, which accounts for the local charge interactions that drive coacervation. We will describe how to obtain explicit forms for values and derivatives of $fEXC{ϕj}$ in Secs. II B and II C as well as  Appendix A (for MC) or  Appendix B (for TM).

The homogeneous free energy F in Eq. (3) corresponds to a Hamiltonian functional $H$⁠,⁹⁰ 

The Hamiltonian is a functional of the position-dependent density fields ϕ_j(r) and auxiliary fields ω_j(r) for all of the species.^90,ρ₀ is the bulk number density. The single chain partition functions Q_A and Q_B, along with the single particle partition functions Q_S and Q_W, are functionals of their respective auxiliary fields iω_j(r).⁹⁰ The first term in the integral is the local excess free energy, to be calculated from MC or TM, and the second term in the integral ensures that the sum of the volume fractions is 1 with a coefficient ζ that is set to be large. In writing this Hamiltonian for inhomogeneous coacervate systems, we introduce an ansatz that a homogeneous value of the excess free energy f_EXC({ϕ_j}) can be used locally at points r, leading to a position-dependent $fEXC({ϕjr})$⁠. This term is thus approximated as a local contribution to the Hamiltonian, based on the observation that coacervation takes place at high salt/polymer concentrations where we expect local charge correlations to dominate.^15,18,52 We do not justify this “local homogeneity” ansatz a priori because electrostatic correlations in MC simulations are longer-range than the field theory grid spacing; we will instead justify a posteriori through comparison to particle-based MD. This theory thus differs from alternative field theoretic formalisms that directly include the Coulomb interaction in the Hamiltonian,^44,59,91 simplifying the SCFT calculation but introducing an assumption that we can use a free energy associated with charge correlations while neglecting how they may be affected by concentration gradients that arise in inhomogeneous polymer systems.^42,90

The goal of SCFT is to find the “saddle point” or sets of fields $ϕj*(r)$ and $ωj*(r)$ where the Hamiltonian is at an extremum,⁹⁰ 

These functional derivatives with respect to the auxiliary field lead to the expression for the density in terms of the propagators,^86–90 

For the polymer components, this becomes^86–90 

For the water and salt components, this becomes^86–90 

Here, we have defined a field W_j = iω_j, and the quantity q_j(r, s; [W_j]) in Eq. (7) is a propagator that represents the probability that a sub-chain of length s has its terminal chain end located at a position r in a field W_j. These are calculated via the diffusion equation, which in one dimension (r → x),^86–90 

We also consider the functional derivatives of the Hamiltonian [in Eq. (5)] with respect to the density fields,^86–90 

The partial derivative is the local excess chemical potential μ_j,EXC(ϕ_i(r)), calculated with respect to the number density ρ_j of species j. This leads to the result^86–90 

In this paper, the value of μ_j,EXC(ϕ_i(r)) is determined from either MC or TM using the set of values ϕ_i(r). A self-consistent scheme is used to calculate the saddle point. This consists of calculating the fields W_j(r) using the species densities ϕ_j(r) via Eq. (11). These fields are then used to calculate the propagator for the polymer species $qjx,s;[Wj]$ in Eq. (9) and subsequently the species densities ϕ_j(r) using Eqs. (7) and (8). These densities can once more be used to calculate W_j(r), and this process repeats until all of the fields have converged to their saddle point values $Wj*(r)$ and $ϕj*(r)$⁠. For this paper, we use the value ζ = 1000k_BT.

We can use simulation to obtain f_EXC so that detailed charge correlation information is incorporated into our SCFT calculations. In this work, we use a standard Monte Carlo (MC) simulation based on the restricted primitive model (RPM) of charged polyelectrolyte systems.^18,23,24,61 The RPM considers all charges to be hard spheres in an implicit solvent, which interact through a Coulomb potential. Polyelectrolyte charges are connected via bond and angle potentials. In this article, we call the hybrid method MC-SCFT; related methods have been used extensively in prior work by the authors,^23,24,61 so we relegate the model in detail to  Appendix A. The key feature of these simulations is that standard Widom insertion methods can be used to determine the excess chemical potential,^61,92,93

Thermodynamic integration of the excess chemical potential enables calculation of f_EXC, which in turn serves as an input to SCFT.

We can also use theory to obtain f_EXC, and in particular, we use the transfer-matrix (TM) theory developed by the authors.^71,72 This model maps the local charge environment around a test chain to a one-dimensional adsorption model.^71,72 This is motivated by the correlations observed in simulation, which exhibit prominent neighbor peaks in an otherwise highly charge screened system.^18,23 Monomer “sites” on the test chain can adsorb oppositely charged species, including both the oppositely charged polyelectrolyte or the counterion and occasionally may remain without an adsorbed species. A transfer matrix M can be written to account for the Boltzmann weight associated with the local environment around each monomer and can be used to calculate a grand canonical partition function that is related to the interaction free energy.⁷¹ With the inclusion of a phenomenological cubic term that accounts for excluded volume, the analytical expression for f_EXC is^71,72

Here, the function Θ is^71,72

We have used the analytical expression for the interaction free energy, where Θ is the largest eigenvalue of M and the quantities A₀, B₀, and $ϵ̃$ can be parameterized directly from molecular simulation.⁷¹ We use this form for this paper, but more generally it is possible to numerically calculate the interaction free energy via Eq. (2) if M is no longer analytically tractable.⁷² See  Appendix B for a more detailed discussion.

Molecular dynamics (MD) simulations were performed on a system consisting of n_A+ polycations and n_A− polyanions, modeled as connected monovalently charged beads with diameter σ. Both polymers have the same degree of polymerization, N. n_S+ cations and n_S− anions are modeled as monovalently charged beads, also with a diameter σ. Water is included as an implicit solvent with a relative dielectric constant, ϵ_r = 78.5. These simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) with a Langevin thermostat.⁹⁴ See Fig. 1 for a schematic of our model. The overall potential energy is given by

U_E models the electrostatic interactions,

where q_i is the charge on bead i, ϵ₀ is the permittivity of free space, and r_ij is the separation between beads i and j. Standard Ewald summation is used to account for the long-range interactions under periodic boundary conditions.⁹² The polymers are bound together with a bond potential U_B,

For this contribution, κ_B is the strength of the bond potential and r₀ is the equilibrium bond distance. In addition to the bond potential, the polymers have a bending potential U_θ,

Here, κ_θ is the strength of the bending potential, θ_{i,i−1,i−2} is the angle between the two bond vectors, and θ₀ is the equilibrium angle between bond vectors. The excluded volume of the beads is modeled by a Lennard-Jones potential U_LJ,

where ϵ_LJ is the depth of the potential well, σ_LJ is the interparticle separation at which U_LJ becomes zero, and r_c is the cutoff distance of this potential.⁹² 

For this simulation, the degree of polymerization is kept at N = 100. We choose parameters consistent with our previous efforts;^18,61 κ_B = 250k_BT, r₀ = 1.05σ, κ_θ = 1.65k_BT, and θ₀ = π. Lennard-Jones parameters to model excluded volume were ϵ_LJ = 10.75k_BT, σ_LJ = σ, and r_c = σ. This mimics the hard-sphere potential in the MC simulations by having the interparticle excluded volume potential at zero until the particles overlap. The bead diameter σ was taken as 4.25 Å.

We use MD simulations and SCFT to capture the physics of this coacervate system, with the latter informed by both MC simulations and TM theory. We can show that all methods, while they provide different levels of resolution and assumptions, exhibit nearly quantitative matching for interfacial properties. Importantly, we motivate our use of MC-SCFT and TM-SCFT by comparing with a completely particle-based simulation. To ensure that we have comparable models, we show that the MC used to inform the SCFT yields the same radial distribution functions when compared to the full-particle MD simulations. Figure 2 shows radial distribution functions from both methods (MC and MD) using n_A+ = n_A− = 6 in a cubic simulation box for ϕ_A = 0.08 and ϕ_A = 0.12 [Figs. 2(a) and 2(b), respectively]. Both MC and MD simulations agree nearly quantitatively. Some differences do occur due to differences in the interaction potentials; the MC uses hard-core interactions that are not possible in MD simulations, where a Lennard-Jones potential is used instead. These disparities are apparent at the contact peak near r/σ = 1, where the differences in these potentials are most pronounced.

The phase behavior of a coacervate with no neutral polymer is shown in Fig. 3(a) (calculated from MC-SCFT), showing the two-phase coacervation regime at low salt ϕ_S and polymer ϕ_A concentrations. Sloped tie lines are observed connecting the large-ϕ_A coacervate phase at a slightly lower ϕ_S than the low-ϕ_A supernatant phase. This is due to the partitioning of salt to the supernatant, which has been demonstrated by the previous computational, theoretical, and experimental work.^{18,30,49,71,72} We also show MD simulation snapshots from phase-separating mixtures along the indicated tie lines, with the orange and blue polymers representing the polycations and polyanions, respectively. This structure becomes more diffuse as salt is increased. This is apparent in Figs. 3(b)–3(g), where we plot the interfacial profiles (ϕ_i versus distance x/σ) that form along these same tie lines [Fig. 3(a)]. These plots include full-particle MD simulations, along with interfacial profiles calculated with both MC-SCFT and TM-SCFT. MD simulations were performed with n_A+ = n_A− = 15 in a rectangular prism simulation box with the x-dimension 3 times the length of the y- and z-dimension. Starting volume fractions for each simulation were taken as the midpoint of the tie lines in Fig. 3(a).

All 3 techniques exhibit nearly quantitative matching in all interfacial profiles shown. This a posteriori justifies the approximations we made in the TM-SCFT and MC-SCFT, which is that we have a separation of length scales between the charge correlations (which are treated as a homogeneous system at each SCFT grid point) and the polymer length scales (which span multiple SCFT grid points).⁶¹ We do note that correlations in Fig. 2 have a similar range to the grid spacing used in our simulations (0.75σ), which is set relatively small to resolve the interfacial structure. However, this correlation length scale seems to still be sufficiently small compared to the length scales of variation in ϕ_A and ϕ_S such that square gradient corrections to f_EXC are unnecessary to match SCFT with MD. The increasingly diffuse interface at higher salt concentrations poses a practical challenge for the full-particle MD, where it is difficult to obtain a good average of the interface; this is why we focus on tie lines far from the critical point. We thus demonstrate that MC-SCFT and TM-SCFT capture the correct interfacial profile, motivating our use of these field-based methods to calculate interfacial properties such as interfacial tension, interfacial width, and the surface excess of salt.

We can use the results of both MC-SCFT and TM-SCFT to determine the interfacial tension of a coacervate. We use the expression^95–99 

We have used tildes to denote normalization of energy scales by k_BT and length scales by σ. The second and third terms in the integrand are square gradient expressions for the polyelectrolyte and the neutral polymer, approximated using the result from the Random Phase Approximation (RPA).^98,99 The first term in the integrand is the free energy of a homogeneous system with the volume fractions at a point x relative to the bulk free energy of the system, given by the expression^96,97

Here, $μi∞$ is the chemical potential of species i in the bulk phase.

We show a calculation of this interfacial tension $γ̃$ for a coacervate without the presence of a neutral polymer, shown in Fig. 4(a). Consistent with the previous experimental and theoretical literature,^{25,26,33,34,44} we plot the decrease in $γ̃$ as a function of salt concentration, where we choose the supernatant salt concentration $ϕSβ$⁠. We will typically denote phases by a superscript α (coacervate) and β (supernatant). Indeed, we demonstrate that both MC-SCFT and TM-SCFT exhibit nearly identical values of $γ̃$⁠, with only a small horizontal offset of the $γ̃$ versus $ϕSβ$ curve. We can also show that our results exhibit the scaling law described by Qin et al.,³³ for the interfacial tension $γ̃/γ̃0∼(1−ϕSβ/ϕScrit)3/2$ as $ϕSβ$ approaches its critical point value $ϕScrit$ (Fig. 5). Here, interfacial tension is normalized by its zero-salt value, $γ̃0$⁠. Qin demonstrated matching between this prediction and with experimental data^25,26,33,44 despite the use of Voorn-Overbeek theory in the derivation of the scaling behavior. Our results are consistent with the Qin prediction, showing a 3/2 slope in a log-log plot in Fig. 5 for both TM-SCFT and MC-SCFT; our results are also consistent with the experimental results from the literature,^25,26,33,44 which are also included in Fig. 5 as open symbols. We note that these results are based on mean-field SCFT models; a different power law will likely be observed as the critical point is approached, which would require beyond-mean-field polymer field theories.

We consider other interfacial properties, such as the interfacial width D. D is defined using the interfacial profile for the polyelectrolyte,⁹⁵ 

Here, the derivative in the denominator is taken at the midpoint concentration between the polyelectrolyte concentrations in the coacervate phase $ϕAα$ and the supernatant phase $ϕAβ$⁠. We plot D/σ as a function of $ϕSβ$ for both MC-SCFT and TM-SCFT in Fig. 4(b). Similar to Fig. 4(a), we again see mainly small differences in this quantity between the two methods. The exception to this comes at large values of $ϕSβ$⁠, which is close to the critical point such that D/σ becomes large. The curves do not become large at the same rate, consistent with Fig. 4(a) which shows that the critical value of $ϕSβ$ appears to be shifted. The physical observation that D increases is consistent with traditional polymer theory,^42,43,98,100 where a decrease in the driving force for phase separation broadens the interface. In non-charged systems, this is usually controlled by temperature; however we are instead weakening the electrostatic interaction by using salt concentration $ϕSβ$ instead.

This system has effectively 3 species, the polyelectrolytes, salt ions, and water; thus, we can calculate the surface excess of salt at the interface as the overall salt concentration is increased. The surface excess of any species i can be calculated via the relationship^96,97

Here, D is defined in Eq. (22) and D^α is the size of the phase α. This definition, however, depends on the choice of interface position; a related value is independent of this choice,^96,97

This value of $Γi(j)$ describes the surface excess of i, using a reference of j ≠ i such that it does not depend on the choice of j ≠ i. This physically represents the deviation in excess ϕ_i [captured by the first term in Eq. (24)] away from any excess ϕ_j calculated from the normalized interfacial structure in the other species j [the second term in Eq. (24)]. We plot $ΓS(A)$⁠, the surface excess of salt, as a function of $ϕSβ$ in Fig. 4(c).^96,97 We once more observe only small differences between MC-SCFT and TM-SCFT and show that the surface excess approaches $ΓS(A)=0$ at the limits of $ϕSβ$ (where there is no salt) and at the critical value of $ϕSβ$ where the system becomes homogeneous. At intermediate values, the surface excess reaches a maximum, indicating that salt preferentially partitions to the interface. The surface excess salt remains quite small, however, with $ΓS(A)∼7×10−4$ indicating that any deviations of ϕ_S from what would be expected from the ϕ_A profile are between one or two orders of magnitude smaller than the overall variations in ϕ_S.

Interfacial tension plays a significant role in the classical theoretical work on the self-assembly of block copolymers.^{40–43,100,101} The form of the interfacial free energy follows that of a simple polymer-polymer interface and is driven by a positive short-range Flory χ parameter.^42,43,90,98 Self-assembled structures formed due to coacervation are driven by long-range electrostatic attractions that are not well represented by χ;¹⁰² instead, we use the analogous interface between a coacervate and a neutral polymer to determine the interfacial tensions that would arise in the coacervate-driven self-assembly. Furthermore, the presence of a neutral polymer has implications for biological coacervates, which are often in crowded macromolecular environments.^103,104 Recent reports from experiment demonstrate that “crowding” can have a significant impact on the formation of bio-inspired coacervates.¹⁰⁴ 

First, we show that the addition of a neutral polymer species significantly affects the phase diagram depending upon the concentration of the neutral polymer [Figs. 6(a) and 6(b)]. To plot the effect of the neutral polymer, we note that it primarily partitions into the supernatant phase. We thus adjust the neutral polymer concentration by fixing its value in this phase $ϕBβ$⁠. The corresponding neutral polymer concentration in the coacervate phase $ϕBα$ is typically negligible, except near the critical point. We thus emphasize that Fig. 6 is the projection of a three-dimensional phase diagram onto the ϕ_S versus ϕ_A plane.

The inclusion of the neutral polymer results in a large extension of the two-phase region at high ϕ_S and low ϕ_A, which extends further to the right as the concentration $ϕBβ$ increases. This trend is observed in both MC-SCFT and TM-SCFT. We interpret this trend as due to the matching of pressure between the supernatant and coacervate phases; the coacervate phase in the absence of a neutral polymer is limited to low ϕ_A and ϕ_S due in part to the high pressure due to excluded volume in the polymer-dense coacervate.¹⁸ This phase would expand except this would dilute the combinatoric entropy of polymer-polymer interactions that primarily drives phase separation.^71,72 With the addition of the neutral polymer, this polymer partitions to the supernatant and applies a counterpressure. When ϕ_A ∼ ϕ_B, the pressures of the two phases balance essentially removing the pressure-based driving force for miscibility; this results in an extended two-phase region. To demonstrate that the changes in phase behavior are driven primarily by the excluded volume of the neutral polymer, we remove the ϕ_B contribution to the excluded volume in the TM calculation [we remove ϕ_B from the second term of Eq. (B5)]. The resulting phase diagram is plotted in Fig. 6(c) and exhibits only minor changes in the binodal as $ϕBβ$ is increased, contrasting significantly with the phase diagrams seen in Figs. 6(a) and 6(b).

These observations are further supported by the examination of the salt partitioning $λ=ϕSα/ϕSβ$⁠, where $ϕSα$ is the volume fraction of salt in the coacervate phase and $ϕSβ$ is the volume fraction of salt in the supernatant phase. This is plotted as a function of $ϕSβ$ in Figs. 6(d) and 6(e), corresponding to parts a and b directly above. In the limit of no neutral polymer, the salt partitions such that the salt concentration in the supernatant is greater than that in the coacervate phase (λ < 1). This is consistent with the previous experiment and simulation,^18,71,72,105 and we attribute this to the excluded volume of the polymer in the coacervate phase. Upon addition of a neutral polymer, the salt partitioning behavior is altered, in different ways depending on the value of $ϕSβ$⁠. At low $ϕSβ$⁠, salt partitioning becomes increasingly even between the two phases (λ → 1). This is likely due to the neutral polymer in the supernatant, which begins to equalize the excluded volume in the two phases as its concentration increases. This occurs below the critical salt concentration of the phase diagram in the absence of a neutral polymer; however at high $ϕSβ$⁠, the opposite trend is observed. Here, larger values of $ϕBβ$ lead to smaller values of λ. We attribute this “switch” to the large extent that the coacervate phase is increasingly pushed to large ϕ_A and ϕ_S, where excluded volume effects become enhanced. We also plot λ versus $ϕSβ$ in Fig. 6(d) when excluded volume of the neutral polymer is removed from the TM calculation. Similar to Fig. 6(c), we find that only minor changes are observed when $ϕBβ$ is increased, consistent with the argument that the excluded volume of the neutral polymer is what drives the marked differences seen in Figs. 6(d) and 6(e).

Our excluded volume argument is consistent with the differences between TM-SCFT and MC-SCFT, particularly with respect to the value of λ. While some aspects of the qualitative trends are similar, there are some differences, namely, the lack of a critical point in the region tested for TM-SCFT and the different magnitudes for the plots of λ. We attribute this to the phenomenological argument for the excluded volume term in Eq. (2), which is chosen to approximate the free energy associated with packing the species.^71,72 We also note that Fig. 6(d) becomes non-monotonic at large values of $ϕSβ$⁠, with an increase in partitioning with increasing $ϕSβ$ at intermediate salt concentrations. This does not occur for the TM-SCFT predictions in Fig. 6(e). We note that this non-monotonicity occurs at a high-ϕ_A and high-ϕ_S region of the phase diagram in Figs. 6(a) and 6(b), where there are also different shapes in the binodal curves. We once more attribute these differences to the phenomenological excluded volume term in TM-SCFT, which should become significant in this region.

We choose a point in the phase diagrams in Fig. 6 to demonstrate that these large shifts in the phase diagram are indeed observed in the full MD simulation. We selected a value of $ϕBβ=0.03$ and average concentrations of ⟨ϕ_A⟩ = 0.025 and ⟨ϕ_S⟩ = 0.12 for all simulation/theory methods (MD, TM-SCFT, MC-SCFT). We plot the resulting volume fraction profiles in Fig. 7, demonstrating that all three methods exhibit near-quantitative matching. Indeed, this choice of salt and polymer concentrations will not phase separate without the presence of the neutral polymer. Snapshots at the bottom of Fig. 7 show the MD simulation corresponding to the interfacial profile. Both show the same system; however, the top snapshot includes the ions, and the bottom snapshot removes the ions to show the polymer components more clearly. Consistent with the expectations from the SCFT calculations, the neutral polymer (green) is almost completely in the supernatant phase and the polycation and polyanion (blue and orange) are completely in the coacervate phase.

We use both MC-SCFT and TM-SCFT to determine how the interfacial tension $γ̃$ changes as $ϕBβ$ is increased. These results correspond to the phase diagrams in Fig. 6 and are shown in Figs. 8(a) and 8(d) for the two different methods. For both sets of results, the interfacial tension monotonically increases with $ϕBβ$⁠, including both above and below the critical $ϕSβ$ at ϕ_B = 0. Interestingly, for small quantities of $ϕBβ$⁠, the surface tension remains very close to $γ̃≈0$ for a large range of $ϕSβ$⁠. This trend is also apparent in the interfacial width D [Figs. 8(b) and 8(e)], which exhibits significantly broader interfaces at small values of $ϕBβ$⁠. After the initial addition of $ϕBβ$⁠, the width of the interface begins to decrease to what appears to be a limiting value of D at high $ϕBβ$⁠. Finally, an increase in $ϕBβ$ drastically increases the surface excess of salt $ΓS(A)$ [Figs. 8(c) and 8(f)], which we attribute to the increase of excluded volume in both the coacervate and supernatant phases. In this case, the salt preferentially partitions to the interface where excluded volume is minimized. This is directly observed in Fig. 7, where the profile shows a significant increase of salt at the interface for all methods (MD, MC-SCFT, and TM-SCFT). The calculated value of $ΓS(A)≈3×10−3$ is similar to the magnitude of the positive deviation in ϕ_S at the interface.

For all interfacial properties in Fig. 8, MC-SCFT and TM-SCFT both capture the same qualitative trends. However, there are differences connected to the lack of critical points in TM-SCFT seen in Fig. 6. This manifests as a decrease in the interfacial tension $γ̃$ and interfacial excess of salt $ΓS(A)$ to 0, along with a sharp increase in D for the MC-SCFT results. Similar to the phase diagram, we primarily attribute these differences to our approximate treatment of the excluded volume in the TM theory.^71,72

In this manuscript, we show that we can use a range of computational and theoretical routes to understand the interfacial properties of polymeric complex coacervates. We can directly simulate interfaces using molecular dynamics, using a restricted primitive model representation that considers all polymeric and salt species as explicit particles. We match these results to SCFT predictions, which capture charge correlations using either a Monte-Carlo based method (MC-SCFT) or a transfer matrix theory of coacervation (TM-SCFT); we demonstrate nearly quantitative matching among the interfacial profiles of all three methods, as well as matching between the pair correlation functions in the MC and MD simulations. We subsequently use both MC-SCFT and TM-SCFT to probe the interfacial thermodynamics of coacervates, demonstrating trends in the interfacial tension γ consistent with previous observations in the literature.^25,26 Similarly, we provide predictions for the interfacial width D and surface excess of salt $ΓS(A)$⁠.

We subsequently consider the effect of adding a neutral polymer due to its relevancy for the interface of self-assembling coacervate block copolymers. The coacervation phase behavior changes drastically with the neutral polymer partitioning strongly to the supernatant phase and “balancing” the excluded volume pressure of the coacervate phase. This results in a large increase in the region of coacervation, seen in both MC-SCFT and TM-SCFT calculations. This change in the phase behavior is commensurate with large changes in the interfacial properties, with a large region of low-γ observed upon initially including the neutral polymer. There are corresponding changes in interfacial width D and significant increases in the surface-excess salt $ΓS(A)$ that are due to the presence of significant excluded volume in both phases.

Nearly quantitative agreement with full-particle MD, as well as field theories informed by both simulation and theory, provides insight into the important features needed to understand the physics of coacervation. First, the separation of length scales between the polymer conformation and the charge correlations is reflected in the SCFT calculations and does not seem to adversely affect the predictive power of the theory. Second, quantitative prediction hinges on the accuracy of the excluded volume model, which plays a large role near the critical point when there is a neutral polymer present; small differences in this term lead to large changes in the observed phase behavior. Finally, the complexity of the phase behavior when the neutral polymer is present poses a challenge for developing intuitive scaling arguments; in contrast to χ-driven phase separation,^42,43,98 it will likely be difficult to distill quantities like γ down to compact scaling expressions.

While the agreement between the different methods is excellent, we note the limitations of our current approach. In particular, all of the polymers considered possessed a high charge density. This is crucial for the transfer matrix theory, which assumes that most of the charges along the chain are matched with a condensed counterion or an oppositely charged chain.⁷¹ It may be possible to extend the theory to low charge densities; however in these limits, we approach regions where field theoretic approaches may be more apt.^{33,44,51,54–58,106} Nevertheless, we use realistic parameters in our model that are roughly consistent with standard polyelectrolytes used in complex coacervates and have exhibited qualitative matching to experiment in a number of circumstances.^18,23,24 We further note that all methods used in this paper (MC-SCFT, TM-SCFT, and MD) in principle have challenges in describing systems where both ϕ_S and ϕ_A are small. For MD and MC, we can easily go to low salt concentrations but do not see a large enough simulation box to explicitly capture the polymer concentrations seen on the low-ϕ_A branch of the binodal curve. We instead rely on the extrapolation of f_EXC from finite concentrations of A to this limit. Similarly, TM does not consider changes in polymer conformations that are known to play a role in polyelectrolyte thermodynamics at low polymer concentrations.^{56,68,106,107} These may limit the accurate determination of the low-ϕ_A branch of the binodal curve, but we expect this to have little quantitative effect on our predictions of the high-ϕ_A branch.

We anticipate that this work will inform the further development of an understanding of coacervates, in particular, the self-assembly of coacervate-based block copolymers. Models of block copolymer self-assembly incorporate surface free energies via interfacial tension γ,⁴⁰ and we have provided an initial picture of how this term would behave in the case of block copolyelectrolyte solution self-assembly.

This material is based upon work supported by the National Science Foundation under Grant No. DMR-1654158. The authors also acknowledge helpful conversations with Sarah L. Perry.

MC simulations were performed on systems with the same species as in the MD simulations, shown in Sec. II D. However, neutral polymer chains, n_B with the same degree of polymerization as the polyelectrolytes was added in order to model its effect on the phase separation and interfacial profiles of complex coacervates. The total potential energy in these simulations includes an electrostatic potential, a bonding potential, a bending potential, and a hard-sphere potential,

where U_E and U_θ are the same potentials, as described in Sec. II D. Potentials used in the MC simulations and the MD simulations differ only by how the bonding and excluded volume potentials are modeled. Bond potentials in MC simulations are described as

where the asterisk over the summation indicates that only bonded species are considered. This describes the bonds between polymer beads as a square well potential, whereas the MD simulations described bonds as a Gaussian potential. Excluded volume is described as a hard sphere potential U_HS,

In comparison to MD simulations, the hard sphere potential is stronger than the Lennard-Jones potential used to model excluded volume. Differences seen in the radial distribution functions shown in Fig. 2 are attributed to differences in how the bond and excluded volume potentials are modeled in MC and MD.

In this manuscript, in addition to being used to calculate radial distribution functions in Fig. 2, MC simulations were used to determine the excess free energy landscape as a function of salt, polyelectrolyte, and neutral polymer concentration. Widom insertion was performed for each species, i, to calculate excess chemical potentials, μ_EXC,i. This was done using joint insertion of charge neutral pairs and chain-end monomer addition to calculate per-monomer μ_EXC,i for the polymeric species. μ_EXC,i allows calculation of the excess free energy from a reference state $ϕS0$⁠, $ϕA0$⁠, and $ϕB0$⁠,

These f_EXC values can be fit with a third-order polynomial using orthogonal distance regression. Simulations to determine chemical potential values in the limit of ϕ_B = 0.00 were performed with n_A+ = n_A− = 6, and simulations to determine chemical potential values with the neutral polymer used n_B = 12.

Previous work by the authors has demonstrated that a molecularly informed transfer matrix theory of complex coacervation nearly quantitatively matches both phase diagrams and the molecular structure for a number of systems with different molecular features.^71,72 An overview of the theory will be included here, but detailed descriptions can be found in previous work.^71,72

Polymers considered in this manuscript are in the large charge density limit, so a majority of monomers can be assumed to be adjacent to a condensed species of opposite charge.^62,65,71 Both the entropic release of counterions and the combinatoric entropy of many chain-chain interactions drive complex coacervation.^67,71,72 The entropic contributions from these effects depend on the concentration of polyelectrolytes and salt ions. These condensed states can be mapped to an adsorption model by considering each monomer as an adsorption site.⁷¹ Possible states for each adsorption site are a counterion, C, a different oppositely charged polyelectrolyte, P′, the oppositely charged polyelectrolyte continuing its run, P, or the monomer can be unpaired, 0. Transfer matrix formulation allows for calculation of the grand canonical partition function,

where ψ is a column vector of ones and M is a matrix of Boltzmann weight for each adsorption state. Raising the matrix M to the Nth power allows enumeration of the Boltzmann weight of each monomer adsorption state to the partition function.

The matrix M is given by

Each matrix entry is the weight associated due to an adsorbed species at monomer s_i given the adsorbed species at monomer s_i−1. Adjacent pairs of monomers are assigned to factors $A=eμ̃S$⁠, $B=eμ̃P$⁠, and $D=e−ϵ̃$⁠, where the tilde denotes normalization by k_BT. $μ̃S$ and $μ̃P$ are the chemical potentials for salt ions and polyelectrolyte monomers. $ϵ̃$ is a fitting parameter which captures the electrostatic driving force for oppositely charged species to condense onto the chain. For this work, we assume that it is constant for all values of ϕ_A, ϕ_B, and ϕ_S; however, more generally, it may be important to consider variations in $ϵ̃$ due to changes in the local electrostatic environment, especially in inhomogeneous polymer systems such as those described in this paper.^108–110 The parameter E is the single-monomer partition function describing the confinement of adsorbed monomers after the initial monomer is adsorbed to the chain. The transfer matrix M can be multiplied by an arbitrary constant without affecting the thermodynamic properties, so we are free to choose E = 1.

The partition function, Eq. (B1), is analytically solvable by assuming the largest eigenvalue dominates,

A free energy of interaction, F_int, can be calculated as

where A and B can be determined using simple expressions for the chemical potentials $μ̃S=μ̃S0+lnϕS=lnA0ϕS$ and $μ̃P=μ̃P0+lnϕA=lnB0ϕA$⁠. $μ̃S0$ and $μ̃P0$ are reference chemical potentials for the salt and polyelectrolyte, and $A0=expμ̃S0$ and $B0=expμ̃P0$⁠.

The excess free energy of the system, F_EXC, can be calculated using Eq. (B4) with a phenomenological expression capturing the excluded volume of the non-water species,^71,72

where κ determines the strength of the excluded volume interaction and Λ = 0.6875 accounts for the smaller excluded volume of the polymer relative to the salt ions. Substitution of Eq. (B4) and the chemical potential expressions into Eq. (B5) gives the excess free energy as

where

For the species in the system, we can write the chemical potentials as

Here we have defined

These expressions for μ_j,EXC can be used in Eq. (11) to incorporate the transfer matrix theory into the SCFT calculation. Parameters used in this manuscript are A₀ = 20.5, B₀ = 12.2, ϵ = 0.0, and κ = 19.0, consistent with previous work.⁷¹