1 Introduction↩︎

Polymers play a crucial role in many technological and biological applications [1]–[3]. A single polymer exhibits complex dynamics that span multiple time and length scales [4]–[9], [9], [10], [10]–[16]. Multiple polymers put together can give rise to entanglements, leading to dynamics such as reptation [16]–[18]. Mixing two different polymeric species considerably complicates the dynamics, which subsequently depend on the underlying interactions between the polymer types, their degree of entanglements, as well as their mechanical properties. In such a problem, the polymers can phase-separate, or diffuse into one another — a scenario termed polymer-polymer interdiffusion [4]–[9], [9], [10], [10]–[15]. Below, we study the effects of entanglement on the spatio-temporal evolution of polymeric droplets undergoing interdiffusion.

Besides their technological value, concepts from polymer physics can help in understanding many biological systems. One of the predominant current frameworks to understand the formation of biological condensates is via liquid-liquid phase separation, which utilizes a similar framework as for polymer phase separation [19]–[23]. For example, phase separation in polymer-fluid mixtures exhibits biologically relevant morphologies, such as gel networks and condensates [24]–[31]. Consequently, the application of ideas from polymer science to biological scenarios has motivated measurements of mechanical properties of biological condensates, improving our understanding of their rheological responses [32]–[35], [35].

While passive processes such as diffusion and phase separation are ubiquitous, many biological processes involve active production of polymers, for example in transcription and translation. One prominent example of such an active source is the nucleolus — a subnuclear condensate synthesizing and exporting ribosomal RNA (rRNA) [36]–[40]. The nucleolus acts as a persistent source, driving an outward flux of biopolymers into the surrounding nucleoplasm. Clearly, once such a source is present, the polymer concentration (of rRNAs in the above example) is also affected by the presence of the source.

Here, inspired by such biological condensates, we use a two-fluid formalism [41]–[50] to study two types of polymer-polymer interdiffusion problems — one of passive droplet diffusion, and another of source-driven diffusion, both depicted schematically in Fig. 1 (a)-(b) respectively. We first review in Sec. 2 the equations of motion of the two-fluid formalism [41]–[50], which yield a reduced set of one-dimensional (1D) equations of motion. Then we analyze the passive droplet diffusion cases in Sec. 3. We show that under certain approximations, the polymer droplet interdiffusion is equivalent to a nonlinear diffusion problem, allowing us to cast the equations of motion in a self-similar form. We use the 1D formalism to obtain full numerical solutions, and compare those with the self-similar approach. We also demonstrate that our results hold in a two-dimensional (2D) case.

Then, we introduce a source and analyze the emerging scaling, the resulting reduced equations of motion, and the corresponding dynamics in Sec. 4. Specifically, we demonstrate that the source introduces an additional flux term, and breaks the self-similar structure of the passive problem. Consequently, instead of using a self-similar ansatz, we perform a perturbation around the expanding droplet’s tip, and show that the local geometry at the front of a source-driven droplet resembles the one observed in the passive diffusion scenario. We compare our analytical predictions against the 1D and 2D full numerical solutions, demonstrating the spatial similarity to the passive case. Finally, we discuss our results and potential implications in Sec. 5.

2 Equations of motion↩︎

We consider a two-component mixture, of polymeric species A and B, referred to as A-mers and B-mers, both of the same density [41]–[50]. With the A-mer and B-mer velocity fields \(\boldsymbol{u}_A(\boldsymbol{r},t)\) and \(\boldsymbol{u}_B(\boldsymbol{r},t)\), respectively [below we suppress the explicit \((\boldsymbol{r},t)\) dependence for readability], the volume fractions \(\phi_A\) and \(\phi_B\) evolve via the continuity equations \[\label{eq:con95clean} \begin{align} \partial_t \phi_A + \nabla\cdot\left(\boldsymbol{u}_A \phi_A\right) &= S \;, \\ \partial_t \phi_B + \nabla\cdot\left(\boldsymbol{u}_B \phi_B\right) &= 0, \end{align}\tag{1}\] where here we introduced an A-mer source \(S\!\equiv\!S\left(\boldsymbol{r},t\right)\). Because we consider a binary mixture, \(\phi_A + \phi_B\!=\!1\), and we denote \(\phi_A\!=\!\phi\).

The quasi-static momentum conservation balances the pressure \(p\), the osmotic pressure \(\boldsymbol{\Pi}\), and the material stresses \(\boldsymbol{\sigma}_{A}\) and \(\boldsymbol{\sigma}_{B}\) as [49], [50] \[{\boldsymbol{0}} = \nabla\cdot\left(-p\boldsymbol{I}-\boldsymbol{\Pi} + \boldsymbol{\sigma}_{A} + \boldsymbol{\sigma}_{B}\right) \;, \label{eq:qs95clean}\tag{2}\] where \(\boldsymbol{I}\) is the identity tensor.

The motion of one component immediately implies the relative motion of the other component. This relative motion is obtained as [42], [45], [49] \[\label{eq:rel95clean} \boldsymbol{u}_A - \boldsymbol{u}_B = -\frac{1}{ n_0\zeta_0\mu}\left[\left(1-\phi\right)\nabla\cdot\left(\boldsymbol{\Pi}-\boldsymbol{\sigma}_{A}\right)+\phi\nabla\cdot\boldsymbol{\sigma}_{B}\right] \;,\tag{3}\] where \(n_0\) is the mixture’s number density (monomers per volume), \(\zeta_0\) is a monomer friction coefficient (assumed to be the same for both A-mers and B-mers), and \(\mu\!\equiv\!\mu\left(\phi\right)\) is a dimensionless friction function that depends on the polymers and their surroundings [9], [10], [42]. The functional form of \(\mu\) is expressed via the friction function of the A-mers and B-mers, \(\mu_A\) and \(\mu_B\) respectively, as \[\label{eq:mu} \mu \!=\! \frac{\mu_A \mu_B}{\mu_A + \mu_B} \;.\tag{4}\] The individual friction functions \(\mu_A\) and \(\mu_B\) depend crucially on the polymer’s lengths, \(N_A\) and \(N_B\), compared to an entanglement threshold \(N_e\). Polymers shorter than \(N_e\) are considered unentangled, implying \(\mu_{(i)}\!=\!\phi_{(i)}\) (for species \(i\)). Entangled polymers are subjected to increased friction, \(\mu_{(i)}\!=\!\alpha_{(i)}\phi_{(i)}\), where \(\alpha_{(i)}\!\equiv\!N_{(i)}/N_e\!>\!1\) [9], [10].

