Mixing protocols determine liquid–liquid phase separation dynamics in polyelectrolyte complex coacervation

Domain: physics
Published: 2026-01-10
Editorial outcome: major_revision

Abstract
Polyelectrolyte complex coacervation underpins many critical biological processes, yet how different initial mixing protocols determine its liquid-liquid phase separation (LLPS) dynamics remains unclear. Using molecular dynamics simulations, we show that when polycations and polyanions are initially randomly mixed, coacervate domain growth exhibits transient t1/2 scaling, driven by polymer network formation. This phase is followed by either t1 scaling due to hydrodynamic pumping or t1/3 scaling from droplet coarsening, depending on the initial mixing degree. Conversely, starting with spatially separated domains of polycations and polyanions - mimicking LLPS in certain marine organisms - leads to rapid coacervate formation, with early-stage growth following distinct t2/3 scaling due to strong electrostatic attraction, followed by continued growth via polymer accumulation. Both protocols yield significantly faster dynamics than systems initialized with preformed polyion pairs, which exhibit classical t1/3 scaling characteristic of droplet coarsening. These findings highlight the profound impact of initial conditions on LLPS dynamics in polyelectrolyte systems.

Introduction
Polyelectrolyte complex coacervation refers to the liquid-liquid phase separation (LLPS) driven by electrostatic association of oppositely charged polyelectrolytes in an aqueous environment, resulting in the formation of a polymer-rich coacervate phase coexisting with a dilute phase1,2. In recent years, this LLPS process has received significant attention due to its relevance to various critical biological phenomena, including the formation of membraneless organelles3–6, the progression of age-related diseases7,8, and underwater adhesion of marine organisms9–11. While substantial advances have been made in understanding the thermodynamics of LLPS over the past decade, the dynamics of LLPS—including the kinetic pathways and rates of phase separation—remains relatively underexplored. LLPS dynamics play a crucial role in the spatiotemporal regulation of intracellular organelles12, and can significantly impact the life survival of marine species in harsh undersea conditions9–11.

Recently, we studied LLPS dynamics in polyelectrolyte complex coacervation starting from an initial state of preformed polyion pairs13. For symmetric mixtures, the growth of the coacervated domains was shown to follow the classical droplet coarsening mechanism14, characterized by a power-law growth of the coacervated domain as R~tβ, with β=1/3. We shall refer to this pathway as the “thermodynamic pathway” on consideration that, thermodynamically, coacervation can be considered a two-step process: the formation of polyion pairs, and the condensation of polyion pairs into a macroscopic coacervate phase15–18. Other factors, such as polyelectrolyte charge asymmetry13,19, intracellular activities20,21, and coupling to elastic deformation of the cellular matrix22–24 can further suppress the LLPS dynamics. Consequently, the thermodynamic pathway often results in relatively slow coacervation dynamics.

While the thermodynamic pathway can describe the relatively slow coarsening dynamics of coacervate droplets starting from dilute solutions, in reality, LLPS often starts with dense polyelectrolyte (PE) mixtures and exhibits significantly faster phase separation dynamics. For example, many experiments trigger LLPS by applying environmental cues such as lights22 and pH25 to concentrated, well-mixed polyelectrolyte solutions, which results in complete phase separation in a few minutes22,25,26. This relatively rapid phase separation can be conceptually understood by the facilitated contacts between oppositely charged polymers during mixing27–29, while a quantitative examination of the underlying kinetic pathway and scaling law is not fully explored. Recent work by Yuan and Tanaka30 demonstrated that this pathway involves the formation of a network structure, with a t1/2 scaling in domain growth, characteristic of viscoelastic phase separation31. In marine organisms such as Sandcastle worms and mussels, even faster LLPS dynamics are observed, whereby these organisms produce proteins that rapidly form adhesive coacervates to attach to wet surfaces32. In this scenario, oppositely charged glue proteins are secreted separately into “tubes” (e.g., byssal threads in mussels) at high concentrations, and LLPS is initiated by fluxing these proteins into a shared region, with the coacervates forming in as little as 30 s32. Beyond biological systems, recent studies found that applying shear flux to a mixture of polycation and polyanion solutions can significantly speed up their coarsening, where a clear phase separation is observed on the order of seconds33. In comparison, without the initial flux, the phase separation of the mixture can only be observed after at least a few minutes27,34–36. Clearly, different initial mixing protocols can significantly influence the phase separation dynamics of dense polyelectrolyte mixtures. However, the fundamental molecular mechanisms driving the dynamics of LLPS remain poorly understood. In particular, the complete morphological pathways under varying initial mixing protocols, as well as their corresponding time evolution, are not yet fully elucidated.

In this work, using molecular dynamics simulation that explicitly accounts for electrostatics and hydrodynamics, we examine the kinetic pathways and the corresponding LLPS dynamics in polyelectrolyte complex coacervation under several representative initial mixing conditions. We show that different initial mixing protocols lead to unconventional LLPS pathways, which can significantly accelerate the LLPS dynamics, compared to that in the traditional droplet coarsening. These findings highlight the critical role of initial mixing conditions in controlling the LLPS dynamics of complex coacervation, with implications in biotechnological processes, smart materials development, and the design of bioinspired adhesives.

Results
Morphology evolution for distinct initial mixing protocols
We investigate the LLPS dynamics starting from two distinct initial mixing scenarios between the polycations and polyanions. In the first scenario, the system consists of an initially well-mixed solution of polycations and polyanions, mimicking experimental conditions where LLPS starts with a dense mixture and is subsequently triggered by environmental cues22,25,26. We refer to this protocol as the well-mixed pathway (Fig. 1b). The second scenario involves an initial system composed of spatially segregated domains of polycations and polyanions, mimicking the “secrete and flux" coacervation process observed in certain marine organisms9–11. We refer to this setup as the flux pathway (Figs. 1c, 3d). For comparison, we include results from the thermodynamic pathway (Fig. 1a).

We first examine the structural evolution under different pathways, as shown in Fig. 1. In contrast to the thermodynamic pathway (Fig. 1a), where dispersed polycation-polyanion pairs coalesce into progressively larger droplets, the well-mixed and flux pathways exhibit drastically different morphology evolution. For systems starting with a well-mixed solution, the most characteristic feature is the formation of a network structure in the early stages of coacervation, which subsequently breaks down in later stages, eventually evolving into a single coacervate droplet (Fig. 1b). We note that the early-stage network formation and the corresponding dynamics have been recently reported by Yuan and Tanaka30. However, the complete pathway, including the breaking of the network and the later stage dynamics, has not been reported. In systems beginning with spatially segregated polycation and polyanion domains, coacervate droplets form rapidly in the contact zone due to the strong local electrostatic attraction, which is followed by the accumulation of dispersed chains into these interfacial droplets at a later stage, as illustrated in Fig. 1c. The morphological evolution in these two pathways is markedly different from the droplet coarsening mechanism in the thermodynamic pathway (Fig. 1a). Despite the differences in morphological evolution, all pathways ultimately converge to the same final state with a single coacervate droplet13. We note that the droplets in our simulations generally have floppy surfaces, indicating they have low surface tension, which is consistent with previous theoretical and experimental studies in coacervate systems37–39. These results also highlight the efficacy of our GCMe-DPD method in effectively simulating the entire LLPS process of coacervation, capturing the full dynamics from the initial state to the final equilibrium state.

The distinct morphological evolution observed in the well-mixed and flux pathways fundamentally alters the scaling behavior of the domain growth dynamics during LLPS, including both the early and later stages of the process, as detailed below.

Domain growth in the well-mixed pathway
In the early stage of the well-mixed pathway, where the polymer morphology takes on a 3D network structure, the domain growth exhibits a power law scaling of R~t0.46, as shown in Fig. 2a. This growth is significantly faster than that observed under the thermodynamic pathway, where the exponent is 1/3 due to the coarsening of dispersed droplets (Fig. 2d). The rapid network coarsening phenomenon in the early stages of the well-mixed pathway aligns closely with a recent finding by Yuan et. al, where they found polyelectrolyte network structures lead to domain growth scaling as R~t1/230. The connected network formed during this early stage arises from the interplay between local charge balance, low surface tension, and polymer connectivity30. Driven by the minimization of surface energy, this spatial network connectivity and hydrodynamics accelerate the coacervation process, resulting in an exponent close to 1/2, akin to the domain evolution of a viscoelastic polymer network31. While the 1/2 scaling appears universal in network-driven phase separation, its theoretical derivation remains elusive, as stated by Tanaka30,31. Future work is needed to understand this intriguing scaling.

At the later stage of the well-mixed pathway, the network breaks and the domain growth starts to deviate from the t1/2 scaling. From numerous independent runs, we generally observe two distinct types of later-stage dynamics: one exhibiting R~t1/3 and the other showing R~t1, as illustrated in Fig. 2b, c. The type of growth in the later-stage depends on the initial degree of mixing of polycations and polyanions. To quantify the initial degree of mixing, we divide the simulation box into Nc small cells of dimensions 10σ0 × 10σ0 × 10σ0, where 10σ0 is roughly the end-to-end distance of the polyelectrolyte chain in a single-species semi-dilute solution at the same concentration. We then characterize the degree of mixing by the magnitude of macromolecular charge imbalance, defined as

λ = 1/Nc Σ i=1 to Nc |np+(i) - np-(i)|, (1)

where np±(i) is the number of positive/negative monomers within cell i. On average, each bin contains 100 charged monomers; therefore, λ/2 characterizes the percentage (%) of the deviation from charge neutrality in each bin. Generally, when the initial mixture is prepared by randomly distributing the polyelectrolytes into the simulation box (as detailed in Supplementary Information), the system yields λ≈16.4 for the initial configurations, corresponding to 8.2% deviation from the charge neutrality per bin. Under this condition, the later stage dynamics obeys R~t0.34 (Fig. 2b), resembling the dynamics associated with droplet coarsening. This droplet-coarsening feature is more clearly observed from Supplementary Movie S1 in the Supplementary Information, where small droplets move and coalesce with each other in the later stage until the formation of a single large droplet. Our results show that the domain growth at a later stage is mainly driven by droplet-coarsening instead of the diffusion-limited aggregation (DLA) of solid-like precipitates, consistent with recent simulation results13 and recent experimental observation of biomolecular condensate formations22,23. The coarsening behavior is due to the collapse and fragmentation of the early-stage network caused by the imbalance of local macromolecular charge, resulting in multiple more charge-balanced droplets (Fig. 2e). As such the later-stage coacervation proceeds via droplet coarsening, similar to that in the thermodynamic pathways shown in Fig. 2d. A full domain evolution under this pathway is shown in Supplementary Fig. S4a and movie S1 at the Supplementary Information.

