{'fraction_non_alphanumeric': 0.032259344496468535, 'fraction_numerical': 0.010634076660582494, 'mean_word_length': 4.737081720919645, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Coalescing Clusters Unveil New Regimes of Frictional Fluid Mechanics

Droplet coalescence is essential in a host of biological and industrial processes, involving complex systems as diverse as cellular aggregates, colloidal suspensions, and polymeric liquids. Classical solutions for the time evolution of coalescing clusters are typically based on tractable limiting physics, such as analytical solutions to the Stokes equation. By combining computational and theoretical analyses, we show that there is an unexplored family of coalescence processes: those governed by highly dissipative coupling to the environment. This leads to new scaling laws characterizing droplet coalescence, as well as new time-invariant parameterizations of the shape evolution of the coalescing system. We demonstrate these effects via particle-based simulations and both continuum and boundary-integral solutions to hydrodynamic equations, which we then understand in the context of a generalized Navier-Stokes-like equation. Our theoretical description of highly frictional coalescence mathematically maps onto Darcy flow in the presence of surface tension effects, opening up exciting avenues of research in applying well-studied fluid dynamical techniques to a broad range of novel systems.

as the theoretical basis for the study of fusion in such liquid-like systems as aggregates of living cells [7][8][9], shedding light on critical biological processes in embryonic development, cancer invasion, and tissue engineering. In the case of simple liquids at low Reynolds number this phenomenon is well-understood, as conformal mapping techniques can be used to find analytical solutions to the relevant Stokes equation [10,11]. The coalescence of more complex liquid or liquid-like systems -ranging from polymers to colloidal suspensions to aggregates of living cells -lacks such a general understanding. In such settings one often uses the analytical simple-liquid theory as a starting point, layered with additional complexities to either explain deviations from the expected behavior (as in the case of viscoelastic effects that result in arrested coalescence of colloidal droplets with internal microstructure [12] or new scaling laws observed in polymer droplet dynamics polymers [13]) or as a tool to infer the effective surface tension in an otherwise difficult-to-measure system (as in the case of observing both coalescence dynamics and surface tension fluctuations in the fusion of nucleoli [14]). In many of these more complex systems -such as colloids in a viscous solvent or nucleoli in a crowded intranuclear environment -we see that the coalescing components are strongly coupled to a viscous, dissipative background. We argue that such systems are better understood at the hydrodynamical level by adding an additional frictional term to a Navier-Stokes-like description [15], just as such systems are often modeled at the microscopic level by considering Langevin (or Brownian) Dynamics rather than Molecular Dynamics [16]. Such a continuum model is initimately linked to Darcy's Law (a mathematical description of fluid flow through porous media), but here in the context of a compact systems governed by surface tension. We show that this fundamentally changes the nature of droplet coalescence, leading both to a new evolution of the overall shape of the coalescing droplet and new scaling laws characterizing the growth of the neck.
To qualitatively demonstrate how two circular droplets coalesce differently as a function of the microscopic dynamics to which their constituents are subject, we first conduct agent-based simulations in which attractive, Lennard-Jones particles in two coalescing droplets evolve in time according to either standard Molecular Dynamics (MD) or Brownian dynamics (BD) -see Methods section. In Fig. 1(a-d) we show both representative snapshots of this process and ensemble-averaged contours representing the boundary of the coalescing droplets. Immediately apparent is how, in the case of MD, as the "neck" of the droplet grows the overall length of the coalescing droplets decreases. This evolution of shapes is in close quantitative agreement with the contour evolution predicted by directly solving the Stokes equation (shown in Fig. 1(e)) , which was done for this geometry in the classic work of Hopper [10,11]. In contrast, in our BD simulations we see that the overall droplet length changes very little even while the neck grows substantially. Strikingly, a similar qualitative trend was recently observed in nucleoli fusion [14]. In Fig. 1(f) we show that this type of shape evolution appears similar to the result of solving a very different hydrodynamic description -that of Darcy's Law, ∇P = −u. This connection is not accidental: in both BD and fluid flow through pores the momentum of the flow is quickly dissipated to its surroundings and indeed, we will later show that Darcy's Law is an extreme, overdamped limit of a more general Navier-Stokes-like description.
To quantify these differences in shape we focus on two measures of the coalescing droplet geometry: the overall length L and the neck radius R n , as schematically defined in Fig. 2(a) and which we measure in units of the radius of the original drops R. Both of these, and especially the scaling of the neck radius, are common experimental observables used to measure coalescence and extract macroscopic parameters (such as the fluid's viscosity or the effective surface tension governing the coalescence), and as such are natural targets of theoretical analysis. We study the detailed scaling of the neck radius below, and first focus on the "shape evolution curve" of the coalescing droplets in Fig. 2(b). This is a time-independent parameterization of the overall coalescence process, and which we define as the decrease in combined droplet length as a function of the neck radius. In this representation, we find our MD simulations are indeed in quantitative agreement with Hopper's result. We also show the shape evolution curve predicted by an earlier, approximate theory due to Frenkel [17], which makes strong assumptions about the coalescing geometry and which is still frequently used to compare with experimental data on complex systems [7,8,18]. Again contrasting with this result, the BD simulations show strong deviations from Hopper's prediction, and are characterized by an initial plateau in the shape evolution curve. To help compare with experiments, we extract the data on nucleoli fusion from Fig. 1d of Ref. [14]. This data also shows an initial plateau, albeit less pronounced than that of our BD simulations (discussed further below). Another simple yet powerful way to quantify the coalescence of two droplets is through the time evolution of the neck radius. For simple liquids this can be understood via the analytical study of the Stokes equation by Hopper [10], by direct numerical simulations of the Navier-Stokes equation [19,20], and by qualitative theoretical approximations that use force balance at the droplet interface [21] -all of which are well-supported by experimental data [22][23][24][25][26] and MD simulations [27,28]. In simple liquid systems the scaling of the neck growth is determined by a viscous regime at early times, in which R viscous n ∼ t or R viscous n ∼ −t log t, depending on the viscosity and geometrical details of the singular contact point [20,24]. These dynamics can give way to an inertiallydominated regime at late times, in which R inertial n ∼ t 1/2 . However, our MD simulations represent extreme small, viscous drops, where momentum diffuses extremely fast compared to the motion of the neck. Because of this we do not observe the crossover to the inertial regime ( Fig. 2(c)), and the dynamics closely follow just the purely viscous coalescence solution of Hopper [10,11].
In contrast, BD shows a much slower growth of the neck radius with time. It also appears to have two dynamical regimes, but with a scaling law (over a limited range) close to t 1/2 at early times crossing over to a scaling law obviously smaller than t 1/2 (approximately t 1/4 , although the measurement is limited to a relatively small dynamical range) at later times. This has no analog in the analyses of the Navier-Stokes equation discussed above. Suggestively, in Fig. 2(c) we also show the solution of Darcy's Law for this geometry, in which the effect of the surface tension is balanced only by friction −ζu; it also displays this smaller scaling law. These initial conditions for these simulations are extremely singular, and as a result the earliest time data is affected by the detailed assumptions one makes for the initial contacting geometry. We characterize the neck width scaling in Darcy's Law as R n ∼ t α , and find α 1/4. In our BD simulations we similarly label the early-time apparent power-law scaling as t β with β ≈ 1/2 to make it clear that this is a new regime with no a priori relation to the inertial regime in simple liquids (despite the similar exponents).
The stress caused by surface tension is the driving force behind coalesence, and it is set by γκ, the product of the surface tension and the boundary curvature. In Fig. 2(d) we plot the maximum curvature κ (located at the neck) versus R n . We find, with no adjustable parameters, that the curvatures from Darcy's Law and from our BD simulations agree almost quantitatively after a characteristic point. The exception to this agreement occurs at early times. At early times the curvature in our Darcy's Law simulation is strongly affected by our choice of initial conditions (which we, in fact, take from analytical solutions to Stokes equation at R n ≈ 0.2; this solution has a larger curvature than that of the actual Darcy's Law solution). The curvature in our BD simulations is limited by the finite size r e of the particles themselves (i.e., the maximum curvature in our particulate simulations is bounded by ∼ 1/r e ). As such, we expect that the crossover point between the two regimes shown in Fig. 2(c) for Brownian Dynamics simulations should be roughly the crossover point in Fig. 2(d) (both labeled with the green star). It is also natural to predict that in the particle-based simulations, a smaller system size will have a greater finite size effect and thus a later crossover point, as shown in Fig. 2(e) (see Methods for more details of estimating this crossover). After this crossover there is a regime in which κ ∼ R −3 n for the Stokes, Darcy's Law, and BD solutions (the origin of the late-time neck width scaling laws). Before this crossover the curvatures for the Darcy's law and BD simulations deviate from the R −3 n power law in opposite directions, resulting in corrections to the scaling laws of neck growth in the early stage that deviating from t α scaling law in opposite directions. Although Darcy's Law appears to be the correct continuum description for our BD simulations, real systems have neither completely vanishing inertia nor completely vanishing viscosity. To account for behaviors which interpolate between these extremes -as, for instance, is seen in the short but distinct plateau in the shape evolution of the coalescing nucleoli -we consider the more general case of Langevin Dynamics (LD) (see Methods). The coupled Langevin equations describing the movement of the damped particles in the discrete systems can be coarse-grained to a continuum form using the Mori-Zwanzig technique [29,30]. Considering only the long-wavelength limit and neglecting the hydrodynamic effect of the solvent, this continuum description becomes Here, the µ∇ 2 u term represents the viscosity originating from the direct interactions between the colloidal particles or biological cells, and −ζu represents the additional dissipative effects of the environment (e.g., the frequent collisions of the solvent molecules on colloidal particles, or the frictional effects on cellular systems from substrates or the extracellular matrix).
This non-dimensionalization trivially gives an overall system-size dependence for the characteristic time of droplet coalescence, which is supported by our discrete simulations shown in Extended Fig. 2.
As emphasized above, we expect that these extreme limits of Eq. 1 may not appropriately capture all the relevant physics of a real system. Starting with the non-dimensionalization in the viscous case, we nondimensionalize the more general equation as Here we have used the Laplace number (La) as the natural coefficient of the inertial term and the inverse Darcy number (Da −1 ) as the natural coefficient of the frictional term (drawing an analogy between ζ in our case and the permeability K in the usual formulation of Darcy's Law via ζ = µ K ). We note that La = Oh −2 , where Oh is the Ohnesorge number typically used to quantify coalescence dynamics [24]. Even in simple liquids with ζ = 0, we see that a larger value of the Laplace number will lead to a larger deviation from Hopper's solution to the Stokes equation. We test this numerically by solving the standard Navier-Stokes equation in our problem geometry with different values of µ, as shown in Fig. 3(a) (see Methods for details). This result is supported by the recent MD simulations of a model with low viscosity in Ref. [31]; in that work there is a clear departure from the shape evolution curve predicted by Hopper, with a growing shape-evolution-curve plateau as their simulated system size increases. This is precisely what one would expect from the above analysis.
In the case of BD, we see that a larger value of the Darcy number (which corresponds to a larger role for the viscosity in Eq. 4 ) should lead to results which increasingly deviate from the solution based on Darcy's Law. This explains the system-size effect seen in the shape evolution curves of Fig. 3(b), in which the smallest simulated systems deviate strongly from the predictions of the continuum theory (see also Extended Fig. 2b). These analyses, together with numerical simulations of LD using various scales of the viscous damping coefficient λ seen in Fig. 3(c), are summarized into a qualitative "phase" diagram of droplet coalescent motifs in Fig. 3(d).
We briefly comment on a crucial underlying physical reason for the differences observed in the shape evolution of droplets under different microscopic dynamics. We see that there are deviations from Hopper's classical prediction, qualitatively evident in the form of a plateau in the shape evolution curves, for both frictionally dominated systems and inviscid simple liquids. Are the explanations for deviations in these quite physically different systems related? We propose that a simple way to think about this feature of droplet coalescence is in terms of the propagation of perturbations. In particular, the time at which the system length begins decreasing is fundamentally limited by the time needed for information about the formation of an initial contact point between the two disks to propagate to the far boundary of the droplets. One can theoretically quantify this time of propagation via the dynamical structure factor for the complex liquid, S(k, z), which measures particle correlations in both space and time. As encoded by the complex poles of S(k, z) (see Ref. [29]), in momentumconserving systems like simple liquids perturbations are propagated via soundwave modes. The Stokes equation is the extreme case when the sound wave speed is infinite so that all the other cases with finite wave speed deviate from it in the direction of elongating the plateau length. This effect is usually hard to capture as the wave speed is large compared to the drop size, but is more obvious in systems with small viscosity because the neck grows more within a given time window due to finite wave speed. This is seen in both experiments [24] and simulations [31]. Weak frictional or other dissipative effects lead to an exponentially decaying coefficient for the sound wave, but for sufficiently strong dissipative effects (as in overdamped systems) perturbations propagate in a completely diffusive mode with x 2 ∝ Dt. This qualitative difference in information propagation is illustrated in Fig. 4, in which we color the particles near the neck region at the beginning of the simulation with red and then take a snapshot at the time point when the length starts decreasing (quantified by L ≈ 0.9L 0 ) in MD and BD separately. Just as the Stokes equation is the limiting case of infinite sound speed, Darcy's Law is the limiting case of infinite friction and, hence, a vanishing diffusion coefficient. Thus, the solution we find via Darcy's Law represents the upper limit of the plateau length in the shape evolution curves. We see that all of our simulation results and the extracted experimental results are bounded by this limiting curve.
It is striking that the continuum description appropriate for highly frictional aggregate liquids connects with Darcy's Law and the study of fluids in porous media. Such systems are well-studied [32,33], but rarely in a setting where surface tension dominates the system behavior. We believe this study may serve as motivation for precise computational fluid dynamical tools to be deployed in a new problem domain. Tests of our results can further be obtained by connecting Darcy's Law to semi-2D flows in Hele-Shaw cells [34]. For example, we would expect that in the relevant Hele-Shaw geometry coalescence of viscous liquid would follow a neck growth scaling law of R n ∼ −t log(t) since, at early stages, the neck radius is smaller than or comparable to the gap width of the Hele-Shaw cell; this makes it more like the standard 3D coalescence of a simple liquid in the viscous regime, a la Hopper's solution. As the neck grows and the system becomes better described by Darcy's Law rather than the Stokes equation, we would expect a new scaling of R n ∼ t α . These two regimes may have been observed in the experiments of Ref. [35], although the interpretation is complicated by their use of a vertical Hele-Shaw cell to study a drop coalescing into a planar interface.It may be possible to tune the crossover point between the regimes of −t log t and t α by changing the ratio of drop radius and the gap between two plates (although this is complicated by the fact that ζ in the Hele-Shaw system is not an independent parameter, but is itself related to µ).
In summary, we have combined particle-based and continuum simulations to investigate the coalescence of liquid droplets evolving according to different microscopic dynamics. We demonstrate experimentally testable qualitative features that can distinguish the overall shape evolution of the system from the theoretical predictions derived in the context of Stokes flow, and further show that we expect new power laws characterizing the neck growth. Thus, even before considering additional complexities that may come from effects as diverse as biochemical regulatory communication pathways in cellular aggregates [36,37] or viscoelastic effects in polymeric coalescence [13], we have shown that strong coupling to a dissipative environment is itself a mechanism for changing the physics of coalescence. We speculate that this is relevant not only for the nucleolar systems mentioned above [14], but also for other mesoscopically large degrees of freedom for which Brownian or Langevin descriptions are frequently used (such as cell spheroids [7][8][9] and colloidal systems [38,39]). The interplay of viscosity, inertia and dissipation may also appear in other physical systems for which information propagation plays a strong role, such as in the critical role of an effective "inertia" in the propagation of turning information in bird flocks [40]. Finally, although our results focus on purely two dimensions and might have some quantitative discrepancies from three-dimensional systems, our qualitative finding of the lag between length decrease and neck growth can also be observed in the 3D simulation results with Brownian Dynamics [18].

Particle-based Simulations
In our discrete simulations, we use a 2D binary Lennard-Jones system A 65 B 35 with U (r) = 4 σ r 12 − σ r 6 , where σ AA = 1.0, AA = 1.0,σ AB = 0.8, AB = 1.5,σ BB = 0.88, BB = 0.5 and r cut = 2.5, as in Ref. [41]. The simulations are run at temperature kT = 0.35, which is in the liquid phase. We simulate three distinct equations of motion: Langevin Dynamics: Brownian Dynamics: using default values of m = 1 and λ = 1 and a time step size of δt = 0.002 in natural Lennar-Jones units. Our MD simulations are done in the canonical (NVT) ensemble with a with Nosé-Hoover thermostat. All of these particlebased simulations are run with the HOOMD-blue [42] simulation package. All the particle-based simulations are run on XSEDE [43].

Continuum simulations of the Navier-Stokes Equation
Our numerical solutions of the Navier-Stokes Equation: ρ Du Dt = ∇P + µ∇ 2 u are obtained using the open-source Basilisk software package [44] with the same parameters for inner and outer fluid: ρ = 1, γ = 1 and µ as written in the figure in the main text.

Boundary Integral Method Simulations of Darcy's Law
Due to the well-known appearance of spurious vorticities (and other numerical artifacts) in continuum simulations of systems controlled by surface tension [45], we solve Darcy's law, ∇P = −ζu, for the coalescing droplet geometry via the Boundary Integral Method [46,47]. For the case of two discs of Brownian fluid 1 with frictional coefficient ζ 1 coalescing in a background Brownian fluid 2 with ζ 2 , we can get the distribution density of Green's function dipoles on the interface as: where I(q(x 0 )) = ∇G(x 0 , x) ·n(x)q(x)dS, C = ζ 2 /ζ 1 , and κ(x 0 ) is the curvature.
For C = 1 this leads to an integro-differential equation that typically requires costly iterative methods to solve; this is especially true in the context of the highly singular initial geometry of two barely-contacting disks in our study. Based on previous theories and experiments [21,26], we expect the fundamental behavior of the system to be similar for a very broad range of innerand outer-friction coefficients and thus focus on the computationally simpler C = 1 case.
We carry out this computation using custom MATLAB 2020b code. We use 1200 points on the interface which we distribute at each update based on the arc-length ds between neighboring marker points. We use a relation such that ds ∝ x 2 + y 2 , to ensure a denser distribution of points near the neck region (where the curvature and, hence, the need for numerical stability, is highest). The interface is advanced using an adaptive stepsize fifth-order Dormand-Prince method [48], and for an initial shape we use r(θ) = 4(a − sin 2 θ) with a = 1.001, which a polar-coordinate representation of Hopper's solution for the Stokes equation [10] at a very early time.

Measuring droplet curvature in simulations
The curvature at the neck is measured based on the average contours of five simulations of Molecular Dynamics and three simulations of Brownian Dynamics. We first parametrize the contours in terms of arc length s, starting from and ending at the left side. Then we use the csaps function in MATLAB to smooth the contours with the cubic spline. The smoothing parameter is carefully tuned to make sure that the neck region is not oversmoothed. One example of the smoothing result is shown in Extended Fig. 3. The curvature is averaged over the five data points closest to the minimum neck radius point for both the upper part and lower part of the neck. Curvature at each point is calculated using the formula κ = x (s)y (s)−y (s)x (s) (x (s) 2 +y (s) 2 ) 3/2 . For the numerical solution of Darcy's law, the curvature is measured in a similar way, but without the contour average and only for the minimum neck radius point.

Estimation of crossover points in neck growth curves
For a discrete system with initial drop size R and particle size r e , we can estimate the crossover point by intersecting the κR vs. R n /R curve obtained from Darcy's law (Fig. 2d), which is κR ≈ 0.4( Rn R −3 , with the upper limit line for the unitless curvature κR. Specifically in our simulation, the initial neck width is about the size of three elements (as shown in Extended Fig. 1(a)), thus the upper limit is κR ≈ R 3re . Together these limiting forms suggest a crossover neck radius of R cross n R = 1.2re R 1/3 . The crossover points in simulations with N equaling 50000, 12500, and 8000 separately are estimated by intersecting the fitting lines for the t α regime and the t β regime, as shown in Extended Fig. 1(b-d).  Extended Fig. 2 The dependence of time scale on system size. (a) Neck growth curves for Molecular Dynamics simulations with different system sizes. The time is normalized with the initial radius R for one droplet in each case. (b) Neck growth curves for Brownian Dynamics simulations with different system sizes. The time is normalized with the cube of the initial radius R for one droplet in each case. The insets are the curves before normalization by system size.

Smoothed Contour
Original Contour Extended Fig. 3 Smoothing of the Contours