The osmotic pressure \(\boldsymbol{\Pi}\) originates from the polymer concentration, i.e., chemical potential, gradients. We utilize the Flory-Huggins (FH) thermodynamic potential [51], [52] \(F_{FH}\!=\!k_B T n_0 f(\phi)\), where \(k_B\) is the Boltzmann constant, \(T\) is the temperature, and \[\label{eq:FH} f(\phi) = \frac{\phi}{N_A} \log \phi + \frac{\left(1-\phi\right)}{N_B} \log\left(1-\phi\right) + \chi \phi \left(1-\phi\right) \;.\tag{5}\] Here, \(N_A\) and \(N_B\) are the number of monomers per A-mer and B-mer respectively, and \(\chi\) is an energy interaction parameter between the two species. As we are interested in the mixing of the two species, we omit surface tension. For simplicity, we consider non-interacting polymers, setting \(\chi\!=\!0\). For the FH potential of Eq. 5 , the resulting osmotic pressure is isotropic, obtained as \(\boldsymbol{\Pi}\!=\!k_B T n_0\boldsymbol{I} \left(\phi\partial_\phi f-f\right)\) [1], [42]. Note that the osmotic pressure divergence \(\nabla\cdot\boldsymbol{\Pi}\) is equivalent to the spatial gradients of the chemical potential, as \(\!\nabla\left[\phi\partial_\phi f-f\right]\!=\phi \nabla\left(\partial_{\phi} f\right)\) (we use \(\boldsymbol{\Pi}\) to be consistent with the existing literature). We plot \(f(\phi)\) for different values of \(N_A\) (in the non-interacting limit \(\chi\!=\!0\), and for monomer surroundings, \(N_B\!=\!1\)) in Fig. 1 (c).

The stresses \(\boldsymbol{\sigma}_{A}\) and \(\boldsymbol{\sigma}_{B}\) evolve via an upper-convected Maxwell (UCM) model [49], [53], \[\label{eq:UCM} \stackrel{\smalltriangledown}{\boldsymbol{\sigma}}_{(i)} = G_{(i)}\left(\boldsymbol{L}_{(i)} + \boldsymbol{L}_{(i)}^T\right) - \frac{1}{\lambda_{(i)}}\boldsymbol{\sigma}_{(i)} \;,\tag{6}\] where \(i\) denotes either A or B (no summation implied), \(\stackrel{\smalltriangledown}{\bullet}\equiv\partial_t \bullet + \boldsymbol{u}_{(i)} \cdot \nabla \bullet -\left[\boldsymbol{L}_{(i)}^T\cdot \bullet + \bullet \cdot \boldsymbol{L}_{(i)}\right]\) is the upper-convected time derivative, \(G_{(i)}\) is the shear modulus, \(\lambda_{(i)}\) is the relaxation time, and \(\boldsymbol{L}_{(i)}\!\equiv\!\nabla \boldsymbol{u}_{(i)}\) is the velocity gradient tensor for the \(i\)th component. While below we focus on a regime where the material stresses \(\boldsymbol{\sigma}_A\) and \(\boldsymbol{\sigma}_B\) could be safely neglected, we retain those for completeness.

To ease the analysis of the above equations of motion, we rescale Eqs. 1 6 with \(\lambda_A\) as a time scale, \(\sqrt{\frac{k_B T \lambda_A}{\zeta_0}}\) as a length scale, and \(n_0 k_B T\) as an energy density scale. Thus, we arrive at a unitless, 1-dimensional (1D) reduced form of the above equations, cast as (\(i=A,B\)) \[\begin{gather} \partial_t \phi + \partial_x\left(u_A \phi\right) = S \;, \tag{7}\\ -\partial_t \phi + \partial_x\left[u_B \left(1-\phi\right)\right] = 0 \;, \tag{8} \\ \partial_x\left(-p-\Pi + \sigma_{A} + \sigma_{B}\right) = 0 \;, \tag{9}\\ u_A - u_B = -\frac{1}{ \mu}\left[\left(1-\phi\right)\partial_x\left(\Pi-\sigma_{A}\right)+\phi\partial_x\sigma_{B}\right] \;, \tag{10}\\ \partial_t \sigma_{(i)} + u_{(i)} \partial_x \sigma_{(i)} = 2\left(\sigma_{(i)}+G\right)\partial_{x}u_{(i)}- \frac{1}{\lambda_{(i)}}\sigma_{(i)} \tag{11}\;. \end{gather}\]

Below we study the passive diffusive dynamics, and source-driven diffusive dynamics, of asymmetric polymer mixtures, where we assume \(N_A\!\gg\!N_B\). We also assume that \(N_B\) is less than the entanglement threshold \(N_e\), implying the B-mers are not entangled throughout the dynamics. To allow the A-mers to transition from a possibly entangled state (if they are sufficiently long) to an unentangled state, we use \(\mu_A\!=\!\phi\left[\alpha \phi + \left(1-\phi\right)\right]\) as the A-mers friction function. This function reduces to \(\mu_A\!=\!\phi\) in the unentangled limit (when \(\alpha\!=\!1\)). If the polymers are sufficiently long to entangle, it allows for a transition \(\mu_A\!\rightarrow\!\phi\) at low volume fractions, and \(\mu_A\!\rightarrow\!\alpha\) as \(\phi\!\rightarrow\!1\). We plot three examples of \(\mu(\phi)\) for different \(\alpha\) values in Fig. 1 (d). As \(\alpha\) increases, \(\mu\) develops an asymmetry, increasing the friction at intermediate \(\phi\) regions and approaches \(1-\phi\) — an effect that is due to entanglements.