On the other hand, a fully locally charge-balanced initial state can be prepared by carefully positioning polycations and polyanions adjacent to each other, resulting in λ=0. Under this optimally mixed configuration, the later-stage dynamics exhibit a much faster scaling with R~t0.93 (Fig. 2c, f). This scaling is characteristic of the hydrodynamic pumping mechanism14, where the connectivity of the dense domains is coupled to a strong inertial effect, leading to a scaling transition to R~t1 at the later stage. Such a transition pathway is illustrated in Fig. 2f. In this case, the network formed during the early stage maintains local macromolecular charge balance and does not break into multiple droplets over time. Instead, the network coarsens in a continuous fashion into larger domain sizes, driven by rapid hydrodynamic pumping (Fig. 2f). In rare occasions when the initially mixture is prepared stochastically, we also observe similar hydrodynamic pumping at the later stages when the local macromolecular charges happen to be nearly balanced. A full domain evolution from early stage networking to later-stage hydrodynamic pumping is shown in Supplementary Fig. S4a at the Supplementary Information.

Domain growth in the flux pathway
For the flux pathway initiated with spatially separated domains, the system exhibits very fast LLPS dynamics with R~t0.65 in the early stage, as shown in Fig. 3a. During this initial phase, oppositely charged polyelectrolytes near the interfacial region quickly flow into each other, leading to the rapid formation of a few droplets in the contact region, with many nearby uncomplexed chains, as illustrated in Fig. 1c. The scaling exponent is reminiscent of the hydrodynamic-flux driven coarsening dynamics in deeply quenched binary liquid mixtures, which was predicted to follow R~t2/340. Although the direct connection between these two systems is not clear, both cases involve strong inertia effects. A detailed derivation of this scaling will necessarily couple electrostatics, hydrodynamics, and polymer conformation with the nonequilibrium phase separation, which is well beyond the scope of this work. Here, we provide a phenomenological understanding of how the electrochemical potential gradient in the flux pathway can drive this accelerating LLPS dynamics.

The electrochemical potential consists of a concentration part and an electrostatic part, both favors bringing the two oppositely charged polyelectrolytes into the contact region. The electrostatic potential Ψ(x,t) can be approximated by considering a 1-dimensional charge density profile in the x direction, and is given by the Poisson equation:

∂2Ψ(x,t) / ∂x2 = - ρch(x,t) / ϵ ϵ0 r, (2)

where ρch(x,t) is the net charge density distribution (including polycations, polyanions, and the corresponding counterions, see detailed numerical calculation in Supplementary Information). Ψ contributes to the electrochemical potential of the polyelectrolyte as μe p(x) = q(x)Ψ(x), where q(x) is the charge density of polyelectrolytes. At the early stage, due to the separation of polycation and polyanion regions, μe p is significantly higher at the two sides than in the interfacial region, as shown in Fig. 3c. Consequently, there is a strong electrochemical potential driving force for the polyions from the two sides to flow into the interfacial region, until Δμe p drops to below 10 kb T (Fig. 3c). Note that this potential (Donnan potential) is a consequence of macromolecular charge separation and therefore is undiminished by screening. Thus, although the nominal Debye screening length (0.34 nm, or 0.48 nm if only counterions are considered; see Supplementary Information) is much smaller than our simulation box, this electrochemical potential driving force will take effect. Indeed, the 2/3 law persists much longer when we increase the systems three times larger in the x direction (see details in section E, Supplementary Fig. S5 in Supplementary Information), indicating the timespan of this scaling regime is likely to be extended in larger systems. In real systems, the two separated domains are much larger than the box size in our simulation; therefore, we expect the t2/3 growth to persist until most of the polyions merge with the clusters in the interfacial region. We believe that this flux-driven LLPS is the key for the rapid production of adhesive coacervates in marine organisms such as mussels and Sandcastle worms, as schematically illustrated in Fig. 3d.

At the later stage of the flux pathway, coacervation proceeds by further adsorbing the uncomplexed polyions into the droplets, as well as the coarsening of the droplets towards complete coacervation, as shown in Fig. 1c. Theoretically, droplet coarsening should result in domain growth with R~t1/3. Our results show that the dynamics in the later stage of the flux pathway has scaling R~t0.39, as shown in Fig. 3b. The slightly faster scaling than the expected t1/3 is probably due to the residual electrochemical driving force for the polyions to migrate towards the interfacial region, as schematically shown in the inset of Fig. 3b, c. A full domain evolution in the flux pathway is given in Supplementary Information, Supplementary Fig. S4b.

The different scaling laws in the different pathways have significant consequences for the macroscopic phase separation time in coacervation. For the flux pathway, if we extrapolate using the early-stage scaling—since we expect this regime to persist for macroscopically large systems—the formation of a coacervate with a domain size of R=0.6 cm will only take about 10 s (based on taking τ0 = 7.8 ps and σ0 = 0.6 nm; see details in Supplementary Information). From the same estimate, in order to reach a state of similar domain size, it will take 47.5 years via the thermodynamic pathway! For the well-mixed pathway, the crossover from the early stage to the later stage is independent of system size14 and the total time it takes to reach a droplet of macroscopic size is dominated by the later-stage dynamics. With hydrodynamic pumping, macroscopic coacervation can form rapidly due to the linear scaling. If the later stage follows droplet coarsening, the rate of coacervation will be similar to the thermodynamic pathway. Note that such extrapolations are predicated on the assumption that there are no other macroscopic correlation lengths that take effect to change the scaling. The validity of this assumption remains to be tested.

Concentration dependence in early-stage dynamics
The early-stage network formation in the well-mixed pathway depends on polyelectrolyte concentration (ϕ)30. To examine this dependence, we quantify the polyelectrolyte concentration in terms of the overlap concentration ϕ*=0:64×6×N=(πR3 ee), where Ree is the polymer root-mean-square end-to-end distance in a salt-free polyelectrolyte solution41. The results discussed so far correspond to ϕ=2.7ϕ*. As shown in Fig. 4(a), LLPS dynamics are slower in more dilute systems (ϕ=1.0ϕ*, 2.0ϕ*). At ϕ=2.0ϕ*, rather than forming a network, oppositely charged polyelectrolytes associate into small droplets, which then coarsen with an exponent of 0.31, resembling spinodal decomposition. At ϕ=1.0ϕ*, local associations of polycations and polyanions lead to charge-asymmetric droplets coexisting with naked polyions (Fig. 4b), further slowing coarsening (Supplementary Fig. S6a)13. On the other hand, if the polymer concentration is significantly lower (ϕ≪ϕ*, but still larger than the coexisting supernatant concentration), the polyions are expected to form charge-balanced pairs first, and the coarsening dynamics should follow R~t1/3. However, direct simulation of LLPS of dilute systems starting with dispersed chains is not computationally feasible within our GCMe-DPD framework, as it would require a much larger system to obtain a statistically reliable domain growth law. In Supplementary Information, we increase the chain length to N=50 with varying concentrations, the results are qualitatively the same. These findings indicate that network-accelerated LLPS requires relatively high polymer concentrations (ϕ ≈ 3ϕ*, Fig. S6c), consistent with the recent finding30. However, in real systems, the molecular weights of the polymers are often quite high, naturally satisfying this requirement since ϕ*~N−242. Thus, LLPS dynamics via the formation of a transient network are expected to be more common in experiments.

Methods
Coarse-grained model and simulation methodology
Our simulations employ the Gaussian Core Model with Smeared Electrostatic (GCMe) Interactions51,52 combined with dissipative-particle-dynamics (DPD) thermostat to capture hydrodynamic interactions53–55. This method, which we refer to as the Gaussian Core Model with Smeared Electrostatics and Dissipative Particle Dynamics (GCMe-DPD), enables efficient modeling of both the equilibrium and nonequilibrium dynamics of polyelectrolyte solution systems. The soft-core Gaussian potentials used in GCMe can be analytically derived from Gaussian-smeared mass and charges51. Compared to other coarse-grained models such as the Lennard-Jones based bead-spring models44,46,47 or other chemically specific models such as MARTINI56,57, our GCMe-DPD has the advantages in capturing the long-range hydrodynamic and electrostatic interactions, both of which are critical for the non-equilibrium LLPS dynamics. In addition, the high computational efficiency of GCMe-DPD allows us to explore the full LLPS kinetics (from coarsening to complete phase separation) of large coacervate systems, as evidenced in our results in both the main text and Supplementary Information.

Simulation details
In our simulations, the polyelectrolyte monomers, counterions, and solvent molecules (water) are represented as coarse-grained beads, with each bead having a Gaussian distribution in mass and charge density around their center. We set the excluded volume interactions to be the same among all species to represent good solvent quality for the polymers. We use nondimensionalized units in which the mass, energy, and length are scaled by the particle mass m0, the thermal energy kbT at room temperature, and the Gaussian-core particle size σ0, all taken to be unity. Each simulation contains an equal number of fully charged polycations and polyanions of the same chain length, accompanied by their corresponding monovalent mobile counterions to maintain overall charge neutrality. Unless otherwise stated, our simulation box has dimensions of 60σ0 × 60σ0 × 60σ0 and contains 4% polyions with chain length N=40, 4% counterions, and 92% water, representing a semi-dilute polyelectrolyte solution. The Bjerrum length lB is set to match the dielectric environment of water at room temperature. The electrostatic strength is lB/b≈1.37 (where b is the charge spacing along the neighbored polymer bead), which falls into the intermediately charged regime for typical polyelectrolytes. Additional simulation details, including choice of interaction parameters, the simulation setup, and analysis methods, are provided in the Supplementary Information and in our original methodology papers51,52.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "A sufficiently detailed understanding with new quantification of the dynamics of temporal pathways of LLPS of polymers, as influenced by initial state of the ensemble of oppositely charged polyelectrolytes in salt-free aqueous solution is reported. The results are original and of significance to a broad audience across various fields. The DPD method of simulation used in correct (including the hydrodynamic pumping) and the methods are appropriate. The system (simulation box) dimensions, polymer chain lengths, DPD parameters for chain structure and charges, the simulation times, are all appropriate and consistent for the right representation of a reasonable physical system. All necessary details are provided appropriately in the manuscript and supplemental information section. The paper is well written. The theoretical model (Poisson's equation) for the 1-D flux-pathway is simple yet provides an explanation within the scope of this study, some theoretical basis for understanding the origins of the dynamical behavior for this particular case. The work carried out in this study supports the conclusions and claims of the authors, and there is no additional evidence needed. This work meets the expected standards in the field, and also in suitable for the journal Nature Communications. The work is relevant to the field of synthetic polyelectrolytes, polymer science and engineering, natural sciences, and to a little extent (though restricted extent and in a qualitative sense) to processes that occur in biology. The noteworthy results of this study are the new dynamical scalings of LLPS of oppositely charged polyelectrolytes towards the formation of polyelectrolyte complexes (PEC) reported for the flux-pathway and for the well-mixed pathway (which is driven by network formation initially itself), and the distinctly different dynamical scalings observed both these pathways for initial and later phases of PEC formation. These results are remarkable and useful to further our understanding of the dynamics of formation of such complex systems that occur in nature. I recommend publication in the present form.", "reviewer_2": "Overall, this is an interesting and much-needed work investigating the dynamics of phase separation in coacervates and why they may differ from thermodynamic phase separation and dynamics of phase separation in neutral polymer systems. In particular, the authors’ main question and approach are good ones; their comparison of the impact of initial condition is relevant to recent experimental literature that observes impact on mixing parameters on coacervate morphology and growth, and also requires the simulation of counter ions (thus use of MD simulation is a good choice).\n\nHowever, there are a number of Major Revisions to address, before the reviewer recommends the article for acceptance to publication. See below.\n\nMajor Comments:\n1. The authors only include a single paragraph for Discussion, and no Conclusion. In addition, the first subsection of the Results is not a good place to describe the Methods. This should be improved.\n2. Methods: For good reproducibility of your work and easy readability, make sure you include all the details of the methods.\nSome items I noticed were:\ni) GCMe-DPD abbreviation is not defined in the main text.\nii) \"We set the excluded volume interactions to be the same among all species to represent good solvent quality for the polymers\". Does this mean the excluded volume interactions for counterions, water molecules, and polymer ions are treated equally? Please explain this choice.\niii) There are many methods to coarse-grain water. Please state which you model chose and why, and similarly for the other components. These choices can affect the results of the work.\niv) It's critical to maintain charge neutrality, yet the counterions are mobile. How can you maintain local electroneutrality, is there a constraint that each ion must always be paired with another?\nv) In brief, the authors should make sure to address what assumptions are made in the simulations. In addition, the authors should briefly address what the method does better (or worse) than other possible approaches. Give some sense of to what degree the results are simplified compared to real life.\n3. Slopes: The authors labeled in most cases 2 different (but very similar) values of the slope as \"bounds\" on each simulation run, but did not clarify:\ni) What range of data do these slopes represent and how were they obtained? They are reported currently with a possibly-misleading amount of certainty, i.e., without error, even though obtained from the data. Uncertainty estimates and/or variational analysis, or at minimum the standard protocol which was used across all datasets to obtain these \"upper\" and \"lower\" bounds, should be included in the Supplemental Information.\nii) In many cases, the values of the bounds are very similar. Are these \"upper\" and \"lower\" boundaries actually necessary? The reviewer found them to be distracting from the beta-values, and in addition, they were barely mentioned/discussed in the work, despite being in nearly every plot. Do these bounds serve some conceptual purpose that is not currently apparent to the reviewer, beyond giving a sense of the uncertainty in beta-value?\niii) Relatedly, beta-values are reported with very tiny errors, of order 0.1 +/- 0.001. This seems underestimated given the error bar size. Have these errors been calculated when accounting for uncertainty in the data? In addition, consider reporting a 95% CI of beta instead and getting rid of the bounds, unless they have an additional purpose.\niv) If there is a utility to the bounds beyond giving a sense of the uncertainty in beta, it would strengthen the work to discuss their meaning. Please discuss how these boundary values should be interpreted, including a comparison to literature.\n4. To report the results, you should be sure that there are no confinement effects, since box size can significantly affect the coarsening rate. It is mentioned that the system size influences the duration of the early stage in the flux pathway, and that \"Unless otherwise stated, the simulation box has dimensions [60 × 60 × 60]\". Did you try a larger box size for any other pathways, or how are you sure that the box size is not affecting other results?\n5. The author’s main point in the Introduction is that various “initial mixing conditions can significantly influence and accelerate the phase separation dynamics”.\ni) Though this is the case, the authors take a strange interpretation of the literature, comparing “slow” coacervate in “hours or days” with “faster” times. In general, this is not consistent with experimental literature. Ref 25,26 do not mention coacervate “forming over hours or even days”. (bottom of page 2). Where is the author finding this claim? In fact Ref 25 mentions “that primary complexes form very fast within the dead time” of stopped-flow mixing.\nii) Some works which show a pathway of formation of primary aggregates during mixing, which is generally attributed to diffusion-limited attraction of oppositely-charged polymers during mixing (and complete phase separation during mixing) are: Trinh et al. Macromolecular Chemistry and Physics 1997, Liu et al. Advances in Colloid and Interface Science 2017, and Liu et al. Soft Matter 2016.\niii) One work which supports your main point in the introduction and the motivation of the work is Edwards et al., Soft Matter 2023 (and in supplemental) for the flux pathway\niv) For some other systems, scientists must wait minutes after mixing to see phase separation, such as poly(vinylpyridine) and poly(methacrylic acid) in Trinh et al.; Salvatore et al. Biomacromolecules 2011; Mohammadi et al. Journal of Colloid and Interface Science 2020; and Sanchez et al. Journal of Colloid and Interface Science 2006.\nIn general, the Introduction and Discussion should account for an understanding of the results in these and/or similar works at the authors' discretion.\n6. Experimental works describing precipitate versus coacervate growth in the same types of chemistries have observed diffusion Mejri and Aschi, Colloids and Surfaces B: Biointerfaces 2022 show that precipitate growth resembles diffusion-limited aggregation, instead of liquid-like growth signatures, but these types would be indistinguishable using a structure-factor based measure of structure used in this work. Still, given the need for more Discussion about growth rates, could the authors comment more on which types of growth mechanisms could correspond to the observed exponents?\n7. The claim made at the end of the paper \"Additionally, ultrafast coacervation can occur in the flux pathway, driven by electrochemical potential gradients from the initially separated polycation and polyanion domains, mimicking the rapid adhesive processes observed in marine organisms\" is quite strong, when compared to the what was found in the Results, namely that \"The slightly faster scaling [0.39 exponent] than the expected t^1/3 is probably due to the residual electrochemical driving force for the polyions to migrate towards the interfacial region, as schematically shown ...\". It would make the paper stronger to discuss this result more. Why is it \"ultrafast\" in the context of literature for theory, simulation, and/or experiment? What is the meaning of the t^(1/3) exponent? You should discuss LSW theory/Ostwald Ripening, and diffusion-limited growth, and explain the 0.39 exponent within this context; and whether/what faster exponents have been previously observed in experiment, simulation, or estimated in theory. Also, it is important to clearly differentiate in the Conclusions what is suggested by the results (e.g., electrochemical potential drives ultrafast phase separation), versus what you have directly observed in them (an increase in domain growth exponent to 0.39 versus the 0.33).\n7. Sentence pg. 12: “On the other hand, if the polymer concentration is significantly lower ... but still larger than the coexisting supernatant concentration), the polyions are expected to form charge-balanced pairs first, and the coarsening dynamics should follow R ~ t^(1/3)”. But, this claim was not investigated. It would be nice to see it using a simulation, or is there a reason you did not include this data?\ni) At the overlap concentration, can you quantify how much the droplets are charge-asymmetric, and how much this “further slows” coarsening?\nMinor comments:\n1. First sentence of the Results is awkward, rephrase.\n2. Lambda values are only labelled on some of the plots.\n3. The letter labels (a), (b), (c), (d) in the Figures are often in a strange location, e.. in the middle of each subplot.\n4. Please double-check the work for sentence-level edits. For example, a grammatical error in the very 1st sentence of abstract can obfuscate your meaning. Check throughout.\ni) In addition, 1st sentence pg. 3: How does a “thermodynamic pathway ... model the coarsening dynamics”? Consider rewriting to clarify what point you are making here.\nii) The proper phrasing in Acknowledgements is “We thank X for their generous allocation of Y”.\n5. In Figure S3 and related discussion, make sure to acknowledge Figure S4.", "reviewer_3": "Manuscript Summary and Key Findings\nIn this manuscript, the authors present a study on the effects of mixing protocols on the dynamics of liquid-liquid phase separation of polyelectrolyte complex coacervates. The authors use molecular dynamics simulations to evaluate the domain growth of the coacervate droplets for two distinct initial mixing conditions. The first is a well-mixed solution of oppositely charged polyelectrolytes. The second scenario involves a system where the polycations and polyanions are separated into domains that are put in contact with each other, the authors called this scenario the flux pathway.\n\nThe authors found that the two initial mixing scenarios produced different domain growth behaviors. The first scenario, the well mixed pathway, showed the formation of a network structure in the early stages that would then break down into droplets which eventually merge into a single droplet. In the flux pathway, coacervate droplets will form at the interface between the two domains, which will then merge and draw in free chains to grow the domains. The well mixed pathway showed early-stage growth with scaling of t1/2, and the flux pathway showed an early-stage domain growth with scaling t2/3.\n\nThe authors also evaluated the polymer concentration effects on the domain growth mechanism and dynamics and found that the network structure that leads to enhanced early-stage growth requires a sufficiently high concentration to observe the enhanced growth dynamics. The authors showed that the initial mixing conditions can be used to accelerate coacervation dynamics, and the implications of system size are discussed for the various initial mixing conditions.\n\nReview\nThis manuscript provides valuable insight into the dynamics of polyelectrolyte coacervation that will be of considerable interest to the coacervate and polyelectrolyte complex communities. The findings in this paper are, to my knowledge, novel and important. I believe the methods used by the authors are valid for the study, and the findings presented within the manuscript are well reasoned and supported by the data. The manuscript is clear and concise, and provides the context necessary to understand the significance and motivations of the study. The references cited in the manuscript are appropriate and many are recent. The information provided in the supplementary materials is appropriately placed as supplemental and provides more details regarding the simulation methods and data analysis. The methods have been sufficiently described to allow for reproduction of the study by other researchers.\n\nThe modeling and simulation methods used in this study was developed and previously published by the corresponding authors, and addresses some key issues known within community. Molecular dynamics simulations of complex coacervates can be complicated by the length and timescale requirements to observe phase separation, as well as the nature of the hard cores of the Lennard Jones potential and the treatment of bulky ions as point charges within the center of the ion bead. The Gaussian Core Model with Smeared Electrostatic Interactions methodology presented by the authors addresses these three potential issues while maintaining the ability to evaluate dynamics and consider hydrodynamic interactions. This methodology could be used by other researchers studying coacervation and other solutions of charged macromolecules, particularly for out of equilibrium and dynamics studies.\n\nIt is my belief that the data from which the scaling exponents are drawn from are high quality and robust. The authors in this manuscript provide well-reasoned explanations with mathematical descriptions for some of the findings in the study. For example, they found two different late-stage domain growth dynamics for the well mixed pathway that they attributed to the initial degree of mixing. The authors provided a method for evaluating the degree of mixing and found that indeed the value for this degree of mixing correlated with the late-stage domain growth dynamics. Another example of the providing rigorous explanations for their findings can be found in their explanation for the driving force that leads to the t2/3 domain growth scaling observed in the flux pathway. The authors mathematically described the electrochemical potential and showed how the potential difference between the interface and the two initial domains decreases over time to lead to a reduction in the domain growth dynamics between the early stage and later stage of the flux pathway coacervation.\n\nThe authors’ discussions and findings regarding the concentration dependence on early-stage domain growth dynamics and the system size implications are very insightful. For example, the authors extrapolated the scaling for macroscopically large systems, and found that the flux pathway can undergo coacervation in a matter of seconds, but found that the baseline thermodynamic pathway would take 47.5 years to form a coacervate domain of the same size.\n\nIt is not entirely clear to me the level of significance the broader audience of Nature Communications would find in this manuscript, though I believe this topic is of considerable interest to the soft materials scientific community. Coacervates represent a potential environmentally friendly and sustainable solution for many materials applications, where these characteristics are becoming increasingly sought after. I believe this manuscript is of high significance and high quality and should be published in Nature Communications. I have a few questions for the authors to possibly consider discussing within the manuscript.\n\nQuestions for the authors\n1. One key finding in this study is the enhanced domain growth as a result of network formation. This network formation is concentration dependent, and I am wondering about the implications of the inhomogeneous polarizability throughout the system that will be more prevalent at higher concentrations. How does the model account for these local changes in polarizability between water dominant regions and polymer dominant regions?\n2. What are the implications of using the same bead size for the salt counterions as the other species, on the dynamics of coacervation? How accurate are the salt ion mobilities captured in this model and could deviations to the salt ion mobilities lead to significant changes in the mechanisms and dynamics of coacervate domain growth?\n3. It seems to me that the interface surface tension between the coacervate domain and the solvent and free chain domain is important for the dynamics of droplet coarsening. Is this the case and how accurately is this interface surface tension being modeled with this coarse-grained model? Would you expect the structuring of water around the coacervate droplets to have an impact on the droplet coarsening dynamics?\n4. I am wondering about the domain size calculations from the structure factor and the error in the domain size as a result of noise in the structure factor profiles. Do you have a sense for how sensitive the domain size calculation is to the choice of q-vector binning? Presumably, larger systems would allow for obtaining structure factor profiles with less noise, and so I am wondering if the systems studied are sufficiently large to get accurate domain sizes from the structure factor."}, "author_response": "Responses to the Comments of Reviewer #1\n“A sufficiently detailed understanding with new quantification of the dynamics of\ntemporal pathways of LLPS of polymers, as influenced by initial state of the ensemble\nof oppositely charged polyelectrolytes in salt-free aqueous solution is reported. The\nresults are original and of significance to a broad audience across various fields.\nThe DPD method of simulation used is correct (including the hydrodynamic pumping)\nand the methods are appropriate. The system (simulation box) dimensions, polymer\nchain lengths, DPD parameters for chain structure and charges, the simulation times,\nare all appropriate and consistent for the right representation of a reasonable physical\nsystem.\nAll necessary details are provided appropriately in the manuscript and supplemental\ninformation section. The paper is well written. The theoretical model (Poisson's\nequation) for the 1-D flux-pathway is simple yet provides an explanation within the\nscope of this study, some theoretical basis for understanding the origins of the\ndynamical behavior for this particular case.\nThe work carried out in this study supports the conclusions and claims of the authors,\nand there is no additional evidence needed. This work meets the expected standards in\nthe field, and also is suitable for the journal Nature Communications. The work is\nrelevant to the field of synthetic polyelectrolytes, polymer science and engineering,\nnatural sciences, and to a little extent (though restricted extent and in a qualitative sense)\nto processes that occur in biology.\nThe noteworthy results of this study are the new dynamical scalings of LLPS of\noppositely charged polyelectrolytes towards the formation of polyelectrolyte\ncomplexes (PEC) reported for the flux-pathway and for the well-mixed pathway (which\nis driven by network formation initially itself), and the distinctly different dynamical\nscalings observed both these pathways for initial and later phases of PEC formation.\nThese results are remarkable and useful to further our understanding of the dynamics\nof formation of such complex systems that occur in nature.\nI recommend publication in the present form.”\nReply:\nWe sincerely thank the reviewer for the high evaluation of our manuscript and the\nsupport for its publication in Nature Communications in its current form.\nResponses to the Comments of Reviewer #2\n“Overall, this is an interesting and much-needed work investigating the dynamics of\nphase separation in coacervates and why they may differ from thermodynamic phase\nseparation and dynamics of phase separation in neutral polymer systems. In particular,\nthe authors’ main question and approach are good ones; their comparison of the impact\nof initial condition is relevant to recent experimental literature that observes impact on\nmixing parameters on coacervate morphology and growth, and also requires the\nsimulation of counterions (thus use of MD simulation is a good choice).”\nWe thank the reviewer for the recognition of the significance and timeliness of our work.\nBelow, we provide a point-by-point response to the reviewer’s comments.\nMajor Comments:\n1. “The authors only include a single paragraph for Discussion, and no Conclusion. In\naddition, the first subsection of the Results is not a good place to describe the Methods.\nThis should be improved.”\nReply:\nWe thank the reviewer for the suggestion on adding more discussions and reorganizing\nthe methodology section.\nIn the revision, following the reviewer’s suggestion, we have moved the methodology\ndescriptions to the end of the manuscript in the “Method” section. In addition, we have\nsignificantly expanded the discussion section (see Page 13&14 in the main text) to\naddress some of reviewer’s later comments and some of the reviewer #3’s comments.\n2. “Methods: For good reproducibility of your work and easy readability, make sure\nyou include all the details of the methods. Some items I noticed were:\ni) GCMe-DPD abbreviation is not defined in the main text.\nReply:\nWe thank the reviewer for carefully reviewing our manuscript.\nIn the updated manuscript, we add the full name of method before the abbreviation (see\nPage 15).\nii) \"We set the excluded volume interactions to be the same among all species to\nrepresent good solvent quality for the polymers\". Does this mean the excluded volume\ninteractions for counterions, water molecules, and polymer ions are treated equally?\nPlease explain this choice.\nReply:\nWe appreciate the reviewer’s suggestions on clarifying our simulation parameters. In\nour simulations, the excluded-volume interactions among all species (counterions,\nwater molecules, and polymer ions) are set to be the same, which represents a good\nsolvent condition for the polymers on the Flory-Huggins level. As the parameters in\ncoacervate systems are very large, this simplification allows us to focus on the effect of\nelectrostatics (and hydrodynamics) on the LLPS dynamics. We note that such\nsimplifications are widely used in coarse-grained simulations of generic coacervation\nsystems in the community to tackle the electrostatic effects ([J. Yuan & H. Tanaka, Nat.\nCommun., (2025)]; [Y. N. Fang, et al., Proc. Natl. Acad. Sci., (2023)]; [B. Yu, et al.,\nMacromolecules, (2023)]; [B. Yu, et al., Macromolecules, (2022)]; [D. Sharon, et al.,\nMacromolecules, (2022)]; [H. Liang, et al., Macromolecules, (2022)]).\nIn the revision, we add the following clarifications on Page S2 in SI: “In this work, we\nfollow the common practice in coarse-grained simulations ([J. Yuan & H. Tanaka, Nat.\nCommun., (2025)]; [Y. N. Fang, et al., Proc. Natl. Acad. Sci., (2023)]; [B. Yu, et al.,\nMacromolecules, (2023)]; [B. Yu, et al., Macromolecules, (2022)]; [D. Sharon, et al.,\nMacromolecules, (2022)]; [H. Liang, et al., Macromolecules, (2022)]) to set the\nexcluded volume interactions the same among all species (𝐴 = 25.2𝑘 𝑇𝜎%&),\nwhich represents good solvent conditions for polymers on the Flory-Huggins level. This\nsimplified setup allows us to focus on the essential LLPS dynamics driven by long-\nrange electrostatic and hydrodynamic interactions.”\niii) There are many methods to coarse-grain water. Please state which you model chose\nand why, and similarly for the other components. These choices can affect the results\nof the work.\nReply:\nWe agree with the reviewer that there are many coarse-grained methods. In this work,\nwe used our developed GCMe-DPD method based on the following reasons:\n1) Our GCMe-DPD method effectively couples long-range hydrodynamic interactions\n(HI) and electrostatic interactions, which are the two major factors for the non-\nequilibrium LLPS dynamics focused on this work. To our knowledge, while other\ncoarse-grained methods commonly used in modeling coacervate systems have their\nmerits, such as Leonard-Jones (LJ) based bead-spring methods ([B. Yu, et al.,\nMacromolecules, (2023)]; [D. Sharon, et al., Macromolecules, (2022)]; [H. Liang, et\nal., Macromolecules, (2022)]), and Martini force field ([M. Tsanai, et al., Chem. Sci.,\n(2021)]; [L. Wang, et al., Nat. Commun., (2025)]), these methods (with or without\nexplicit solvents) fail to conserve local momentum and thus do not rigorously capture\nhydrodynamic interactions. Given the importance of HI in non-equilibrium LLPS\ndynamics, our GCMe-DPD offers clear advantages.\n2) The choice of using Gaussian soft-core potentials in the coarse-grained excluded\nvolume and electrostatic potentials is due to their theoretical tractability, i.e., they can\nbe analytically derived respectively from Gaussian-smeared mass and charges and are\nmore amenable to analytical manipulations, as documented in our previous publication\n[B. Ye, S. Chen & Z.-G Wang, JCTC, (2024)].\n3) The use of soft-core potentials and efficient implementation in GPUs make our\nGCMe-DPD method feasible to model large coacervate systems that are necessary for\nexploring long-time LLPS dynamics. Note that our simulation systems are the largest\n(in terms of the number of polymer beads) coacervate systems in the literature on\nstudying the LLPS dynamics. For example, in this work, GCMe-DPD can capture the\nfull LLPS pathway of a system containing 64800 polymer beads, with explicit\nhydrodynamics and electrostatics.\n4) While some effects, such as the explicit hydration of the ions and the related specific\nion effects, and the explicit polarization of the solvent molecules, are missing in this\nmodel, we believe the model captures the leading order effects due to both electrostatic\ncorrelations and hydrodynamic interactions, and the level of treatment is consistent with\nthe majority of existing models on coacervation phase equilibrium and dynamics, as\nmentioned above.\nIn revision, we add the following rationales in Page 15 in the main text:\n“The soft-core Gaussian potentials used in GCMe can be analytically derived from\nGaussian-smeared mass and charges [B. Ye, S. Chen & Z.-G Wang, JCTC (2024)].\nCompared to other coarse-grained models such as the Lennard-Jones based bead-spring\nmodels ([B. Yu, et al., Macromolecules, (2023)]; [D. Sharon, et al., Macromolecules,\n(2022)]; [H. Liang, et al., Macromolecules, (2022)]) or other chemically specific\nmodels such as MARTINI ([M. Tsanai, et al., Chem. Sci., (2021)]; [L. Wang, et al., Nat.\nCommun., (2025)]), our GCMe-DPD has the advantage in capturing the long-range\nhydrodynamic and electrostatic interactions, both of which are critical for the non-\nequilibrium LLPS dynamics. In addition, the high computational efficiency of GCMe-\nDPD allows us to explore the full LLPS kinetics (from coarsening to complete phase\nseparation) of large coacervate systems, as evidenced in our results in both the main\ntext and Supplemental Material.”\niv) It's critical to maintain charge neutrality, yet the counterions are mobile. How can\nyou maintain local electroneutrality, is there a constraint that each ion must always be\npaired with another?\nReply:\nIn our simulations, we only maintain the global charge neutrality by including an equal\namount of small counterions for the polymers to the whole simulation box. Like all the\ncommon MD simulations, we do not impose any constraints on the small ions and allow\nthem to follow their trajectories in the NVT ensemble. On small length scales (within\nthe Debye screening length), deviation of charge neutrality is possible due to\nfluctuations. However, beyond the Debye screening length, charge neutrality will\nnaturally be met in an ensemble-average sense.\nIn revision, we add the following clarifications on page 15 in the main text: “Each\nsimulation contains an equal number of fully charged polycations and polyanions of the\nsame chain length, accompanied by their corresponding monovalent mobile\ncounterions to maintain overall charge neutrality.”\nv) In brief, the authors should make sure to address what assumptions are made in the\nsimulations. In addition, the authors should briefly address what the method does better\n(or worse) than other possible approaches. Give some sense of to what degree the results\nare simplified compared to real life.”\nReply:\nWe thank the reviewer for the valuable suggestion of clarifying the assumptions,\nadvantages, and limitations of our simulation methodology. We have already listed the\nadvantages of our method in the reply to their previous comment. Indeed, as we focus\non the essential physics of LLPS dynamics on length and time scale well beyond the\nresidue level, our simulation methodology has the following major simplifications and\nthe corresponding limitations: 1) our model treats the excluded volume interaction\namong all coarse-grained beads (including polymers, counterions, and water)\nidentically, representing a good solvent condition for the polymers; As such, our model\ndoes not include specific chemical incompatibility arising from short-range interactions\namong different species (which are not difficult to include but introduces additional\neffects beyond those addressed in our manuscript). 2) We use a uniform dielectric\nconstant to model electrostatics; therefore, the possible local dielectric inhomogeneity\nis not accounted for. In principle, all these effects can be captured by all-atom\nsimulations of specific systems, but they are computationally prohibitive to explore the\nLLPS dynamics in such large coacervate systems.\nWe note that these simplifications are common practices in coarse-grained simulation\nin coacervate systems ([J. Yuan & H. Tanaka, Nat. Commun., (2025)]; [Y. N. Fang et\nal., Proc. Natl. Acad. Sci., (2023)]; [B. Yu, et al., Macromolecules, (2023)]; [B. Yu, et\nal., Macromolecules, (2022)]; [D. Sharon, et al., Macromolecules, (2022)]; [H. Liang,\net al., Macromolecules, (2022)]). As this work concerns the scaling law of LLPS\ndynamics on length and time scales well beyond the residue level, we thus believe that\nincorporating atomistic details will not qualitatively change the conclusions of this\nwork.\nIn response to the reviewer’s suggestion, we add these discussions in the “Discussion”\nsection (See page 13&14):\n“As the focus in this study is using coarse-grained MD simulation to elucidate the\nscaling behavior in LLPS dynamics of polyelectrolyte complex coacervations on length\nand time scale well beyond the residue level, our model follows the common practices\nof coarse-grained simulations ([J. Yuan & H. Tanaka, Nat. Commun., (2025)]; [Y. N.\nFang et al., Proc. Natl. Acad. Sci., (2023)]; [B. Yu, et al., Macromolecules, (2023)]; [B.\nYu, et al., Macromolecules, (2022)]; [D. Sharon, et al., Macromolecules, (2022)]; [H.\nLiang, et al., Macromolecules, (2022)]) to incorporate some necessary simplifications.\nFor example, our simulation systems treat the excluded volume interaction among all\ncoarse-grained beads (including polyions, counterions, and water) identically,\nrepresenting a good solvent condition for the polymers at the Flory-Huggins level. As\nsuch, our model does not include the solvent quality effects or specific ion effects\narising from specific short-range interactions among different species…”\n“… Furthermore, we use a homogeneous dielectric constant to model electrostatics,\nwhich does not account for possible local dielectric inhomogeneity arising from\npolymer aggregations. Accounting for this local dielectric effect requires better\ntheoretical understanding and new methodology development in the construction of\ncoarse-grained models. In principle, all these effects can be captured by all-atom\nsimulations of specific systems. However, such calculations are at present\ncomputationally prohibitive for exploring the long-time behavior of large coacervate\nsystems and are well beyond the scope of this work. As this work concerns the scaling\nbehavior in LLPS dynamics of generic coacervate systems on length and time scales\nwell beyond the residue level, we believe that incorporating atomistic details will not\nqualitatively change our conclusions.”\n3. “Slopes: The authors labeled in most cases 2 different (but very similar) values of\nthe slope as \"bounds\" on each simulation run, but did not clarify:\ni) What range of data do these slopes represent and how were they obtained? They are\nreported currently with a possibly-misleading amount of certainty, i.e., without error,\neven though obtained from the data. Uncertainty estimates and/or variational analysis,\nor at minimum the standard protocol which was used across all datasets to obtain these\n\"upper\" and \"lower\" bounds, should be included in the Supplemental Information.\nii) In many cases, the values of the bounds are very similar. Are these \"upper\" and\n\"lower\" boundaries actually necessary? The reviewer found them to be distracting from\nthe beta-values, and in addition, they were barely mentioned/discussed in the work,\ndespite being in nearly every plot. Do these bounds serve some conceptual purpose that\nis not currently apparent to the reviewer, beyond giving a sense of the uncertainty in\nbeta-value?\niv) If there is a utility to the bounds beyond giving a sense of the uncertainty in beta,\nit would strengthen the work to discuss their meaning. Please discuss how these boundary values should be interpreted, including a comparison to literature.”\nReply:\nWe thank the reviewers’ comments on the use of bounds of the slopes, and we agree\nwith the reviewers that such bounds are unnecessary and should be removed.\nIn revision, we have removed the “upper” and “lower” boundaries from all figures.\nInstead, we now present the uncertainty in the domain growth data using shaded areas\nto represent the standard deviations across the ten independent simulations in each plot\n(see Fig. 2(a)-(d), 3(a), 3(b), 4(a) in the updated main text and Fig. S1, S7 in SI).\nii) “Relatedly, beta-values are reported with very tiny errors, of order 0.1 +/- 0.001.\nThis seems underestimated given the error bar size. Have these errors been calculated\nwhen accounting for uncertainty in the data? In addition, consider reporting a 95% CI\nof beta instead and getting rid of the bounds, unless they have an additional purpose.”\nReply:\nThe original error analysis of beta values comes from fitting the domain growth data\npoints after taking the average of ten individual simulations. Indeed, such fitting\nunderestimates the errors from independent runs.\nIn response to the reviewer’s suggestion, in addition to plotting the standard deviations\nof independent runs, we now report the 95% confidence interval (CI) of beta from\nfitting all the individual runs (see Fig. 2(a)-(d), 3(a), 3(b), 4(a) in the updated main text\nand Fig. S1, S7 in SI).\n4. “To report the results, you should be sure that there are no confinement effects, since\nbox size can significantly affect the coarsening rate. It is mentioned that the system size\ninfluences the duration of the early stage in the flux pathway, and that \"Unless otherwise\nstated, the simulation box has dimensions [60 × 60 × 60]\". Did you try a larger box size\nfor any other pathways, or how are you sure that the box size is not affecting other\nresults?”\nReply:\nWe thank the reviewer for the suggestion on further discussing the finite-size effects in\nother pathways. In the thermodynamics pathway, the scaling law of t1/3 is independent\nof the system size once the box size exceeds the correlation length, which is widely\nestablished theoretically [Siggia, Phys. Rev. A, (1979)] and has been reproduced by our\nprevious work [S. Chen & Z.-G. Wang, PRL, (2023)]. Theoretically, the well-mixed\npathway should only have one correlation length and the scaling behavior should not\ndepend on system size (since the correlation length is determined by surface tension\nand is not related to system size).\nIn response to the reviewer’s suggestion, besides performing one large simulation of\nthe flux pathway (by setting 𝑙 = 180𝜎 ), we performed one additional simulation by\nincreasing the system size to 80 × 80 × 80𝜎& (see Fig. R1(a)) in the well-mixed pathway.\nFigure R1(a) shows that the characteristic scaling behavior for the well-mixed pathway\nremains relatively unchanged compared to original systems, including the occurrence\nof 𝑡(/* scaling at early stage, and the time to transition to 𝑡(/& (Fig. R1(a)). However,\nthe later stage in the well-mixed pathway has intrinsically large fluctuations (see Fig.\n2(b) in the main text), due to the small number of droplets at later time and the stepwise\ndroplet coarsening dynamics. Good statistics require many independent runs of large\nsimulations, which is too expensive to perform (Simulating a single system of this size\nrequires 1 month of running time on the high-performance cluster (HPC) with two\nRTX4090 GPUs). Note that our original system sizes are already the largest among\nexisting literature on LLPS dynamics of coacervate systems. Overall, our results\nindicate that the scaling behavior in the well-mixed pathway is independent of system\nsize (Fig. R1(a)), as expected.\nFor the flux pathway, we performed additional simulations by increasing the system\nsize three times along the flux direction. We found the 2/3 scaling persists much longer\nwith the increased box size (Fig. R1(b)), indicating that the scaling can be extrapolated\nfurther.\nIn revision, we add the results of large systems on well-mixed pathway, together with\nthe original results of the flux pathway, to section E of SI (see page S7&S8).\nFigure R1: Domain evolution with different simulation box sizes. (a) Well-mixed\npathway; (b) Flux pathway. The dashed lines mark the transition from the early stage\nto the later stage.\n5. “The author’s main point in the Introduction is that various “initial mixing conditions\ncan significantly influence and accelerate the phase separation dynamics”.\ni) Though this is the case, the authors take a strange interpretation of the literature,\ncomparing “slow” coacervate in “hours or days” with “faster” times. In general, this is\nnot consistent with experimental literature. Ref 25,26 do not mention coacervate\n“forming over hours or even days”. (bottom of page 2). Where is the author finding this\nclaim? In fact Ref 25 mentions “that primary complexes form very fast within the dead\ntime” of stopped-flow mixing.\nReply:\nWe thank the reviewer for carefully reviewing our manuscript and pointing out that our\ninitial statement is not consistent with experimental literature. We agree that our initial\nstatement regarding “slow” coacervation occurring over “hours or days” is incorrect.\nIn the revised manuscript, we have removed this claim (See page 2 in the main text).\nii) Some works which show a pathway of formation of primary aggregates during\nmixing, which is generally attributed to diffusion-limited attraction of oppositely-\ncharged polymers during mixing (and complete phase separation during mixing) are:\nTrinh et al. Macromolecular Chemistry and Physics 1997, Liu et al. Advances in\nColloid and Interface Science 2017, and Liu et al. Soft Matter 2016.\niii) One work which supports your main point in the introduction and the motivation of\nthe work is Edwards et al., Soft Matter 2023 (and in supplemental) for the flux pathway\niv) For some other systems, scientists must wait minutes after mixing to see phase\nseparation, such as poly(vinylpyridine) and poly(methacrylic acid) in Trinh et al.;\nSalvatore et al. Biomacromolecules 2011; Mohammadi et al. Journal of Colloid and\nInterface Science 2020; and Sanchez et al. Journal of Colloid and Interface Science\n2006.\nIn general, the Introduction and Discussion should account for an understanding of the\nresults in these and/or similar works at the authors' discretion.”\nReply:\nWe thank the reviewer for providing these highly relevant references.\nIn revision, we make the following modifications in the introduction section by\ndiscussing and citing the relevant work pointed by the reviewer (see page 3 in the main\ntext):\n“For example, many experiments trigger LLPS by applying environmental cues such\nas lights [D. S. Lee, N. S. Wingreen & C. P. Brangwynne, Biophys. J, (2021)] and pH\n[M. G. Last, S. Deshpande & C. Dekker, ACS Nano, (2020)] to concentrated, well-\nmixed polyelectrolyte solutions, which results in complete phase separation in a few\nminutes ([D. S. Lee, N. S. Wingreen & C. P. Brangwynne, Biophys. J, (2021)]; [D.\nBracha, et, al., Cell, (2018)]; [M. G. Last, S. Deshpande & C. Dekker, ACS Nano,\n(2020)]). This relatively rapid phase separation can be conceptually understood by the\nfacilitated contacts between oppositely charged polymers during mixing ([C. K. Trinh\n& W. Schnabel, Macromol. Chem. Phys., (1997)]; [X. Liu, J. P. Chapel & C. Schatz,\nAdv. Colloid Interface Sci., (2017)]; [X. Liu, et al., Soft Matter, (2016)]), while a\nquantitative examination of the underlying kinetic pathway and scaling law is not fully\nexplored. Recent work by Yuan and Tanaka...”\n“In marine organisms such as Sandcastle worms and mussels, even faster LLPS\ndynamics are observed, whereby these organisms produce proteins that rapidly form\nadhesive coacervates to attach to wet surfaces [R. J. Stewart, C. S. Wang & H. Shao,\nAdv. Colloid Interface Sci., (2011)]. In this scenario, oppositely charged glue proteins\nare secreted separately into “tubes” (e.g., byssal threads in mussels) at high\nconcentrations, and LLPS is initiated by fluxing these proteins into a shared region,\nwith the coacervates forming in as little as 30 seconds ([R. J. Stewart, C. S. Wang & H.\nShao, Adv. Colloid Interface Sci., (2011)]). Beyond biological systems, recent studies\nfound that applying shear flux to a mixture of polycation and polyanion solutions can\nsignificantly speed up their coarsening dynamics, where a clear phase separation is\nobserved on the order of seconds [C. E. Edwards, et al., Soft Matter, (2023)]. In\ncomparison, without the initial flux, the phase separation of the mixture can only be\nobserved after at least a few minutes ([C. K. Trinh & W. Schnabel, Macromol. Chem.\nPhys., (1997)]; [D. Salvatore, et al., Biomacromolecules, (2011)]; [P. Mohammadi, et\nal., J. Colloid Interface Sci., (2020)]; [C. Sanchez, et al., J. Colloid Interface Sci.,\n(2006)]). Clearly, different initial mixing protocols can significantly influence the phase\nseparation dynamics of dense polyelectrolyte mixtures. However, the fundamental\nmolecular mechanisms driving the dynamics of LLPS remain poorly understood. In\nparticular, the complete morphological pathways under varying initial mixing protocols,\nas well as their corresponding time evolution, are not yet fully elucidated.”\n6. “Experimental works describing precipitate versus coacervate growth in the same\ntypes of chemistries have observed diffusion show that precipitate growth resembles\ndiffusion-limited aggregation, instead of liquid-like growth signatures, but these types\nwould be indistinguishable using a structure-factor based measure of structure used in\nthis work. Still, given the need for more Discussion about growth rates, could the\nauthors comment more on which types of growth mechanisms could correspond to the\nobserved exponents?”\nReply:\nWe thank the reviewer for the suggestion to further discuss which of the growth\nmechanisms could correspond to the observed scaling exponents. Indeed, both the\ndiffusion-limited aggregation and liquid-like coarsening could lead to a scaling\nexponent of 𝑡(/&. In our simulations of the different pathways, we found that the small\ndroplets will move and coarsen with each other (see movie S1 in SI) at the later stage\nof the LLPS, with the formation of a single droplet at the end, demonstrating the domain\ngrowth is driven by liquid-like coarsening mechanisms. This liquid-like growth\nsignature is consistent with the results in our previous work [S. Chen & Z. -G. Wang,\nPRL, (2023)] and recent experimental observation of biomolecular condensate\nformations ([Y. Zhang, et al., PRL, (2021)]; [D. S. Lee, N. S. Wingreen & C. P.\nBrangwynne, Biophysical Journal, (2021)]).\nIn response to the reviewer’s suggestion, we clarify the liquid-like domain growth\nmechanism in Page 8 in the main text:\n“This droplet-coarsening feature is more clearly observed from movie S1 in the\nSupplemental Material, where small droplets move and coalesce with each other in the\nlater stage until the formation of a single large droplet. Our results show that the domain\ngrowth at later stage is mainly driven by droplet-coarsening instead of the diffusion-\nlimited aggregation (DLA) of solid-like precipitates, consistent with recent simulation\nresults [S. Chen & Z. -G. Wang, PRL, (2023)] and recent experimental observation of\nbiomolecular condensate formations ([Y. Zhang, et al., PRL, (2021)]; [D. S. Lee, N. S.\nWingreen & C. P. Brangwynne, Biophysical Journal, (2021)]).”\n7. “The claim made at the end of the paper \"Additionally, ultrafast coacervation can\noccur in the flux pathway, driven by electrochemical potential gradients from the\ninitially separated polycation and polyanion domains, mimicking the rapid adhesive\nprocesses observed in marine organisms\" is quite strong, when compared to the what\nwas found in the Results, namely that \"The slightly faster scaling [0.39 exponent]\nthan the expected t^1/3 is probably due to the residual electrochemical driving force for\nthe polyions to migrate towards the interfacial region, as schematically shown ...\". It would\nmake the paper stronger to discuss this result more. Why is it \"ultrafast\" in the context\nof literature for theory, simulation, and/or experiment? What is the meaning of the\nt^(1/3) exponent? You should discuss LSW theory/Ostwald Ripening, and diffusion-\nlimited growth, and explain the 0.39 exponent within this context; and whether/what\nfaster exponents have been previously observed in experiment, simulation, or estimated\nin theory. Also, it is important to clearly differentiate in the Conclusions what is\nsuggested by the results (e.g., electrochemical potential drives ultrafast phase\nseparation), versus what you have directly observed in them (an increase in domain\ngrowth exponent to 0.39 versus the 0.33).”\nReply:\nWe thank the reviewer for the insightful suggestion to provide a deeper discussion on\nthe accelerated growth dynamics regarding the flux pathway. In the flux pathway, a\nclear 𝑡*/&scaling occurs at the early stage of the LLPS, and this scaling persists longer\nas the system size increases as we discussed in our manuscript (see Fig R1(b) in this\nresponse letter and Fig. S5 in SI). As such, we believe that in real physical systems\nwhere the separated polymer solution domain is much larger, the LLPS dynamics are\nsignificantly accelerated by the persisted 𝑡*/&scaling. This is what we meant by the\nultrafast dynamics.\nThe 0.39 slope at the later stage of the flux pathway will occur much later in real\nsystems. Indeed, both the Lifshitz-Slyozov-Wagner (LSW) theory on Ostwald ripening\nand diffusion-limited growth predict that the later stage of phase separation dynamics\nfollows 𝑅~𝑡(/& . This theoretical framework does not involve the long-range\nelectrochemical potential difference relevant in the flux pathway. The residual\nelectrochemical potential gradient is expected to accelerate LLPS dynamics, leading to\nan intriguing exponent of around 0.39. However, we do not find any experimental\nstudies reporting these accelerations. It would be highly desirable for future\nexperiments to validate our predictions.\nIn response to the reviewer’s suggestion, we add the following modifications in the\ndiscussion section in the main text (see page 13&14):\n“Interestingly, the 𝑡*/& scaling observed in the early stage of the flux pathway persists\nmuch longer in a larger system, as shown in Fig. S5. We thus believe that in real, larger\nsystems where the separated polymer solution domain is much larger, the 𝑡*/& scaling\nwill persist to macroscopic length scale and significantly accelerate the coarsening;\nwe refer to this coarsening process as ultrafast LLPS dynamics. In later stage, the 𝑡$.&,\nscaling is slightly larger than the classical 𝑡(/& scaling predicted by the Lifshitz-\nSlyozov-Wagner (LSW) theory on Ostwald ripening or diffusion-limited growth. This\nsmall acceleration is probably due to the residual electrochemical gradient, which\nprovides additional driving force to facilitate the migration of polyions towards the\ngrowing domains. It would be highly desirable for future experiments to validate these\npredictions.”\n8. “Sentence pg. 12: “On the other hand, if the polymer concentration is significantly\nlower ... but still larger than the coexisting supernatant concentration), the polyions are\nexpected to form charge-balanced pairs first, and the coarsening dynamics should\nfollow R ~ t^(1/3)”. But, this claim was not investigated. It would be nice to see it using\na simulation, or is there a reason you did not include this data?\nReply:\nWe thank the reviewer for the valuable suggestion to provide further data to support\nour claim about the coarsening dynamics at lower polymer concentration. Indeed, this\nstatement is a theoretical expectation based on the well-established understanding from\nour previous work ([S. Chen & Z.-G. Wang, PRL, (2023)]; [S. Chen & Z.-G. Wang,\nPNAS, (2022)]) and other studies on dilute, opposite-charged, polyelectrolyte systems\n[K. I. Tainaka, J. Phys. Soc. Jpn., (1979); D. Priftis, et al., Langmuir, (2012)]. In the\nvery dilute polyelectrolyte solution containing a few polycations and polyanions, the\nstrong electrostatic attraction promotes the initial formation of charge-neutral polyion\npairs ([S. Chen, P. Zhang & Z.-G. Wang, Macromolecules, (2022)]; [S. Chen & Z.-G.\nWang, PNAS, (2022)]; [S. Adhikari, M. A. Leaf & M. Muthukumar, JCP, (2018); [K.\nT. Delaney & G. H. Fredrickson, JCP, 2017]) Starting with many of these small neutral\npairs, the coarsening dynamics will follow the classical 𝑡(/&scaling, as demonstrated in\nour previous work [S. Chen & Z.-G. Wang, PRL, (2023]].\nHowever, direct simulation of LLPS starting with dilute, dispersed chains is not\ncomputationally feasible, as the corresponding scaling law of LLPS will require an\nextremely large simulation box to contain a reasonable amount of polyion pairs solvated\nin a significant amount of water.\nIn revision, we clarify this by adding the following text to page 12 in the main text.\n“However, direct simulation of LLPS of dilute systems starting with dispersed chains\nis not computationally feasible within our GCMe-DPD framework, as it would require\na much larger system to obtain a statistically reliable domain growth law.”\ni) At the overlap concentration, can you quantify how much the droplets are charge-\nasymmetric, and how much this “further slows” coarsening?”\nReply:\nIn response to the reviewer’s question, we performed a quantitative analysis of the\ncharge asymmetry within single droplets at different polymer concentrations. Figure\nR2(a) shows that the droplets exhibit a larger degree of charge asymmetry at\noverlapping polymer concentration (𝜙 = 1.0𝜙∗), compared to the case of 𝜙 = 2.0𝜙∗.\nThe charge asymmetry slows down the coarsening dynamics, as shown in Fig. R2(b)\nand Fig. 4(a) in the main text.\nIn revision, we add the above discussion and Fig. R2(a) into page S9 of SI.\nFigure R2: (a) The distribution of the number of droplets with different charge\nasymmetry degrees 𝑎 in the systems with 1.0𝜙∗ and 2.0𝜙∗. Here, 𝑎 is defined as 𝑎 =\n./0 (3!,3\")\n, where 𝑛 and 𝑛 are respectively the number of polycation and polyanion\n( *\n.67(3!,3\" )\nbeads in a single droplet. (b) Time evolution of the domain growth in the early stage of\nthe randomly distributed well-mixed system with concentration 1.0𝜙∗, 2.0𝜙∗, and\n2.7𝜙∗, with chain length 𝑁 = 40.\nMinor Comments:\n1. First sentence of the Results is awkward, rephrase.\nReply:\nWe thank the reviewer for pointing this out.\nIn the updated manuscript, we have rephrased it as: “In this study, we employ molecular\ndynamics simulations that explicitly incorporate electrostatics and hydrodynamics to\ninvestigate the kinetic pathways and LLPS dynamics of polyelectrolyte complex\ncoacervation under several representative initial mixing conditions.”\n2. Lambda values are only labelled on some of the plots.\nReply:\nIn the updated manuscript, we have labeled all plots in the well-mixed pathway with\ntheir corresponding 𝜆 values to consistently show the initial mixing degree (See Fig. 2\nin the main text).\n3. The letter labels (a), (b), (c), (d) in the Figures are often in a strange location, e.. in\nthe middle of each subplot.\nReply:\nIn the updated manuscript, we have now repositioned all subplot labels (a), (b), (c), etc.,\nto the top-left corner of the pictures throughout the manuscript to ensure a uniform data\npresentation.\n4. Please double-check the work for sentence-level edits. For example, a grammatical\nerror in the very 1st sentence of abstract can obfuscate your meaning. Check throughout.\ni) In addition, 1st sentence pg. 3: How does a “thermodynamic pathway ... model the\ncoarsening dynamics”? Consider rewriting to clarify what point you are making here.\nii) The proper phrasing in Acknowledgements is “We thank X for their generous\nallocation of Y”.\nReply:\nWe thank the reviewer for carefully reviewing our manuscript and raising these\ngrammatical errors in our main text.\nIn the updated manuscript, we have rephrased these sentences in the abstract,\nintroduction, and acknowledgement.\n5. In Figure S3 and related discussion, make sure to acknowledge Figure S4.\nReply:\nIn the updated manuscript, we have referenced Fig. S5 (the original Fig. S4) to discuss\nthe system size dependence of the early-stage scaling regime in both the well-mixed\nand flux pathway.\nResponses to the Comments of Reviewer #3\n“Manuscript Summary and Key Findings\nIn this manuscript, the authors present a study on the effects of mixing protocols on the\ndynamics of liquid-liquid phase separation of polyelectrolyte complex coacervates. The\nauthors use molecular dynamics simulations to evaluate the domain growth of the\ncoacervate droplets for two distinct initial mixing conditions. The first is a well-mixed\nsolution of oppositely charged polyelectrolytes. The second scenario involves a system\nwhere the polycations and polyanions are separated into domains that are put in contact\nwith each other, the authors called this scenario the flux pathway.\nThe authors found that the two initial mixing scenarios produced different domain\ngrowth behaviors. The first scenario, the well mixed pathway, showed the formation of\na network structure in the early stages that would then break down into droplets which\neventually merge into a single droplet. In the flux pathway, coacervate droplets will\nform at the interface between the two domains, which will then merge and draw in free\nchains to grow the domains. The well mixed pathway showed early-stage growth with\nscaling of t1/2, and the flux pathway showed an early-stage domain growth with scaling\nt2/3.\nThe authors also evaluated the polymer concentration effects on the domain growth\nmechanism and dynamics and found that the network structure that leads to enhanced\nearly-stage growth requires a sufficiently high concentration to observe the enhanced\ngrowth dynamics. The authors showed that the initial mixing conditions can be used to\naccelerate coacervation dynamics, and the implications of system size are discussed for\nthe various initial mixing conditions.\nReview\nThis manuscript provides valuable insight into the dynamics of polyelectrolyte\ncoacervation that will be of considerable interest to the coacervate and polyelectrolyte\ncomplex communities. The findings in this paper are, to my knowledge, novel and\nimportant. I believe the methods used by the authors are valid for the study, and the\nfindings presented within the manuscript are well reasoned and supported by the data.\nThe manuscript is clear and concise, and provides the context necessary to understand\nthe significance and motivations of the study. The references cited in the manuscript\nare appropriate and many are recent. The information provided in the supplementary\nmaterials is appropriately placed as supplemental and provides more details regarding\nthe simulation methods and data analysis. The methods have been sufficiently described\nto allow for reproduction of the study by other researchers.\nThe modeling and simulation methods used in this study was developed and previously\npublished by the corresponding authors, and addresses some key issues known within\ncommunity. Molecular dynamics simulations of complex coacervates can be\ncomplicated by the length and timescale requirements to observe phase separation,\nas well as the nature of the hard cores of the Lennard Jones potential and the treatment of\nbulky ions as point charges within the center of the ion bead. The Gaussian Core Model\nwith Smeared Electrostatic Interactions methodology presented by the authors\naddresses these three potential issues while maintaining the ability to evaluate dynamics\nand consider hydrodynamic interactions. This methodology could be used by other\nresearchers studying coacervation and other solutions of charged macromolecules,\nparticularly for out of equilibrium and dynamics studies.\nIt is my belief that the data from which the scaling exponents are drawn from are high\nquality and robust. The authors in this manuscript provide well-reasoned explanations\nwith mathematical descriptions for some of the findings in the study. For example, they\nfound two different late-stage domain growth dynamics for the well mixed pathway\nthat they attributed to the initial degree of mixing. The authors provided a method for\nevaluating the degree of mixing and found that indeed the value for this degree of\nmixing correlated with the late-stage domain growth dynamics. Another example of the\nproviding rigorous explanations for their findings can be found in their explanation for\nthe driving force that leads to the t2/3 domain growth scaling observed in the flux\npathway. The authors mathematically described the electrochemical potential and\nshowed how the potential difference between the interface and the two initial domains\ndecreases over time to lead to a reduction in the domain growth dynamics between the\nearly stage and later stage of the flux pathway coacervation.\nThe authors’ discussions and findings regarding the concentration dependence on early-\nstage domain growth dynamics and the system size implications are very insightful. For\nexample, the authors extrapolated the scaling for macroscopically large systems, and\nfound that the flux pathway can undergo coacervation in a matter of seconds, but found\nthat the baseline thermodynamic pathway would take 47.5 years to form a coacervate\ndomain of the same size.\nIt is not entirely clear to me the level of significance the broader audience of Nature\nCommunications would find in this manuscript, though I believe this topic is of\nconsiderable interest to the soft materials scientific community. Coacervates represent\na potential environmentally friendly and sustainable solution for many materials\napplications, where these characteristics are becoming increasingly sought after. I\nbelieve this manuscript is of high significance and high quality and should be published\nin Nature Communications. I have a few questions for the authors to possibly consider\ndiscussing within the manuscript.”\nWe thank the reviewer for carefully reading our manuscript, and we are grateful to\nknow that the reviewer found our work of high significance and would support the\npublication of our manuscript in Nature Communications. Below, we provide a point-\nby-point response to the comments.\nQuestions for the authors\n“1. One key finding in this study is the enhanced domain growth as a result of network\nformation. This network formation is concentration dependent, and I am wondering\nabout the implications of the inhomogeneous polarizability throughout the system that\nwill be more prevalent at higher concentrations. How does the model account for these\nlocal changes in polarizability between water dominant regions and polymer dominant\nregions?”\nReply:\nWe thank the reviewer for raising this important issue, which is indeed common to most\ntheories and essentially all coarse-grained molecular dynamics simulations of polymer\nsolution systems that treat the solvent as a dielectric continuum. As the reviewer\ncorrectly points out, the inhomogeneous polarizability throughout the system can\nbecome more pronounced at higher concentrations, which requires a better theory to\naccount for the local polarization. While computational methods for including this local\nchange in polarizability do exist ([D. Grzetic, et al. ACS Macro Lett., (2019)]; [A.\nGavrilov., JCP, (2020)]; [A. Gavrilov & E. Kramarenko, JCP, (2021)]), accounting for\nthis effect within the context of polyelectrolyte complex coacervation presents\nsignificant computational challenges and is beyond the scope of this work. We admit\nthe limitation of the dielectric continuum approach. However, we note the general\nsuccess of this dielectric continuum approach for capturing electrolyte and\npolyelectrolyte systems even at fairly high salt concentrations [R. Triolo, J. Grigera &\nL. Blum, J. Phys. Chem., (1976)], as well as its success in capturing the network\nformation in LLPS [Yuan & Tanaka, Nat. Commun, (2025)]. As such, we expect our\ncentral conclusions on the network-formation dynamics to hold.\nIn response, we added the following discussion to the Discussion section (See page 14\nin the main text):\n“Furthermore, we use a homogeneous dielectric constant to model electrostatics, which\ndoes not account for possible local dielectric inhomogeneity arising from polymer\naggregations. Accounting for this local dielectric effect requires better theoretical\nunderstanding and new methodology development in the construction of coarse-grained\nmodels. In principle, all these effects can be captured by all-atom simulations of specific\nsystems. However, such calculations are at present computationally prohibitive for\nexploring the long-time behavior of large coacervate systems and are well beyond the\nscope of this work. As this work concerns the scaling behavior in LLPS dynamics of\ngeneric coacervate systems on length and time scales well beyond the residue level, we\nbelieve that incorporating atomistic details will not qualitatively change our\nconclusions.”\n“2. What are the implications of using the same bead size for the salt counterions as the\nother species, on the dynamics of coacervation? How accurate are the salt ion mobilities\ncaptured in this model and could deviations to the salt ion mobilities lead to significant\nchanges in the mechanisms and dynamics of coacervate domain growth?”\nReply:\nThe reviewer has raised an excellent question about the implications of the salt\ncounterion mobilities for our results. We agree that the simplified setup (commonly\nused in coarse-grained simulations) of using the same bead size for the salt ions as other\nspecies does not account for the corresponding specific ion effects. While the coarse-\ngrained approach may not capture the precise mobility of specific ions, it can capture\nthe difference (on the order of magnitude) in length and time scales between small ions\nand polymer aggregates. In our systems, the length and time scales for polymer phase\nseparation are well beyond the ion size and diffusion time. For example, Figure. R3\nshows the mean-squared-displacement (MSD) of an ion with size 𝜎 reaches the free\ndiffusive regime at the time on the order of 10𝜏 , while the relevant timescale of the\nLLPS dynamics investigated in this work is on the order of 10&𝜏 ~108𝜏 . As such, we\nbelieve the precise value of ion mobility is inconsequential to the qualitative behavior\nin the LLPS dynamics.\nIn revision, we add the following discussion to the SI section (See page S3 in the SI):\n“The use of the same bead size for the ions as for other species may not represent the\nprecise mobility of specific ions. However, we expect the salt ion mobilities will not\nfundamentally change the polymer LLPS pathways and dynamics, as the corresponding\nlength and time scales are well beyond the ion size and diffusion time.”\nFigure R3: The mean square displacement ⟨𝑟#(∆𝑡)⟩ of the counterions in the\nthermodynamic pathway.\n“3. It seems to me that the interface surface tension between the coacervate domain and\nthe solvent and free chain domain is important for the dynamics of droplet coarsening.\nIs this the case and how accurately is this interface surface tension being modeled with\nthis coarse-grained model? Would you expect the structuring of water around the\ncoacervate droplets to have an impact on the droplet coarsening dynamics?”\nReply:\nThis is an interesting question regarding the influence of interfacial tension on the\ndroplet coarsening dynamics. As the reviewer correctly points out, the interfacial\ntension between the coacervate phase and the supernatant is a major driving force for\nthe coarsening dynamics. Theoretical and experimental studies have demonstrated that\nthe surface tension of the electrostatic-driven coacervates is very low ([J. Qin, et al.,\nACS Macro Lett., (2014)]; [P. Zhang & Z. -G. Wang, Macromolecules, (2021)]; [D. J.\nAudus, et al., PRL, (2021)]). Such low surface tension is reasonably captured by coarse-\ngrained simulations ([J. Yuan & H. Tanaka, Nat. Commun., (2025)]; [Y. N. Fang, et al.,\nProc. Natl. Acad. Sci., (2023)]; [B. Yu, et al., Macromolecules, (2023)]; [B. Yu, et al.,\nMacromolecules, (2022)]; [D. Sharon, et al., Macromolecules, (2022)]; [H. Liang, et\nal., Macromolecules, (2022)]) similar to ours. As such, we believe our systems\nrepresent the correct coacervation dynamics driven by electrostatic correlations of\npolyelectrolytes.\nRegarding the water structure around the coacervate droplets, we expect this would\ninfluence the droplet coarsening dynamics by changing their induced electrostatic\npotential at surface ([L. Shi et al., Nature, (2025)]; [S. Pullanchery et al., Science,\n(2021)]; [J. K. Beattie & A. M. Djerdjev, Angew. Chem. Int. Edn, (2004)]). However,\nexploring this effect is beyond the scope of this work for the following two reasons: 1)\ncapturing the characteristics of water structure would require all-atom simulations,\nwhich is nearly computationally unaffordable for simulating the LLPS dynamics with\nthe large length and time scale in this study; 2) The underlying water orientational\nordering at oil-water interface and the induced surface electrostatic potential is currently\nunder debate and will require very accurate polarizable potentials ([S. Pullanchery,\net al., Science, (2021)] ;[Zhao et al., Angew. Chem. Int. Edn, (2025)]).\nIn response, we add the following modification to the main text and SI.\n“We note that the droplets in our simulations generally have floppy surfaces, indicating\nthey have low surface tension, which is consistent with previous theoretical and\nexperimental studies in coacervate systems ([J. Qin, ACS Macro Lett., (2014)]; [P.\nZhang & Z.-G. Wang, Macromolecules, (2021)]; [D. J. Audus, et al., PRL, (2021)]).”\n–– see page 5&6 in the main text.\n“Besides, our model does not capture the possible water orientational order and the\ncorresponding induced electrostatic potential at coacervate surfaces ([L. Shi et al.,\nNature, (2025)]; [S. Pullanchery et al., Science, (2021)]; [J. K. Beattie & A. M. Djerdjev,\nAngew. Chem. Int. Edn, (2004)]), which could influence the coarsening dynamics. ” –\n– see page 14 in main text.\n“4. I am wondering about the domain size calculations from the structure factor and the\nerror in the domain size as a result of noise in the structure factor profiles. Do you have\na sense for how sensitive the domain size calculation is to the choice of q-vector binning?\nPresumably, larger systems would allow for obtaining structure factor profiles with less\nnoise, and so I am wondering if the systems studied are sufficiently large to get accurate\ndomain sizes from the structure factor.”\nReply:\nThis is an interesting technical question. In principle, all the q-vectors used in structure\nfactor calculations should be consistent with the system dimensions, which is defined\nas (𝑛 ∙2𝜋/𝑙 , 𝑛 ∙2𝜋/𝑙 , 𝑛 ∙2𝜋/𝑙 ) , where 𝑛 , 𝑛 , 𝑛 are integers and 𝑙 , 𝑙 , 𝑙 are\nthe box dimensions. In data analysis, the contributions of these wavevectors are usually\nsummed to a reduced wavenumber with the same magnitude q. In our work, we use\nthe minimal possible q-bin width with ∆𝑞 = 2𝜋/𝑙 (since the corresponding length\nscale of ∆𝑞 cannot exceed the simulation box size), to maximize the resolution of the\nstructure factor. As such, our calculations already have the highest possible resolution\nof S(q). Here, we also compare the case when we increase the bin width by a factor of\n2. Figure R4 shows that the domain size and scaling law in the thermodynamics\npathway remain quantitatively unchanged. This demonstrates that our system (has the\nsmallest size of q-bin) is sufficiently large to get statistically reliable results.\nIn response, we have added the following clarification to page S6 in SI:\n“The q-space is defined as (𝑛 ∙2𝜋/𝑙 , 𝑛 ∙2𝜋/𝑙 , 𝑛 ∙2𝜋/𝑙 ) , where 𝑛 , 𝑛 , 𝑛 are\nintegers and 𝑙 , 𝑙 , 𝑙 are the box dimensions. In our data analysis, we employed the\nminimal possible q-bin width with ∆𝑞 = 2𝜋/𝑙 , to maximize the resolution of the\nstructure factor. We have shown that within the uncertainties of the error bar, increasing\nthe q-bin size by a factor of 2 has no effect on the domain growth behavior; see Fig. S3”\nFigure R4: Time evolution of the domain growth that obtained from different q-vector\nbinning strategy with ∆𝑞 =2𝜋/𝑙 and 4𝜋/𝑙 , respectively, in an individual simulation\nof thermodynamic pathway."}, {"round": 2, "reviewer_comments": {"reviewer_3": "I'd previously recommended the manuscript for publication and I continue to recommend publication of the manuscript. The authors answered my questions and remarked on these answers within the manuscript. I believe that the authors have also satisfactorily responded to the questions and points brought up by the other reviewers.\n(Remarks on code availability)\nThe authors provide the simulation files needed to run the simulations and they provide notes to aid someone in attempting to run the code. I was not able to run the code as I do not have LAMMPS installed, but the simulation files appear to be complete. The authors provide instructions for implementing the custom pair style needed to implement the simulation model and an example simulation script file with the input structure file."}, "author_response": ""}]

Decision letter