As we would like to focus on the role of entanglement and the passive or source-driven cases, we will neglect the roles of the material stresses \(\boldsymbol{\sigma}_A\) and \(\boldsymbol{\sigma}_B\). To estimate their relative importance, we note that polymeric stresses \(\boldsymbol{\sigma}_{(i)}\) are usually characterized by the elastic moduli \(G_{(i)}\), and relaxation times \(\lambda_{(i)}\), as explained above [45], [49], [53] [see also Eq.@eq:eq:UCM ]. While precise measurements are challenging, biological condensates yield maximal storage moduli of \(\sim\!1\!-\!10 \;\text{Pa}\) at the low-frequency limit, and maximal viscosity of \(\sim\!10\!-\!100 \;\text{Pa} \cdot \text{s}\) [32]–[35], [35]. Comparing these estimates with the energy density scale \(n_0 k_B T \simeq \text{KPa}\) at \(T\!=\!300~K\) (using \(n_0\!=\!1 ~\text{mM}\) as a low estimate for the monomer density), clearly indicates that under typical conditions the osmotic pressure is the dominant source of stress.

The relaxation times of biopolymers may also be sufficiently fast for our interests. Biological polymers exhibit relaxation times of \(\lambda\!\simeq\!10 \;\text{ms}\) [32], [35]. If we compare this to the expected (linear) diffusion time \(\tau_D\!\simeq\!R^2 \zeta_0 / k_B T\) where \(R\) is a typical droplet radius \(R\!\simeq\!5 \;\mu \text{m}\), and \(\zeta_0\!\simeq\!6\pi\eta a\) with \(\eta\!=\!1 \;\text{mPa}\cdot\text{s}\) and \(a\!\simeq\!1 \;\text{nm}\) being the monomer (protein) typical size, we get \(\tau_D\!\simeq\!100 \;\text{ms}\), again confirming that polymer relaxation occurs at shorter times than those of interest here. Interestingly, a reminiscent time-scale separation was observed in [35] for trapped droplets composed of different biopolymers and probed with optical tweezers. These separations of stress scales and time scales allow us to confine our analysis to cases in which these two assumptions hold, and the material stresses \(\boldsymbol{\sigma}_A\) and \(\boldsymbol{\sigma}_B\) can be safely neglected. While velocity gradients and polymer conformations near interfaces could induce mechanical stresses, considering the osmotic pressure captures the dominant driving force needed to describe the emerging nonlinear transport. As we show below, the simplifying assumptions above allow us to study polymer-polymer interdiffusion and its dependence on two central aspects — the degree of entanglement and the absence or presence of a source.

3 Droplet diffusion↩︎

We first consider a droplet of A-mers placed in a B-mers medium, depicted schematically in Fig. 1 (a). As the species are assumed to be non-interacting (\(\chi\!=\!0\)), entropic forces cause the two species to diffuse into one another. This process is characterized by an effective nonlinear diffusion equation.

To see this description, we utilize our 1D formalism, and add Eq. 7 8 . We define \(v\!\equiv\! u_A\phi+u_B(1-\phi)\), which can be interpreted as the mean velocity, and obtain \(\partial_x v\!=\!0\) in the absence of a source (\(S\!=\!0\)), implying \(v\!=\!0\) due to isotropy. We define the relative speed \(w\!\equiv\!u_A- u_B\), and using \(u_A\!=\!v+(1-\phi)w\), we obtain that in the absence of a source \(u_A\!=\!(1-\phi)w\). We then use Eq. 10 with Eq. 7 (neglecting the material stresses \(\sigma_A\) and \(\sigma_B\)), to write \[\label{eq:droplet951d} \partial_t \phi = \partial_x\left[\frac{\phi(1-\phi)^2}{\mu}\partial_x\Pi\right] \;,\tag{12}\] which corresponds to a nonlinear diffusion equation, with \[\label{eq:osmotic951d} \Pi\!=\!-\left[\log(1-\phi)+\left(1-N_A^{-1}\right)\phi\right] \;,\tag{13}\] where we set \(N_B\!=\!1\) and \(\chi\!=\!0\), and \[\label{eq:mu95exp} \mu = \phi\left(1-\phi\right)\frac{1-\phi+\alpha\phi}{1-\phi+\phi\left(1-\phi+\alpha\phi\right)}\tag{14}\] [see Eq. 4 ].

In the unentangled limit, we set \(\alpha\!=\!1\), reducing \(\mu_A\!=\!\phi\), and \(\mu_B\!=\!1-\phi\), and from Eq. 14 \(\mu\!=\!\phi(1-\phi)\) [see also Fig. 1 (d) for a plot of \(\mu(\phi)\) for different \(\alpha\)]. Using the 1D osmotic pressure of Eq. 13 , and \(N_A\!\gg\!1\), Eq. 12 becomes \[\label{eq:drop951d95unentangled} \partial_t \phi \!\simeq \!\partial_x\left[\frac{(1 + N_A \phi)}{N_A}\partial_x\phi\right] \underset{\phi\gg N_A^{-1}}{\rightarrow} \partial_x\left(\phi \partial_x \phi\right)\;.\tag{15}\] For sufficiently low volume fractions \(\phi N_A\!\ll\!1\), the above equation reduces to a linear diffusion equation. However, supposing that the initial volume fraction within the droplet is such that \(\phi\!\gg\!N_A^{-1}\), Eq. 15 approaches a nonlinear diffusion equation, as indicated above.

Using \(\ell\) as a typical length scale, Eq. 15 suggests \(\ell^2\!\sim\!\phi t\). We combine this scaling relation with fact that the total volume fraction \(\mathcal{I}\!\equiv\!\int_{-\ell}^{\ell}\phi dx\) must be conserved throughout the dynamics, implying \(\ell \phi\!\sim\! \mathcal{I}\). Put together, these scaling relations suggest \(\ell\!\sim\!\left(\mathcal{I}t\right)^{1/3}\), and \(\phi\!\sim\!\left(\mathcal{I}^2/t\right)^{1/3}\).

We now use the self-similarity ansatz, \(\phi(x,t)\!\equiv\!\left(\mathcal{I}^2/t\right)^{1/3}\Phi_u\left(\xi\right)\) (where \(\Phi_u\) denotes the unentangled \(\Phi\)), with the self-similar coordinate \(\xi\!\equiv\!x / \left(\mathcal{I}t\right)^{1/3}\). Substituting into the approximate nonlinear diffusion equation, Eq. 15 , yields an ordinary differential equation (ODE), \[\label{eq:drop951d95unentangled95ODE} \left(3\Phi_u\Phi_u'+\xi \Phi_u\right)'\!=\!0 \;.\tag{16}\] Integrating once, assuming a symmetric droplet profile \(\Phi_u'(0)\!=\!0\), we arrive at \(\Phi_u\!=\!-\frac{\xi^2}{6}+c\), where \(c\) is an integration constant. The obtained solution is a parabola that intersects the \(\xi\) axis \(\Phi_u(\xi_*)\!=\!0\), at \(\xi_*\!=\!\sqrt{6c}\). Integrating over the domain and demanding the volume fraction conservation \(\int_0^{\xi_*}\Phi_u\!=\!\frac{1}{2}\) gives \(c\!=\!\left(\frac{3}{32}\right)^{1/3}\).

We simulate the 1D system of Eqs. 7 11 , with a finite value of \(N_A\!=\!10^4\) and setting \(\alpha\!=\!1\) (using \(N_B\!=\!1\), \(\chi\!=\!0\) , and in the absence of \(\sigma_A\) and \(\sigma_B\)). We initialize the distribution of \(\phi\) as a Gaussian of typical length of \(2\), and normalize it such that at its peak, \(\phi\!=\!\phi_{\text{max}}\), i.e., \(\phi(x,t=0)\!=\!\phi_{\text{max}}\exp(-x^2/8)\), and follow its diffusion. We then solve Eqs. 7 11 , first by calculating the velocity profiles \(v\) and \(w\), followed by integrating \(\phi\) in time [54], [55].

We show in Fig. 2 (a) the \(\phi\) profiles obtained at different times during the diffusion dynamics. Then, we apply the self-similar transformation, and re-plot the solution in terms of the self-similar function \(\Phi_u\), and variable \(\xi\) in Fig. 2 (b). We also plot the theoretical solution \(\Phi_u\) on top. As time advances, the profiles approach the predicted solution, demonstrating the validity of the approximate nonlinear diffusion equation Eq. 15 for finite \(N_A\) values, and the emerging scaling relations above. In Fig. 2 (c) we plot trajectories of \(\phi\) profiles obtained for three different \(\mathcal{I}\) values (modified by changing \(\phi_{\text{max}}\)). Rescaling those profiles as shown in Fig. 2 (d), reveals the scalings and analytical solution applies to those different parameter values and at different times.

We repeat a similar analysis for the case of entangled A-mers. For \(\alpha\!\gg\!1\), we obtain \(\mu\!\rightarrow\!1-\phi\) [see Fig. 1 (d)], and Eq. 12 reduces to \[\label{eq:drop951d95entangled} \partial_t \phi \!\simeq \!\partial_x\left[\frac{(1 + N_A \phi)}{N_A}\phi \partial_x\phi\right] \underset{\phi\gg N_A^{-1}}{\rightarrow} \partial_x\left[ \phi^2 \partial_x \phi\right]\;.\tag{17}\] In this case, the scaling relations are \(\ell^2\!\sim\!\phi^2 t\) and \(\phi \ell \!=\! \mathcal{I}\), yielding \(\ell\!\sim\!\left(\mathcal{I}^2t\right)^{1/4}\), and \(\phi\!\sim\! \left(\mathcal{I}^2/t\right)^{1/4}\). Using the ansatz \(\phi(x,t)\!\equiv\!\left(\mathcal{I}^2/t\right)^{1/4} \Phi_e\left(\xi\right)\) (here with a subscript \(\Phi_e\) to denote the entangled case) with \(\xi\!\equiv\!x / \left(\mathcal{I}^2t\right)^{1/4}\), we reduce Eq. 17 to \[\label{eq:drop951d95entangled95ODE} \left(4\Phi_e^2\Phi_e'+\xi \Phi_e\right)'\!=\!0 \;.\tag{18}\]

Using a similar procedure as described above we obtain \(\Phi_e\!=\!\frac{1}{2}\sqrt{\frac{4}{\pi}-\xi^2}\). We show in Fig. 3 numerical solutions for a 1D system in the entangled case in similar fashion to Fig. 2 (with \(N_A\!=10^4\) and \(\alpha\!=\!10^4\)). Specifically, Fig. 3 (a) shows an exemplary trajectory at different times, and the inset shows a cusp occurring once \(\phi\!\simeq\!\alpha^{-1/2}\) — a “tail” where the simplification to \(\phi^2\partial_x\phi\) is no longer valid, and a different scaling law emerges. This tail signifies a transition between an entangled state to the unentangled state, characterized by a different diffusion scaling. Then, Fig. 3 (b) shows the collapse of the same trajectory using the suggested scalings, with the theoretical solution on top.

The agreement between the full numerical solution and the self-similar solution is partial — the droplet obtained from the numerical solutions of Eqs. 7 11 spreads further, and the ODE solution predicts higher central concentrations and lower concentrations at the rims. We suspect this partial agreement is due to the tails referred to above. As these tails move material away from the central droplet, we anticipate that as their relative contribution decreases, the solutions should approach the self-similar solution of Eq. 18 . In Fig. 3 (c) we show solutions with different \(\mathcal{I}\) (obtained by changing the \(\phi_{\text{max}}\) values), and in Fig. 3 (d) we show their rescaled counterparts. As the initial maximal concentration \(\phi_{\text{max}}\) is increased, \(\mathcal{I}\) increases, and the relative contribution of the tail to the dynamics decreases. The rescaled solutions indeed confirm that as the tail contributions become negligible (higher \(\phi_{\text{max}}\) and \(\mathcal{I}\) values) the solutions approach the expected self-similar form.

In two dimensions (2D) we can repeat the procedure above, to approximate the problem as a nonlinear diffusion equation, which could be cast in a self-similar framework (see App. 6). Using \(\boldsymbol{r}\) as the position vector, and assuming isotropy, we use the radial distance \(r\!\equiv\!|\boldsymbol{r}|\) to get \(\phi(r,t)\!\equiv\!\left(\mathcal{I}/t\right)^{1/2} \Phi_u(\xi)\), and \(\xi\!\equiv\!r / \left(\mathcal{I}t\right)^{1/4}\), and \(\phi (r,t)\!\equiv\!\left(\mathcal{I}/t\right)^{1/3} \Phi_e(\xi)\), and \(\xi\!\equiv\!r / \left(\mathcal{I}^2 t\right)^{1/6}\) for the unentangled and entangled cases, respectively. Using the transformations we obtain \[\begin{gather} \xi^2 \Phi_u + 4 \xi \Phi_u \Phi_u' = 0 \;\quad \text{for} \;\alpha\!=\! 1 \;, \tag{19} \\ \xi^2 \Phi_e + 6 \xi \Phi_e^2 \Phi_e' = 0 \;\quad \text{for} \;\alpha\!\gg\! 1 \;.\tag{20} \end{gather}\] Solving those in a similar fashion to the 1D case yields \(\Phi_u\!=\!\frac{1}{2}\left(\frac{1}{\sqrt{\pi}} - \frac{\xi^2}{4}\right)\), and \(\Phi_e\!=\!\frac{1}{2}\sqrt{\left(\frac{2}{\pi}\right)^{2/3} - \frac{2\xi^2}{3}}\). We show our numerical solutions to Eqs. 1 3 (neglecting \(\boldsymbol{\sigma}_A\) and \(\boldsymbol{\sigma}_B\)), and the comparison to the expected self-similar solutions in Fig. 4.

Analysing the passive diffusion problems, we note the similarity in the attained spatial form between the 1D and 2D cases, i.e., between Figs. 2, 3, and  4. Specifically, the unentangled cases vary with \(\xi^2\), and the entangled cases vary as \(\sqrt{\xi_*^2-\xi^2}\). Those similar spatial profiles hint that dimensionality plays a rather minor role in terms of the observed spatial profiles.

Generically, the suggested simplifications reduce both the 1D and 2D cases to a nonlinear diffusion problem with \(\partial_t \phi \!=\!\nabla\cdot \left(\phi^{m-1}\nabla \phi\right)\), where \(m\) quantifies the degree of nonlinearity that directly depends on \(\alpha\) — the unentangled case implies \(\alpha\!=\!1\) leading to \(m\!=\!2\), while the entangled case has \(\alpha\!\gg\!1\) leading to \(m\!=\!3\). Interestingly, similar equations appear in the context of diffusion in porous media [56], [57]. The diffusive self-similar length-scale and amplitudes depend on the nonlinearity \(m\), and the dimensionality \(d\) of the problem. A passive droplet diffusion problem in \(d\) dimensions, with a nonlinearity of power \(m\), suggests the scaling \(\phi\!\sim\!\mathcal{I}^{\gamma_1}t^{\delta_1}\) and \(\ell\!\sim\!\mathcal{I}^{\gamma_2}t^{\delta_2}\), with \[\label{eq:scaling95diffusion} \begin{align} &\gamma_1\!=\!2/p \;, & \, \delta_1\!=\!-d/p \;, & \\ &\gamma_2\!=\!(m-1)/p \;,& \, \delta_2\!=\!1/p \;,& \end{align}\tag{21}\] with \(p\!\equiv\!2+d(m-1)\). Transitioning from the unentangled to entangled scenario slows the dynamics down — \(m\) increases causing \(p\) to increase as well. This parallels with polymer reptation in entangled scenarios, absent from the unentangled case [9], [10], [16]. So, while reptation is a single polymer mode-of-motion, it affects the emerging diffusion on a coarse-grained scale. In both cases, approximating the equations of motion in the limit of sufficiently large \(N_A\) and \(\phi\), and using the corresponding \(\alpha\) values (while neglecting the material stresses \(\sigma_A\) and \(\sigma_B\)) allowed us to reformulate the interdiffusion problem as a self-similar droplet diffusion problem. The emerging equations were simple enough to solve analytically, and resulted in good agreement with both our 1D and 2D numerical solutions.

4 Source-driven diffusion↩︎

The passive droplet diffusion scenario investigated above may become modified in certain biological cases. Some biological condensates serve as effective sources, e.g., in the synthesis of RNAs and proteins, injecting polymers at a certain rate into the system [36]–[40]. We now modify the formalism above to account for these scenarios, treating the synthesis centers as external polymer sources, depicted schematically in Fig. 1 (b).

We employ a similar numerical solution method as for the passive case discussed in Sec. 3 to obtain numerical solutions for the source-driven diffusion. Here, we set \(\phi(x,t=0)\!=\!0\), and use a localized Gaussian as our source term, \(S(x,t)\!=\!\frac{f(t)}{Z} \exp(-x^2/2)\), with the normalization factor \(Z\!\equiv\!\sum\exp(-x^2/2)\;\Delta x\) over our integration domain (here \(\Delta x\) is the resolution of the discretized domain). We show two exemplary trajectories for the unentangled \(\alpha\!=\!1\) case in Fig. 5 (a), and \(\alpha\!\gg\!1\) in Fig. 5 (c), for a constant-rate, localized source \(f(t)\!=\!Q\), as obtained from our numerical simulations of Eqs. 7 11 . The A-mers accumulate near the origin, increasing \(\phi\) locally, hence increasing the osmotic pressure/chemical gradients, driving outward diffusion. The volume fraction profile evolves from its early time behavior, until it converges to the long time profile. While there are intricate spatial features at early times, and near the origin, e.g., the presence of a “tip”, in both cases the volume fraction distributions near the edges, or droplet’s “front shape", attained at long times are reminiscent of their spatial features in the passive droplet expansion (see Figs. 2-3).

To characterize the volume fraction redistribution dynamics, we measure two typical droplet radii, denoted by \(\ell_{in}\) and \(\ell_{out}\), defined here, respectively, as the half-width of the \(\phi\!>\!0.1\) domain for \(\ell_{in}\), and \(\phi\!>\!0.01\) domain for \(\ell_{out}\) [as annotated in Fig. 5 (a) and (c)]. We plot their time-dependencies for various source rates \(Q\), as shown in Fig. 5 (b) and (d) for the unentangled and entangled cases, respectively. At very early times \(\ell_{in}\) approximately increases according to a linear diffusion equation response \(\sim\!t^{1/2}\), while \(\ell_{out}\) diffuses according to the passive droplet cases analyzed in Sec. 3, i.e., following Eq. 21 . At longer times, the two radii converge into similar nonlinear scalings. These early-time, spatially-dependent scalings indicate the volume fraction redistributes heterogeneously throughout the droplet.

To try and understand these spatial and temporal behaviors, we revisit the continuity equations Eq. 7 8 . Upon addition, and using \(v\!\equiv\!\phi u_A + (1-\phi) u_B\), we obtain \(\partial_x v\!=\!S\). We then use this together with \(w\!=\!u_A-u_B\), to express \(u_A\), and substitute back into Eq. 7 , leading to the evolution equation \[\label{eq:source951d} \partial_t \phi = \partial_x\left[\frac{\phi\left(1-\phi\right)^2}{\mu}\partial_x \Pi - v \phi \right] + S \;,\tag{22}\] Compared to the passive droplet equivalent, Eq. 12 , the additional source term not only contributes by adding more A-mers into the system, but also modifies the average velocity \(v\) and contributes a flux term. We expect a similar effect happens in higher dimensions as well, as a consequence of a material source near the origin (see App. 7).

We now suppose that the source is spatially-localized, and takes the general form \(S(x,t)\!=\!\delta(x) f(t)\). Combining the unentangled and entangled cases, in the limit of \(\phi\!\gg\!N_A^{-1}\), Eq. 22 reduces to \[\label{eq:source951d95simplified} \partial_t \phi + v \partial_x \phi \!\simeq\! \partial_x\left(\phi^{m-1}\partial_x \phi \right) + \delta(x)f(t)\left(1-\phi\right) \;,\tag{23}\] where we have used the fact that \(\partial_x v\!\equiv\!\delta(x)f(t)\). Irrespective of the nonlinearity exponent \(m\), the presence of the advection term \(v\partial_x \phi\) suggests transforming the above into a reduced ODE [e.g., Eqs. 16 , 18 ,  19 and  20 ] may involve translation in addition to rescaling. Additionally, combining the source together with the incompressibility \(\phi_A + \phi_B\!=\!1\) renders the source’s amplitude to depend on the local concentration \(\phi\), as \(\delta(x)f(t)\left(1-\phi\right)\).

With these features in mind, we consider again the dominant balance of the osmotic pressure with the concentration time-variation \(\ell^2 \!\sim\!\phi^{m-1} t\), together with the increase in mass \(\phi \ell \!\sim\!\int_0^t f(t')dt'\). To make progress, we assume \(f\) has a temporal power-law variation \(f\!=\!Q t^{q-1}\), allowing us to obtain the scaling relations for any power-law increase \(q\), and to generalize for more complicated time-variations (potentially approximating those as piecewise power-law functions). Taking \(Q t^{q-1}\) for our source, we obtain \(\phi \ell \!\sim\! Q t^q\), from which we recover \(\phi(x,t)\!\equiv\!Q^{\tilde{\gamma}_1}t^{\tilde{\delta}_1} \Phi(\xi)\) and \(\xi\!\equiv\!x/Q^{\tilde{\gamma}_2}t^{\tilde{\delta}_2}\) with \(p\!\equiv\!2+d(m-1)\) as defined above, and \[\label{eq:scaling95source} \begin{align} &\tilde{\gamma}_1\!=\!2/p \;,& \, \tilde{\delta}_1\!=\!(2q-d)/p \;, &\\ & \tilde{\gamma}_2\!=\!(m-1)/p \;,& \, \tilde{\delta}_2\!=\!\left[1+q(m-1)\right]/p \;.& \end{align}\tag{24}\]

Comparing Eq. 24 with the scaling relations in the passive droplet diffusion case, Eq. 21 , reveals that only the time dependence, via the \(\tilde{\delta}\) exponents, are modified by the presence of a source, in which case \(Q\) plays a role reminiscent of \(\mathcal{I}\) in the passive case. These scaling exponents are the ones shown in the late-time regimes in Fig. 5 (b) and (d), signifying the droplet’s expansion before saturation.

Combining the scaling and taking into account the advection term, we propose to use the 1D transformation \(\phi\!=\!Q^{\tilde{\gamma}_1}t^{\tilde{\delta}_1} \Phi(\eta)\) and \(\eta\!=\!\left[x - \frac{1}{2 q}Q t^{q}\right]/ Q^{\tilde{\gamma}_2}t^{\tilde{\delta}_2}\!\equiv\!\xi+\xi_s\), where we defined \(\xi_s\!\equiv\!-\frac{Q^{1-\tilde{\gamma}_2} t^{q-\tilde{\delta}_2}}{2q}\) as the rescaled source coordinates (see App. 8). Using these in Eq. 23 (setting aside the source term, which we treat as a boundary condition to our equations), we obtain \[\begin{gather} \left[3\Phi_u\Phi_u'+(1-2q)\eta\Phi_u\right]'+3 q \eta \Phi_u' = 0 \quad \text{for} \;\alpha\!=\!1 \;, \tag{25} \\ \left[4\Phi_e^2\Phi_e'+(1-2q)\eta\Phi_e\right]' + 4q\eta\Phi_e' = 0 \quad \text{for} \;\alpha\!\gg\!1 \;. \tag{26} \end{gather}\] Clearly, in both cases the emerging equations are cast as ordinary differential equations, however in the self-similar coordinates \(\eta\) the source is moving back to \(\xi_s\), and the boundary conditions themselves are time-dependent. So, while the equations are ordinary differential equations, the domain of integration is time dependent. This explicit time dependence breaks the self-similar structure.

Instead of solving the ODEs of Eq. 25 26 with a time-dependent domain and boundary conditions, we perform an asymptotic expansion near the front \(\eta_*\), where \(\Phi(\eta_*)\!=\!0\). We anticipate that while the solutions’ features near the source may vary, the droplets’ fronts will attain a similar diffusion-dominated, spatial dependence as their passive counterparts, as already demonstrated in Fig. 5 (a) and (c). Substituting the expansion \(\tilde{\Phi}\simeq a_1 y^{1/(m-1)} + a_2 y^{2/(m-1)}\) (where \(y\!\equiv\!\eta_*-\eta\)) in Eq. 25 26 , using \(m\!=\!2\) and \(m\!=\!3\) for \(\Phi_u\) and \(\Phi_e\) respectively, collecting the leading orders, and demanding that those vanish, lead to \(a_1\!=\!\frac{2}{3}\eta_*\), and \(a_1\!=\!\sqrt{\frac{3}{2}\eta_*}\) for \(m\!=\!2\) and \(m\!=\!3\) respectively.

We rescale the numerical solutions obtained for the unentangled and entangled cases, at different times and for different amplitudes \(Q\), according to Eq. 24 , as shown in Fig. 6 (a) and (b). For high \(Q\) values, our numerical integration did not advance far enough to transition into the long-time scaling regime, but the other solutions converge onto similar curves. The obtained curves vary especially near the origin, but as the front is approached, the different solutions converge onto a similar profile. In the insets in Fig. 6 (a) and (b), we show the rescaled solutions in \(\eta\) coordinates, shifted to the front \(\eta_*\), and plot the leading order expected behavior discussed above. Clearly, the well-developed fronts all agree with our predicted leading-order behavior.

The constant source with \(q\!=\!1\) is a special case in which the source \(\xi_s\) moves with the same time-dependence as the rescaled coordinates \(\xi\). We now show that similar spatial profiles arise for \(q\!=\!2\), i.e., with a time-dependent source. We plot the rescaled numerical solutions at various times and for different \(Q\) values in Fig. 6 (c) and (d), showing an explicit linear front geometry for the unentangled case in panel (c). In panel (d) we include an inset plotting \(\Phi_e^2\) versus \(\eta-\eta_*\), showing that \(\Phi_e\!\sim\!\sqrt{\eta_*-\eta}\), similar to the \(q\!=\!1\) case above.

Finally, we perform 2D numerical simulations with \(q\!=\!1\), and compare our results in a similar fashion in Fig. 7. Panels (a) and (b) show exemplary radially-averaged \(\phi\) dynamics. In 2D, taking \(q\!=\!1\) reveals an interesting case in which the source amplitude is time-independent [see \(\tilde{\delta}_1\) in Eq. 24 ], and \(\tilde{\delta}_2\!=\!1/2\) independently of \(m\), implying the unentangled and entangled time-dependencies are similar. We show this scaling in the insets of panels (a) and (b) respectively (here probing only a single length-scale \(\ell_{in}\), as the half-width of the \(\phi\!>\!0.01\) domain). Then, panels (c) and (d) (and the inset therein) show the rescaled solutions, revealing similar spatial features at the droplets’ front, as the 1D counterparts (see Fig. 6).

5 Discussion↩︎

Inspired by biocondensate diffusive and mixing dynamics, we used a two-fluid formalism to analyze polymer-polymer interdiffusion problems. We presented the equations of motion and demonstrated that under certain approximations the interdiffusion of a polymer droplet could be cast as a nonlinear diffusion equation, similar to diffusion in porous media [56], [57]. We highlighted the fact that the entanglement parameter \(\alpha\) can dramatically alter the nonlinearity exponent \(m\), changing the anticipated scaling and spatial structure of the solutions. The results obtained from our full numerical simulations in both 1D and 2D corresponded well with the anticipated self-similar analytical solutions for the two \(\alpha\) limits considered for unentangled case \(m\!=\!2\) (corresponding to \(\alpha\!=\!1\)), and the entangled case \(m\!=\!3\) (corresponding to \(\alpha\!\gg\!1\)).

Then, we introduced a source into the equations of motion — akin to a continuous synthesis of polymers, prevalent in biological contexts. We demonstrated that for this case there are multiple temporal scaling laws that characterize the spread at different time regimes, and that while the source alters the time-dependence of the solutions and breaks the self-similar structure, the local geometry near the front bears similarity to the “passive diffusion” cases. Using both 1D and 2D numerical solutions, we demonstrated this similarity persists for various \(q\) values, and in different spatial dimensions \(d\).

We suspect that the local front geometry could affect polymer conformations [58]. The sharp interface of the droplets hints that the velocity profile should vary dramatically there, and should converge to a constant sufficiently far away. This local variation at the rim induces strain-rate gradients, and induce compressive forces on the polymers, causing them to adopt folded conformations. We suspect such effect could partially explain the compression of ribosomal RNA conformations as they drift away from the nucleolus [38], [58].

In the above work we neglected elastic and viscous effects that could modify the shape of the formed droplets, and cause the solutions to deviate from the reduced self-similar solutions. Such coupling between osmotic pressure and mechanics lies at the heart of gel physics [24], [26], [29], [31], [47]–[50]. While the source-driven polymer interdiffusion problem considered above is reminiscent of a driven Cahn–Hilliard system [59], [60], the presence of material stresses may introduce additional temporal dependencies and spatial patterns, especially near the center and edges of the expanding droplets.

It may be interesting to consider a source-driven problem for phase separating mixtures. In such cases, the osmotic pressure would oppose the outward flux, which would eventually lead to droplet formation. We suspect this scenario could also be relevant for biological phase separation [61], [62]. Overall, we believe the presented two-fluid framework, the analysis above, and future works would deepen our understanding of polymer interdiffusion physics and biophysical diffusion processes.

Acknowledgments. — A.M. acknowledges support from the James S. McDonnell Foundation Postdoctoral Fellowship Award in Complex Systems. We acknowledge support from the NSF Grant DMS/NIGMS 2245850 and from the Princeton Center for Complex Materials (PCCM), an NSF-supported Materials Research Science and Engineering Center under award DMR-2011750. The authors also thank N. Wingreen and C. P. Brangwynne for helpful comments and discussions.

6 Two-dimensional passive droplet diffusion↩︎

In 2D we assume that the flow is radial, implying the radial velocities \(u^r_A\) and \(u^r_B\) obey \(\nabla\cdot\left(\phi u_A^{r} + \left(1-\phi\right)u_B^r\right)\equiv\!\nabla\cdot \boldsymbol{v}\!=\!0\), implying \(v^r\sim \frac{v_0}{r}\). From here we use similar simplifications as those used for the 1D case, neglecting the material stresses \(\boldsymbol{\sigma}_A\) and \(\boldsymbol{\sigma}_B\), so that \(w^r\equiv u_A^r - u_B^r\!\simeq\!-\frac{\left(1-\phi\right )}{\mu}\partial_r \Pi\) (\(\boldsymbol{\Pi}\!\equiv\! \boldsymbol{I} \Pi\), assumed to be isotropic). Expressing \(u_A^r\!=\!v^r + \frac{(1-\phi)^2}{\mu}\partial_r \Pi\), and using this back in Eq. 1 (without a source) yields \[\partial_t \phi\!\simeq\! \frac{1}{r}\partial_r \left(r\phi^{m-1}\partial_r \phi\right)-v^r\partial_r\phi \;,\] where we used \(\nabla\cdot \boldsymbol{v}\!=\!0\). While \(v^r\simeq \frac{v_0}{r}\) could contribute significantly near the origin, as we are examining a droplet expansion, we assume that at sufficiently late times when the droplet has extended, this term can be neglected, yielding a similar nonlinear diffusion equation as in the 1D case.

From here we proceed with the usual scaling arguments. The nonlinear equation implies \(\ell^2\sim \phi^{m-1} t\), and the integral conservation law implies \(\phi \ell^2\!\sim\mathcal{I}\), yielding \(\phi\sim\left(\mathcal{I}/t\right)^{1/m}\), and \(\ell\sim \mathcal{I}^{\frac{m-1}{2m}}t^{\frac{1}{2m}}\), as obtained from Eq. 21 (with \(d\!=\!2\)).

7 Higher-dimensional problems with a source↩︎

In higher spatial dimensions, we add Eq. 1 to obtain \[\label{eq:incomp95higherD} \nabla\cdot\boldsymbol{v}\!=\!S \;.\tag{27}\] Now, using \(\boldsymbol{w}\) and \(\boldsymbol{v}\) to express \(\boldsymbol{u}_A\!=\!\boldsymbol{v}+(1-\phi)\boldsymbol{w}\), we can rewrite \[\partial_t \phi = \nabla \cdot \left[\frac{(1-\phi)^2\phi}{\mu}\nabla\cdot\boldsymbol{\Pi} - \boldsymbol{v} \phi \right] + S \;,\] analogous to Eq. 22 , where here we again neglect the material stresses \(\boldsymbol{\sigma}_A\) and \(\boldsymbol{\sigma}_B\). We can then express this also as \[\partial_t \phi + \boldsymbol{v}\cdot\nabla\phi = \nabla \cdot \left[\frac{(1-\phi)^2\phi}{\mu}\nabla\cdot\boldsymbol{\Pi}\right] + S(1-\phi) \;,\] analogous to Eq. 23 , where we used Eq. 27 . This explicitly shows that the addition of a source implies an advection term, and a modification to the source amplitude, as was shown in 1D above.

8 Transformation into self-similar variables for the source-driven problem↩︎

Here we provide detailed calculations that show how using the transformations \(\phi\!=\!Q^{\tilde{\gamma}_1}t^{\tilde{\delta}_1} \Phi(\eta)\) and \(\eta\!=\!\left[x - \frac{1}{2 q}Q t^{q}\right]/ Q^{\tilde{\gamma}_2}t^{\tilde{\delta}_2}\) proposed above in Eq. 22 result in ordinary differential equation.

To begin we evaluate \[\begin{gather} \partial_t \eta = -\frac{1}{2} Q^{1-\tilde{\gamma}_2} t^{q-1-\tilde{\delta}_2} - \tilde{\delta}_2 t^{-1} \eta\;, \\ \partial_x \eta = Q^{-\tilde{\gamma}_2}t^{-\tilde{\delta}_2} \;. \end{gather}\] Taking similar derivatives for \(\phi\), we have \[\begin{gather} \partial_t \phi = \tilde{\delta}_1 Q^{\tilde{\gamma}_1}t^{\tilde{\delta}_1-1} \Phi + Q^{\tilde{\gamma}_1}t^{\tilde{\delta}_1} \Phi' \partial_t \eta \;, \\ \partial_x \phi =Q^{\tilde{\gamma}_1}t^{\tilde{\delta}_1} \Phi' \partial_x \eta . \end{gather}\]

Using the above relations in the left-hand side of Eq. 22 (focusing on \(x\!>\!0\) for simplicity), results in \[\partial_t \phi + \frac{1}{2} Qt^{q-1}\partial_x \phi = Q^{\tilde{\gamma}_1}t^{\tilde{\delta}_1-1}\left(\tilde{\delta}_1\Phi - \tilde{\delta}_2\eta \Phi'\right) \;.\] Similarly, using the above in the nonlinear diffusion term yields \[\partial_x\left(\phi^{m-1}\partial_x \phi\right) = Q^{m\tilde{\gamma}_{1}-2\tilde{\gamma}_{2}}t^{m\tilde{\delta}_{1}-2\tilde{\delta}_{2}}\left[\left(m-1\right)\Phi'^2 + \Phi \Phi''\right]\Phi^{m-2} \;.\] Comparing the powers of these two expression reveals \[\begin{gather} \left(m-1\right)\tilde{\gamma}_1 \!=\! 2\tilde{\gamma}_{2}\;, \\ \left(m-1\right)\tilde{\delta}_1 \!=\!2\tilde{\delta}_{2}-1 \;, \end{gather}\] both obeyed by Eq. 24 . Using \(m\!=\!2\) and \(m\!=\!3\) yield Eqs. 25 26 respectively.

References↩︎