Activity drives self-assembly of passive soft inclusions in active nematics

Domain: chemistry
Published: 2026-02-27
Editorial outcome: minor_revision

Abstract
Active nematics are out-of-equilibrium systems in which energy injection at the microscale drives emergent collective behaviors, from spontaneous flows to active turbulence. While the dynamics of these systems have been extensively studied, their potential for controlling the organization of embedded soft particles remains largely unexplored. Here, we investigate how passive droplets suspended in an active nematic fluid self-organize under varying activity levels and packing fractions. Through numerical simulations, we uncover a rich phase diagram featuring dynamic clustering, activity-induced gelation, and a novel activity-driven deformability-induced phase separation regime where activity stabilizes dense droplet assemblies. We find that droplet deformability plays a key role in enabling this regime, as it allows droplets to absorb the stress exerted by the surrounding active fluid. Crucially, we demonstrate that temporal modulation of activity enables precise control over structural morphological transitions. Our results suggest new routes to design adaptive smart materials with tunable microstructure and dynamics, bridging active nematics with applications in programmable colloidal assembly and bio-inspired material design.

Introduction
Active nematics are a class of non-equilibrium systems that consume energy at the microscopic scale, driving them far from equilibrium and enabling the emergence of spontaneous flows and complex macroscopic behaviors [1–4]. At low activity levels, energy injection leads to local deformations of the active nematics, resulting in ordered flow patterns [4–6]. However, as the rate of energy injection increases, deformations strengthen, accompanied by a corresponding rise in the magnitude of active flows. This drives the system into a chaotic regime known as active turbulence, characterized by the continuous creation and annihilation of topological defects [2, 7–14].

The ability to control the chaotic active nematics dynamics is central to exploiting their potential for applications, as for instance the self-organization, assembly, and transport of suspended passive particles [15–27]. In the turbulent regime, active nematics create a unique environment where fluid dynamics and defect-mediated transport can be exploited to orchestrate the motion and organization of passive inclusions [19, 28–30]. Many-body active interactions provide an effective avenue for creating new soft composite materials, such as colloidal crystals and self-quenched glasses. Similarly, the behavior of passive soft inclusions, such as lipid vesicles or polymeric droplets, in active nematics could enable novel applications in drug delivery, microfluidics, and programmable self-assembly. Unlocking this potential requires understanding the interplay between active flows, topology, and passive inclusions. In particular, it is essential to consider how spatial confinement and topological constraints introduced by suspended particles affect emergent behaviors at different activity levels and packing fractions.

The collective organization of two-dimensional colloidal inclusions in passive liquid crystals is strongly influenced by the liquid crystal anchoring [31–36]. In the presence of anchoring (either tangential or homeotropic), each suspended particle introduces a local topological charge of +1, which the surrounding nematic absorbs by generating a total compensating charge of −1. This is typically achieved through the formation of two −1/2 topological defects, so that inclusions mutually interact through weak quadrupolar interactions [31, 37]. These interactions drive various self-assembled configurations, including linear chains of colloids [38, 39], particle aggregates [40, 41], and defect-stabilized colloidal gels [42]. Recent studies suggest that activity can overcome energy barriers and render certain excitations effectively gapless [43]. This opens intriguing possibilities for designing new classes of topological soft materials, where active forces could be harnessed to reshape defect landscapes and guide structure formation beyond equilibrium constraints.

While active droplets suspended in passive fluids [44–54] and colloidal suspensions in passive liquid crystals [55–63] have been extensively studied, the case of passive inclusions suspended in active nematics has received far less attention. Recent research demonstrated that a single passive inclusion suspended in an active nematic exhibits anomalous diffusivity, deviating significantly from Brownian motion due to the interaction between bulk defects and the inclusion-induced nematic deformations [26, 64]. This behavior is strongly dependent on the ratio of the droplet radius to the active length scale, with smaller droplets experiencing significant drifts driven by active flows.

However, whether this behavior can be harnessed for material self-assembly remains an open question, requiring systematic investigation of the dynamical and morphological properties in dense droplet suspensions. In this study, we explore the full phase diagram of a system of multiple passive soft droplets, modeled as multiphase fields, suspended in an active nematic fluid. We systematically vary both the activity level and the droplet packing fraction, revealing a rich diversity of dynamical and morphological regimes.

At intermediate activity levels, we find a dynamic clustering behavior also observed in hard colloidal suspensions [65]. At higher packing fractions, the clustered phase gives way to a gel-like phase where droplets organize into a percolating network spanning the entire system. Crucially, when activity exceeds the threshold for active turbulence, the system undergoes a striking transition: the suspended droplets completely cluster, exhibiting tunable diffusive and morphological properties. This regime resembles motility-induced phase separation (MIPS) observed in systems of self-propelled particles [1, 66–73], but with a fundamental distinction. In classical MIPS, phase separation arises from the particles’ own activity, whereas in our system, it is driven by the motility of the active solvent and by the deformability of the passive soft inclusions. Accordingly, we name this regime Activity-Driven Deformability-Induced Phase Separation, or shortly active-DIPS.

Results
Morphological Characterization. To investigate the collective behavior of passive soft droplets suspended in an active nematic fluid, we employ a multiphase field approach. Each droplet is represented by a distinct phase-separating field ϕ_i, where i uniquely labels each droplet. The orientational state of the active nematic fluid is described by the Q-tensor, defined as Q = S(nn−I/3). Here, S represents the nematic order parameter, quantifying local nematic order, and n is the head–tail symmetric unit vector indicating the local alignment direction of nematic molecules.

We numerically solve the continuum active nemato–hydrodynamic equations coupled to the conserved dynamics of the suspended passive droplets. This model has previously been validated for reproducing the dynamics of active suspensions confined on rigid and deformable shells, active liquid interfaces, and predicting patterns in three-dimensional active suspensions (see the Methods Section for a detailed description of the numerical model).

Non-equilibrium energy injection is introduced into the Navier-Stokes equation via the active stress tensor σ_act = −ζQ, where ζ is the activity parameter. In the following, we refer to the dimensionless activity number defined as A = R_d^2ζ/κ_F which quantifies the relative strength of active stress compared to the elastic stress of the liquid crystal introduced by the inclusions. Here, R_d is the droplet radius and κ_F is the elastic constant of the active liquid crystal. [Refer to the Methods for further details on the numerical model and the definition of the reference time-scale τ and energy-scale ϵ_def.]

At low activity levels (A ≲ 3), the system remains in a quiescent state dominated by passive effects. When the liquid crystal exhibits either homeotropic or tangential anchoring at the droplet interface, the director field wraps around each droplet in a way that imposes a topological charge of +1. To satisfy global topological constraints the liquid crystal compensates by generating two −1/2 defects in the surrounding medium. These defects are not rigidly attached to the droplet interface but are energetically stabilized at a distance ℓ_p = κ_F/W – the so-called penetration length where W denotes the anchoring strength at the droplet surface – and, in isolated droplets, they are typically located at antipodal positions (Fig. 1a(i)).

This configuration induces a quadrupolar short-range interaction between droplets, too weak to drive droplet attraction. As a result, the emergent structure depends sensitively on the droplet packing fraction Φ. Specifically, at low packing fraction and activity number, droplets remain near their initial positions in a quiescent stationary state (Supplementary Movie 1).

The elastic interactions mediated by the defects induced around the suspended droplets become significant when the system is sufficiently active to sustain spontaneous flows, or when the packing fraction is high enough to reduce the average inter-droplet distance to the penetration length of the liquid crystal. In these regimes, defect-mediated forces can influence droplet organization and dynamics, contributing to the emergence of collective behavior. Specifically, when two droplets approach each other, anchoring-induced deformations in the surrounding liquid crystal stabilize the formation of defect bonds. These bonds arise when droplets share two (or more) of their surrounding topological defects, which causes them to adhere (Fig. 2a,b). Such droplet dimers act as aggregation centers, facilitating the growth of larger clusters through the absorption of isolated droplets or the merging with other clusters. In particular, when the packing fraction becomes sufficiently high (Φ ≳ 0.21), the typical distance between droplets becomes small enough for elastic interactions to dominate. These short-range attractive forces are strong enough to drive a global rearrangement of the droplets, leading the system to evolve into a stationary percolating network (Fig. 1a(ii) and Supplementary Movie 2), even in the absence of activity. This behavior is indeed consistent with previous observations in colloids or droplets dispersed in passive nematics, which showed chain-like assemblies—linear colloidal chains [38, 39], particle aggregates [40, 41] and defect-stabilized colloidal gels [42].

Introducing activity into the system profoundly alters its morphology. A sharp transition occurs when the activity level surpasses the threshold for spontaneous flow (A ≳ 3). In this regime, droplets are advected by active jets and exhibit a characteristic stop-and-go motion [64]. This persistent movement causes droplets to interact with nearby droplets, forming defect-mediated bonds that lead to the emergence of small droplet clusters. Interestingly enough, at intermediate activity levels (A ≲ 200), cluster growth is arrested due to the disruptive effects of active jets, which undermine the structural stability of the clusters. This behavior is prevalently observed at low packing fractions and intermediate values of activity (Fig. 1a(iii) and Fig. 1b). Consequently, the system exhibits highly dynamic behavior, with clusters continuously assembling and disassembling over time (Supplementary Movie 3).

At higher packing fractions, activity accelerates the formation of a percolating network even at intermediate activity levels. As a result, the system undergoes a transition from dynamic clustering to a gelated state as the packing fraction increases (Fig. 1b). Specifically, when the packing fraction is sufficiently large, active jets mobilize droplets into chain-like structures that eventually intersect, forming a passive network of droplets that encapsulates regions of active nematic fluid. Notably, the droplet network confines the active liquid crystal into domains smaller than the overall system size, effectively suppressing unstable infrared active modes, leading the system to settle into a stationary state over long timescales (Supplementary Movie 2).

To characterize the system morphologically, we show in Fig. 2c-f the droplet radial distribution function g(r) (see the Methods Section) for different morphological regimes. In the gelation phase, Fig. 2c displays a sharp, well-defined first peak corresponding to the nearest-neighbor droplet distance, followed by a series of secondary peaks, reflecting short-range positional correlations typical of an extended, percolated droplet network lacking long-range order. By comparison, in the clustering phase, Fig. 2d reveals a dominant first peak corresponding to the first-neighbor contribution, along with a smaller secondary peak, consistent with the finite size of the cluster. The two-dimensional structure factor S_ϕ = ⟨ϕˆ*ϕˆ ⟩ (with ϕ = ∑_i ϕ_i the global concentration field) further exhibits broken rotational symmetry, reflecting the local arrangement of droplets within the clusters.

When the activity is increased beyond A ≳ 200 at low packing fractions (Φ ≲ 0.15), the active flow becomes large enough to prevent clustering, leading the system to transition into a new regime characterized by chaotic droplet motion (Supplementary Movie 4) and the absence of structures (isotropic concentration spectrum and flat g(r) in Fig. 2e). We identify this state as active turbulence, where the continuous active injection destabilizes the liquid crystal, resulting in deformed patterns with a characteristic length scale of ℓ_act = κ_F/|ζ|. In this regime, the liquid crystal relaxes the elastic stress caused by activity-induced deformations by nucleating pairs of topological defects. This process leads to the proliferation of topological defects throughout the system, further contributing to the active turbulence dynamics.

Remarkably, at very high activity levels, we observe droplets reorganizing into clusters once again, despite the presence of strong active jets (Supplementary Movie 5). While defect bonds are primarily responsible for cluster stability at lower activity levels, at high activity, the pressure exerted by the active liquid crystal becomes sufficiently large to deform the suspended droplets and stabilize the clusters (Fig. 1a(iv) and Fig. 1b). Accordingly, we refer to this regime as Activity-Driven Deformability-Induced Phase Separation, or shortly active-DIPS. As a result, these clusters remain fully separated, grow over time, and retain their constituent droplets without loss. This occurs because the steric constraints imposed by the surrounding active nematic matrix effectively trap the droplets in place, deforming and arranging them in hexatic fashion as shown by the g(r) and 2D structure factor in Fig. 2f.

The existence of this fully clustered state at high activity is one of the most significant findings reported in this article. Notably, this behavior is largely independent of the packing fraction of the system, highlighting the robustness of steric trapping at large activity. Finally, as illustrated in the phase diagram in Fig. 1b, a transition from a percolating network to a fully separated phase also occurs at large packing fractions.

Kinetic Characterization. Fig. 3a presents the mean squared displacement (MSD) ⟨∆r^2(t)⟩ of droplets across different phases. In the quiescent and gelation phases at low activity, droplets reach a stationary configuration at long times (Fig. 3a(i-ii)), while in the three higher-activity regimes, they exhibit long-time diffusive behavior (Fig. 3a(iii-v)). Notably, this diffusive dynamics emerges independently of packing fraction, as demonstrated by the behavior of the long-time MSD exponents versus activity (Fig. 3e), obtained by performing a statistical average over independent simulations for different packing fractions (see the Methods Section).

In general, the system undergoes a transition from sub-diffusive behavior at low activity to diffusive behavior at higher activity, with the packing fraction influencing the specific activity required for this transition; higher packing fractions reduce the available space for the active nematic instability, necessitating greater activity to generate self-sustained flows and sustain droplet diffusive motion. Notably, near the transition between gelated and clustered states, we observe significant fluctuations in the MSD exponents (Fig. 3f), evidenced by peaks in the standard error of the mean. These fluctuations suggest phase coexistence during the transition to the fully clustered phase at high activity, whereas the transition to the turbulent state at lower packing fractions appears continuous.

Interestingly, the kinetic energy of the droplets, E_d = ⟨m r˙^2/2⟩, and of the fluid, E_kin = ⟨ρv^2/2⟩, exhibit distinct behaviors as the activity level and packing fraction are varied. For a fixed packing fraction, the fluid’s kinetic energy increases monotonically with activity (Fig. 3g), consistent with the progressive amplification of active stresses. However, the dependence on packing fraction is more nuanced. At low packing fractions (e.g., Φ=0.09), and moderate activity (A=432), the system is in the active turbulent regime, with droplets dispersed throughout the domain. In this configuration, the active liquid crystal can develop instabilities over a limited length scale, roughly set by the inter-droplet distance, which restricts the formation of large-scale vortical structures.

By comparison, at the same activity but higher packing fraction (e.g., Φ=0.35), the system enters the active-DIPS phase, where droplets cluster and leave extended regions of active fluid. This phase-separated configuration allows the LC to exploit the full system size to develop active instabilities, resulting in stronger and more coherent flows. Consequently, the fluid kinetic energy can be higher at larger packing fractions, contrary to the intuition that crowding suppresses flow (see also Section Turbulent Flows in active-DIPS of the SI). At higher activity levels, even low packing fraction systems undergo phase separation into the active-DIPS regime. In this case, the lower packing fraction allows for a larger volume of active fluid, which supports stronger flows and thus higher fluid kinetic energy. In contrast, the behavior of the droplet kinetic energy E_d is non-monotonic with increasing activity (Fig. 3h). At low packing fractions, the system transitions from a turbulent regime–where droplets are advected independently by active flows–to a phase-separated state where droplets form a single cluster. In this clustered state, droplet motility decreases, as the diffusivity of the larger aggregate becomes smaller than the typical scale of uncorrelated active jets. Conversely, at high packing fractions, the system transitions directly from a gel-like state with minimal droplet motion to a phase-separated state where the cluster is still advected by the flow. However, due to its larger size, the cluster’s motility remains low, and the overall droplet kinetic energy is reduced.

Clustering Dynamics. The dynamics of active-DIPS reveals a complex interplay between active flows and passive droplet organization. As illustrated in Fig. 4a, the process begins with the formation of small, transient clusters (panels a(i-ii)) similar to those observed at lower activity levels. However, at these elevated activity levels, the system exhibits a qualitatively different behavior as clusters progressively merge into larger and more stable structures (panels a(iii-iv)). These emergent superstructures serve as nucleation centers that incorporate smaller clusters and individual droplets through an ongoing coarsening process. Notably, these macroscopic clusters eventually reach a stationary state where the internal droplets are organized into hexatic configurations.

The stabilization mechanism of these clusters differs fundamentally from the defect-mediated bonding observed at lower activities. In this high-activity regime, where the interstitial space between droplets becomes depleted of liquid crystal, hydrodynamic effects dominate the collective behavior. Cluster stabilization arises from the pressure gradient that develops as a consequence of asymmetric forces acting on peripheral droplets. While interior droplets experience balanced interactions with their neighbors, peripheral droplets are subjected to strong active stresses from the surrounding turbulent fluid. This imbalance establishes a radially decreasing hydrodynamic pressure profile P_h toward the cluster center, as clearly evidenced in Fig. 4b-c.

The dynamics of active-DIPS results from a density-dependent mobility reduction, analogously to conventional MIPS. However, in our case, passive droplets are transported by active flows and feature decreased mobility upon collisions and elastic interactions–which eventually lead to local accumulation. This process is complemented by a feedback mechanism in which accumulating droplets progressively exclude the active fluid from the cluster interior. The latter, in turn, exerts compressive stresses on the cluster boundaries, which are transmitted inward through the compression of the droplet network.

The temporal evolution of the system’s organization is quantitatively captured through spectral analysis. The concentration spectrum ⟨ϕˆ*ϕˆ ⟩ (Fig. 4d) reveals two important features: the development of Bragg peaks at large wavenumbers reflects the emergence of hexatic order within clusters, while the growing intensity at small wavenumbers indicates the progressive increase in cluster size. Furthermore, the kinetic energy spectrum E_k = ⟨vˆ* · vˆ ⟩ (Fig. 4e) maintains the characteristic scaling of active turbulence, demonstrating that the active fluid retains its chaotic nature despite the presence of the separated cluster. From these spectra, we extract the characteristic length scales of the concentration field ℓ_ϕ and velocity field ℓ_v as the first moments of their respective distributions. As shown in Fig. 4f, ℓ_v < ℓ_ϕ throughout the evolution, confirming that the turbulent vortices remains smaller than the emergent droplet structures. Such a length-scale separation highlight the hierarchical organization of the system, where microscopic active turbulence coexists with macroscopic droplet ordering. This is further exemplified by the evolution of the probability density function (PDF) of the droplet velocity shown in Fig. 4g as clustering progresses. At early times, the distribution is approximately Gaussian, aside from heavy tails that reflect intermittent behaviors typical of passive particles suspended in active turbulent nematics [74]. As clustering advances, the Gaussian behavior progressively disappears, leading to a pronounced peak around the mean velocity. This shift indicates that droplet velocities become increasingly correlated–consistent with the emergence of collective motion, where clustered droplets behave as a single body rather than as independent particles advected by the turbulent flow.

We next examine the role of surface tension in the stability of active-DIPS, as illustrated in Fig. 5. In panel (a), we perform a experiment starting from a clustered state and vary the surface tension σ. As surface tension rises, the cluster undergoes melting, which is accompanied by a sharp decrease in the average hydrodynamic pressure gradient ⟨|∇P_h|⟩R_d^3/κ_F (Supplementary Movie 6). This behavior is further detailed in the inset, showing the temporal evolution of the pressure gradient for two representative values of surface tension: low tension maintains a high pressure gradient which stabilizes the cluster, whereas high surface tension makes the droplets stiffer so that they are not able to absorb external pressure eventually leading to melting.

The phase diagram in Fig. 5b further reveals a competition between activity and surface tension; high activity stabilizes active-DIPS, while high surface tension promotes melting and drives the system toward an active turbulent state, with the transition boundary scaling approximately as σ ∼ A^1/2. Our data suggest that this scaling can be rationalized by introducing an active capillary number, defined as Ca_act = ζaR_d^2/σ^2, where a is the bulk compressibility modulus. The observed transition line σ ∼ A^1/2 ∼ ζ^1/2 corresponds to constant Ca_act, confirming that the phase behavior is controlled by the competition between active forcing and interfacial stresses that resist deformation and regulate droplet compressibility.

Discussion. Understanding how clustered and active-DIPS states respond to time-dependent activity changes is crucial for developing robust self-assembly strategies. To investigate this dynamic response, we conducted systematic numerical experiments that probe the system’s behavior under different activity quenching rates. Fig. 6 illustrates the results of this investigation. A disordered configuration with packing fraction Φ = 0.18 at the initial time is evolved by increasing activity until A=1000, to obtain a fully clustered state through active-DIPS (left panel). Starting from this configuration, we then proceed to quench the activity at various rates and observe three possible regimes.

At rapid activity quenching rates (−τdA/dt ≳ 29), the system maintains its clustered configuration through persistent defect-mediated bonds, with the interstitial nematic forming a hexatic matrix of defects between droplets (Supplementary Movie 7). More gradual activity quenching leads to different morphological behaviors. In particular, at intermediate quenching rates (3.7 ≲ −τdA/dt ≲ 29), we observe the fully separated cluster melt, as decreasing active pressure renders peripheral droplets vulnerable to the advective action of turbulent flows. This creates an intriguing coexistence regime where the initial cluster persists alongside a dispersed droplet suspension (Supplementary Movie 8). At slower quenching rates (−τdA/dt) ≲ 3.7, complete disorder is ultimately restored, with the system adiabatically returning to its initial configuration through complete melting (Supplementary Movie 9).

These rate-dependent pathways are summarized in the phase diagram of Fig. 6 (right panel), where the activity number A is plotted against normalized time t/τ, establish temporal activity modulation as a powerful tool for designing reconfigurable activity-driven composites. The ability to select between clustered, partially clustered, or melted states through control of activity dynamics suggests new possibilities for designing materials with tunable structural memory.

From a broader perspective, our work demonstrates how active nematics can orchestrate the behavior of embedded soft particles through the interplay of hydrodynamic flows and topological constraints. The observed transitions–from turbulence-driven clustering to gelation and active-DIPS–significantly broaden the toolkit for non-equilibrium self-assembly by harnessing the chaotic dynamics of active turbulence. We anticipate these findings will stimulate experimental efforts using both biological and synthetic active nematics, particularly to explore how droplet softness, anchoring conditions, and confinement influence the reported phenomenology. Our model conceptually aligns with quasi-2D experimental systems such as active nematic layers confined at a water-oil interface or soap films, where soft inclusions like LLPS droplets or hydrogel particles can be embedded and interact with active flows. These setups offer a promising route to experimentally realize and test the mechanisms uncovered in this study. The principles revealed in this study could ultimately guide the design of adaptive smart materials capable of programmable structural transitions.

Methods
The model

We consider a suspension of isotropic droplets in a quasi two-dimensional active nematic liquid crystal. Each droplet is described using a multiphase approach, modeled by phase fields ϕ_i(r,t), with i = 1,...,N_d, where N_d represents the total number of droplets in the system. The dynamics of the system is characterized by the incompressible two-dimensional velocity field v(r,t), constrained by ∇ · v = 0, and the nematic Q-tensor Q(r,t) = S(nn−I/3). Here, the principle eigenvector n – the so-called director field – defines the 3D direction of alignment of the liquid crystal, while its magnitude S quantifies the degree of local nematic order.

Evolution equations. The dynamics of the fields is governed by the following set of coupled PDEs,

D_t ϕ_i = M∇^2µ_i, (2)

D_t Q = S(W,Q) + γ^−1 H, (3)

ρ D_t v = ∇·(σ_hydro + σ_ϕ + σ_LC + σ_act), (4)

where the operator D_t = ∂_t + v·∇ denotes the material derivative. Eq. (2) is the Cahn-Hilliard equation for the conserved concentration fields ϕ_i, with M the mobility parameter, and µ_i = δF/δϕ_i the chemical potentials of each droplet, F being the free energy of the system introduced below. Eq. (3) is the Beris-Edwards equation ruling the dynamics of the Q-tensor. The operator

S(W,Q) = (ξD + Ω)(Q + I/3) + (Q + I/3)(ξD − Ω) − 2ξ(Q + I/3)Tr(QW), (5)

denotes the co-rotational derivative and defines the dynamical response of the LC to straining and shearing. Here, D = (W + W^T)/2 and Ω = (W − W^T)/2 are the symmetric and anti-symmetric part of the velocity gradient tensor W = ∂_β v_α, respectively. The flow alignment parameter ξ determines the aspect-ratio of the LC molecules and the dynamical response of the LC to an imposed shear flow. Here, we choose ξ = 0.7 to consider flow-aligning rod-like molecules. The coefficient γ in Eq. (3) is the rotational viscosity measuring the relative importance of advection with respect to relaxation, and H = −δF/δQ + (I/3)Tr(δF/δQ). Finally, Eq. (4) is the Navier-Stokes equation for the incompressible velocity field (∇·v = 0) with constant density ρ. Here, the stress tensor has been divided in: (i) a hydrodynamic contribution σ_hydro = −P_h I + η∇v where P_h is the hydrodynamic pressure ensuring incompressibility, and η = 5/3 is the shear viscosity; (ii) an interface contribution

σ_ϕ = ∑_{i=1}^{N_d} [f_i δ_{αβ} − ∂_β ϕ_i (δF/δ(∂_α ϕ_i))], (6)

where F ≡ ∫_V dr f is the free energy; (iii) a LC contribution

σ_LC = −ξH(Q + I/3) − ξ(Q + I/3)H + 2ξ(Q + I/3)Tr(QH) + Q H − H Q − ∂_α Q_γµ ∂(∂_β Q_γµ), (7)

accounting for the elastic backflow; and (iv) the active stress σ_act = −ζQ, where ζ is the activity parameter, positive for extensile systems.

Free Energy. The equilibrium properties of the system are defined by the free energy

F = ∫_V dr f = ∫_V dr (f_ϕ + f_LC + f_anch), (8)

where f_ϕ is the emulsion free energy density, f_LC is the liquid crystal free energy (accounting for both the isotropic-nematic transition and elastic deformations), and f_anch is the anchoring energy at the droplet surfaces. In particular,

f_ϕ = ∑_{i=1}^{N_d} a ϕ_i^2(ϕ_i − ϕ_0)^2 + k_ϕ/2 ∑_{i=1}^{N_d} (∇ϕ_i)^2 + ∑_{i<j} ϵϕ_i^2 ϕ_j^2, (9)

for a,k_ϕ > 0, describes phase separating fields with two minima at ϕ_i = ϕ_0 and ϕ_i = 0 that, at equilibrium, arrange into circular 2D droplets with surface tension σ = √8ak_ϕ ϕ_0^3/9, and interface thickness ξ_ϕ = √2k_ϕ/a. The interaction term proportional to ϵ > 0 describes soft repulsion pushing droplets apart when overlapping.

The LC free energy density

f_LC = A_0/2 (1 − χ(ϕ)/3) TrQ^2 − χ(ϕ)/3 TrQ^3 + χ(ϕ)/4 (TrQ^2)^2 + κ_F/2 (∇Q)^2, (10)

Here, the bulk constant A_0 > 0, and χ(ϕ) is a temperature-like parameter that drives the isotropic-nematic transition that occurs for χ(ϕ)>χ_cr = 2.7 [75], where we have denoted with ϕ = ∑_i ϕ_i the global concentration field. We choose χ = χ_0 + χ_s ϕ, with χ_s > 2.7 and χ < 0 to confine the LC outside of the droplet where ϕ=0 (see Fig. 1a). The last term in Eq. (10) proportional to the elastic constant κ_F captures the energy cost of elastic deformations in the single elastic constant approximation.

Finally, the anchoring contribution to the free energy,

f_anch = W ∑_{i=1}^{N_d} [∂_α ϕ_i Q_αβ ∂_β ϕ_i], (11)

dictates the preferred orientation of the liquid crystal at the droplet interfaces. For W < 0, homeotropic anchoring is enforced at the droplet surfaces.

Simulation protocol

We integrate the dynamics of the hydrodynamic fields in Eq. (2)-(4) in a square grid of size L_g = 320 with a predictor-corrector hybrid lattice Boltzmann approach [44, 76–78] consisting of solving Eq. (2,3) with a finite-difference algorithm implementing first-order upwind scheme and fourth-order accurate stencils for space derivatives, and the Navier-Stokes equation (Eq. (4)) through a predictor-corrector LB scheme on a D2Q9 lattice. For technical details on the lattice Boltzmann methods the reader is referred to [77]. The simulation is initialized by placing droplets within the two-dimensional domain. Each droplet concentration field, ϕ_i, is defined as a circular region of radius R_d, taking the value ϕ_i = ϕ_0 inside the droplet and ϕ_i = 0 outside. The nematic order parameter Q is initialized with S = 0 inside the droplets, while in the surrounding active fluid it is assigned a small magnitude and a random in-plane orientation. The velocity field is initially set to zero.

The key control parameters of the system are: (i) the droplets’ packing fraction Φ = N_d πR_d^2/L_g^2; and (ii) dimensionless activity number

A = R_d^2 |ζ| / κ_F = R_d^2 / ℓ_act^2, (12)

where ℓ_act = √κ_F/|ζ| is the active length-scale. Throughout the paper, results are presented in terms of the dimensionless numbers Φ and A. Furthermore, other physical quantities are normalized as follows: distances are normalized with respect to the droplet radius R_d; energies are normalized with respect to the defect energy ϵ_def = −πκ_F ln(R_d/ℓ_c), and times are normalized with respect to the time-scale of elastic distortions of the liquid crystal τ = γR_d^2/κ_F.

Simulation parameters. We have varied two control parameters: (i) the packing fraction Φ was varied between 0.02 and 0.44, by considering systems comprising a different number of droplets, namely N_d = 5, up to N_d = 100; (ii) the activity number A was varied between 0 (ζ = 0) and 2×10^3 (ζ = 9×10^−2). The other model parameters are chosen as follows. We fixed the droplet radius to R_d = 12, the mobility parameter to M = 0.1, and the rotational viscosity to γ = 1. Unless otherwise stated, the free energy parameters are a = 0.07, ϕ_0 = 2.0, k_ϕ = 0.1, ϵ = 0.05, A_0 = 0.2, κ_F = 0.01, and W = −0.04. To confine LC outside of the droplets, the parameters χ_0 = 3.0 and χ_s = −0.25 were selected. With the chosen parameters, the liquid crystal coherence length is ℓ_c = √κ_F/A_0 ≈ 0.25. This length scale characterizes how rapidly the orientational order decays near a defect core, effectively setting the characteristic defect size to approximately one grid point in the simulation. Unless otherwise stated, simulations were performed at a fixed system size of L_g = 320, where we used periodic boundary conditions along both directions.

Measurements

Radial distribution function analysis. The radial distribution function g(r) measures the probability of finding a pair of droplets separated by a distance r relative to an ideal homogeneous distribution. It is defined as

g(r) = 1/(ρ_0 N_d) ⟨∑_{i̸=j} δ(r−r_ij)⟩, (13)

where N_d is the number of droplets, ρ_0 = N_d/L_g^2 is the average number density, and r_ij is the separation distance between i- and j-th droplets computed using the minimum-image convention to account for periodic boundaries. The angular brackets ⟨·⟩ denote a time average over all time steps. Normalization by ρ_0 N_d ensures that g(r) approaches 1 in a completely homogeneous system.

Mean squared displacement analysis. The mean squared displacement (MSD) at time t is defined as

⟨∆r^2(t)⟩ = 1/N_d ∑_{i=1}^{N_d} |r(i)(t) − r(i)(t_0)|^2, (14)

where N_d is the number of droplets, r(i)(t) is the unwrapped position vector of the i-th droplet’s center of mass at time t, and t_0 = 0 is the reference time. The angular brackets ⟨·⟩ denote an average over all droplets at a given time. Short- and long-time exponents of the MSD are extracted by fitting power laws to the data.

Statistical analysis of long-time MSD exponent. To investigate the long-time dynamical behavior, we analyze the MSD across a range of packing fractions and activity levels. Specifically, we consider three packing fractions: Φ = 0.09, 0.15, 0.22. For each packing fraction, ten distinct activity levels are considered. At each combination of packing fraction and activity level, we perform n = 5 independent simulations.

For every simulation, the MSD is computed as a function of time. A power-law fit of the form MSD(t) ∼ t^α is applied in the long-time regime to extract the MSD exponent α. The exponents obtained from the five independent simulations are then used to compute the mean and the standard error of the mean (SEM) plotted in Fig 3e and Fig. 3f, respectively.

Spectral analysis. The 2D structure factor of the concentration field ϕ(r,t) is defined in Fourier space as

S_ϕ(k,t) = ϕˆ*(t)ϕˆ(t), (15)

where ϕˆ(t) denotes the Fourier transform of the concentration field, ϕˆ(k,t) = ∫ ϕ(r,t)e^{−ik·r} d^2r. Accordingly, we define the concentration spectrum ⟨ϕˆ*ϕˆ ⟩ by integrating over circular shells of constant wavenumber magnitude q = |k|

⟨ϕˆ*ϕˆ ⟩(q,t) = ∫ S_ϕ(k,t) δ(|k|−q) d^2k. (16)

This spectrum quantifies the spatial correlations among droplets across scales, identifying the dominant inter-droplet spacing and characteristic structural organization. The characteristic cluster size ℓ_ϕ is then defined as the inverse of the first moment of the spectrum

ℓ_ϕ(t) = ∫ ⟨ϕˆ*ϕˆ ⟩ dq / ∫ q ⟨ϕˆ*ϕˆ ⟩ dq. (17)

Similarly, the 2D structure factor of the kinetic energy is defined from the velocity field v(r,t) as

E_k(t) = vˆ_k*(t) · vˆ_k(t), (18)

where vˆ(k,t) = ∫ v(r,t)e^{−ik·r} d^2r. We then define the kinetic energy spectrum ⟨E_k⟩ by radial integration:

⟨E_k⟩(q,t) = ⟨vˆ* · vˆ ⟩(q,t) = ∫ E_k(t) δ(|k|−q) d^2k, (19)

and the characteristic flow structure size ℓ_v as

ℓ_v(t) = ∫ ⟨E_k⟩ dq / ∫ q ⟨E_k⟩ dq. (20)

Intermittency analysis. The intermittent dynamics of droplets are characterized by the distributions of the x-component of their instantaneous velocities. Droplet velocities are obtained from the temporal derivatives of their center-of-mass positions, and the probability density function (PDF) of the x-component, v_x, is used to assess deviations from Gaussian behavior. PDFs are collected over all droplets across different time intervals to reveal the statistical features of velocity fluctuations.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In the paper “Activity drives self-assembly of passive soft inclusions in active nematics”, Akduman & Sariyar et al report on\nsimulations of passive droplets embedded within an active nematic. They map a phase diagram of the droplet assemblies\nover various activity levels and packing fractions. They uncover two mechanisms by which the passive droplets interact – at\nlow activity levels the droplets adhere through liquid-crystal-mediated interaction; at high activity levels the droplets cluster\nthrough activity-driven pressure. They thoroughly characterize the morphological and kinetic behavior of the droplets.\nOverall, this is intriguing work and I support its publication provided the points below are addressed.\nMajor comments:\n1. The authors use compressible droplets as the passive constituents in this paper. I found this to be an unintuitive choice\nwhich they explain at the end of the Results section: “…This explains the crucial role of the droplet compressibility in the\ncluster stability [our observations indicate that increasing surface tension and thus reducing compressibility significantly\ndestabilizes droplet clusters.” This seems like an important result and should be supported with data. If the proposed\nmechanism is correct, I would expect the gradients in Figures 4b and 4c to become shallower as surface tension increases.\nAdditionally, this behavior contrasts with MIPS, which supports stable cluster formation in large, incompressible colloids.\n2. Figure 3g shows the fluid kinetic energy vs activity for different packing fractions. In the text the authors comment that “the\nfluid’s kinetic energy increases monotonically with activity…” for a given packing fraction. This seems true. However, the\nauthors do not comment on the behavior at different packing fractions. Based on the previous paragraph (“higher packing\nfractions reduce the available space for the active nematic instability, necessitating greater activity to generate self-sustained\nflows and sustain droplet diffusive motion”). I would expect there to be some frustrated states where the fluid kinetic energy is\nlow because there is not enough space for it to undergo the bend instability to generate flows. So, for a given activity (active\nlength scale) I would have expected that the fluid energy would decrease with packing fraction. If there is any effect, it seems\nto be the opposite of this intuition (looking at activity = 0.7, for example). Can the authors discuss the behavior of the fluid\nkinetic energy at different packing fractions, especially with how it relates to their explanation of droplet MSDs?\n3. I found \"inverse MIPS” to be confusing terminology. Naively, I would expect “inverse MIPS” particles to cluster in regions\nof high activity rather than low; or perhaps that the motility induces dissolution. Additionally, classic MIPS is a dry theory, so\nit is also confusing that these are droplets in a fluid. After reading the paper, I understand the term, but it was not clear by the\nend of the abstract. I suggest something like “fluid-driven MIPS” to highlight the major differences from classic MIPS better.\nMinor comments:\n1. The ordering and labeling of figure panels is often unclear. In particular, the labels in Figures 1d–e, 3g–h, 4f, and 5f are\ndifficult to find. Please revise to ensure that figure panels follow a consistent alphabetical (a–z) order in the figure layout,\ncaption, and main text.\n2. What do the colors represent in Figures 1a, 1b, 1d, and 1e?\n3. The concept of “defect bonds” was difficult to understand on first reading. I recommend adding a schematic of an\nunbonded droplet in Figure 2, with the surrounding defects highlighted in orange, as in Figures 2a and 2b.\n4. Figure 2c and Movie S3 are not discussed in the main text. Please include commentary explaining their relevance.\n5. In Figures 4d and 4e, it is difficult to discern temporal evolution. Consider reducing the number of time points shown, or\nusing semi-transparent lines and removing the dots to improve clarity.\n6. I see a shift upward in the kinetic energy spectra at long time in 4e. Is this just noise? Or is it related to the smaller space\nthat the fluid has to mix in?\n7. The schematic in Figure 5f is difficult to interpret, particularly the angle depictions on the far right. Please consider\nclarifying or simplifying this panel.\n8. If available, please include citations for colloids embedded in passive liquid crystal transitions from quiescence to\ngelation, if available. I would expect that the first few rows of the phase diagram in Figure 1c to have been seen previously in\nboth experiment and theory? If not seen previously, is it because it the region of phase space has not been explored, or is\nthere something fundamentally different about these simulations -- either the droplet instead of colloid or the small amount of\nactivity.", "reviewer_2": "The paper describes the dynamical behavior of passive droplets in active nematics. The large dataset exploring different\npacking fraction and activity has been created using Lattice-Boltzmann simulation. The authors show a very rich dynamics\nwith interesting phase changes. The novelty of the paper concerns the effects of the active particles on the passive droplets\nshowing results that have theoretically and possibly applicative impacts. I believe that this is a nice work that could be\naccepted for publication on NATCOMM, however I have a couple of aspects that I believe should be considered.\n*) The periodic box size used in the simulations can be small for some cases. E.g. we see a typical cluster size of the order\nof 130 with the half-box size at 160. Can the authors perform one simulation for this case with a larger box in order to assess\nthe sensitivity to the box size? It should be possibile that the clustering dynamics can feel in this case the finite periodic\nlength scale of the simulation.\n*) One of the most peculiar results of turbulent flows is the multiscale nature and the intermittency induced by sparse strong\nstructure. Recently it has been shown intermittency using the velocity increments pdf in active turbulence (Yang, Zhu\nNatCom24). Can the authors provide an intermittency analysis on the passive droplets behavior in active turbulence?", "reviewer_3": "The authors present a numerical study of passive soft droplets embedded in an active nematic fluid, revealing multiple\nmorphological and dynamical regimes, including a novel inverse motility-induced phase separation state. Their phase\ndiagram, kinetic analysis, and demonstration of temporal control over morphological transitions are significant contributions.\nHowever, the manuscript remains largely descriptive, offering limited physical interpretation of the mechanisms driving the\nobserved transitions. A deeper discussion, either through scaling arguments or through comparison of key parameters such\nas activity, droplet size, and deformability. would significantly enhance the impact and generality of the work.\nFor example, the transitions between quiescent, gel, percolated, and clustered states appear to depend on balances\nbetween active stresses, ‘softness’/ elasticity of droplets, and packing, but these are not quantitatively explored. The onset of\nspontaneous flow at activity 2×10^-4 is noted, yet its physical significance is not clear. Typically, in MIPS papers, phase\ndiagrams are related to the Peclet number to help the reader with intuition. If this is not defined, then the active regimes\nshould be described as to why this is the threshold, and how it sets the scale for the rest of the activity regime. Additionally,\nto relate the novelty here, it would be helpful to express activity and relate it to the inclusions in terms of a characteristic time\nor length, how many times a droplet traverses its own diameter per unit time, and to compare directly to droplet size and\ndefect spacing across the different regimes. Similarly, the role of softness is acknowledged but not characterized. Therefore,\nseveral foundational questions remain: Is deformability a source of symmetry-breaking essential for defect bonding or\nclustering? How does it modify the effective interaction between droplets? Is this phase diagram specific to soft inclusions, or\nmight similar behavior arise with rigid particles and have not been studied in these activity regimes? Do the droplets exhibit\ninternal flows that couple to the surrounding nematic? How does the droplet-to-defect size ratio influence the nature and\nonset of the various transitions? The distinction between inverse MIPS and classical MIPS requires clarification. Could the\nobserved behavior be interpreted as conventional MIPS, viewed from the co-moving frame of the active medium?\nAddressing these points more systematically would provide important physical intuition and broaden their relevance across\nactive and soft matter systems.\nWhile the authors focus on the inverse MIPS regime at high activity and moderate packing, I would encourage greater\nattention to the intermediate percolated network phase, which seems capable of temporally pinning topological defects and\nsuppressing active turbulence. This regime could be particularly promising for achieving structural control in active\nmaterials.\nadditional comments\n1. I would like to emphasize this point again because I find it to be one of the most interesting and potentially impactful\naspects of the study. The choice to model the inclusions as soft, deformable droplets is central to the novelty, but the\nimplications are not fully explored. The only mention I could find in the text is: This explains the crucial role of droplet\ncompressibility in cluster sta-bility [our observations indicate that increasing surface tension–and thus reducing\ncompressibility–significantly destabilizes droplet clusters]. This should be elaborated quantitatively. What key physical\nphenomena arise specifically due to deformability? Can the authors change the surface tension or viscosity of the droplets?\nWould similar clustering or inverse MIPS behavior be observed with rigid particles, or is the deformation necessary? What is\nthe internal velocity field? A more in-depth discussion of how softness contributes to defect trapping, organization, or phase\nseparation would significantly improve the manuscript and clarify the broader relevance of the results. Additionally in figure\n1-2-3 the deformations are not visible but in figure 4 they are. Are the conditions the same?\n2. While the authors focus on inverse MIPS, the percolated network regime at intermediate activity and high packing fraction\nis also highly intriguing. It appears capable of temporarily pinning defects and suppressing turbulence. Do the droplets\ndeform? I encourage the authors to elaborate on the structure and dynamics of this regime and its potential relevance to\nmaterials with embedded topological control.\n3. Each droplet introduces a +1 charge, leading to two +½ defects. This generates interactions. Please provide more\nquantitative context: what is the typical distance between a droplet and its defects? Are defect locations stable over time, or\ndo they reorganize? How do the defects compensate for the charge of the deformed droplet?\n4. In Fig. 2a, the defects appear significantly smaller than the droplets. Is this scale difference consistent across activity\nlevels? How defect size or separation change with droplet size?\n5. What is the flow in the vesicles? Do the droplets respond purely to external flows, or do they exhibit internal velocity or\nfield gradients? This could influence how they deform or interact with defects. I think the paper would benefit from showing\nhow a single soft particle behaves at different activity levels, specifically in terms of its deformation, the flow fields inside and\noutside the particle, and its mean squared displacement.\nMinor:\n1. The inverse MIPS terminology is compelling but potentially confusing. In classical MIPS, phase separation is driven by\nself-propelled particle motility. Here, the droplets are passive, and it is the surrounding fluid that is active. The authors\nshould, in my personal opinion, name the effect differently to distinguish between the options: a clustering regime\nviewed from a passive-particle perspective or, as the name implies, the breaking of MIPS at high activity. And incorporate the\ndeformability of the particles in the terminology.\n2.It would be helpful to include average velocity fields for each phase. Visualize how droplets interact with active flows and\nwhether there are correlations between flow structures and droplet clusters. Maybe in Fig. 1 or Fig. 3.\n3. Figure 1 would benefit from clearer labeling and a more intuitive layout. Presenting representative snapshots for all\nphases in a single row or panel with labels, with consistent activity and packing fraction annotations.\n4. The model is 2D, but many soft biological systems, such as vesicles and LLPS condensates, are inherently 3D, and\nexperimental nematics fields with +-1/2 defects are 2D. The dimensional mismatch must be discussed. What real systems\ndoes this model best represent? For example, might LLPS droplets be in 2D planar confinement?\n5. Instead of plotting full MSD curves for every regime (Fig. 3a), a more compact representation showing the short-time and\nlong-time exponents vs activity would be more informative. As done for the late times, the MSDs could be collapsed on the\nsame figure for easy comparison with boundaries between 1-2. This would make it easier to pinpoint transitions\nbetween diffusive and turbulent regimes."}, "author_response": "NCOMMS-25-27868A\nActivity drives self-assembly of passive soft inclusions in active nematics\nResponse to Reviewers’ Letters\nReviewer 1\nIn the paper “Activity drives self-assembly of passive soft inclusions in active nematics”, Akduman\n& Sariyar et al report on simulations of passive droplets embedded within an active nematic.\nThey map a phase diagram of the droplet assemblies over various activity levels and packing\nfractions. They uncover two mechanisms by which the passive droplets interact – at low activity\nlevels the droplets adhere through liquid-crystal-mediated interaction; at high activity levels the\ndroplets cluster through activity-driven pressure. They thoroughly characterize the morphological\nand kinetic behavior of the droplets. Overall, this is intriguing work and I support its publication\nprovided the points below are addressed.\nReply: We sincerely thank the reviewer for their careful reading of our manuscript and for their positive\nassessment of our work. We are grateful for the constructive comments and suggestions, which have\nhelped us improve the clarity and quality of our manuscript. Below, we address each point raised in\ndetail.\nReviewer Comment1.1 — The authors use compressible droplets as the passive constituents\nin this paper. I found this to be an unintuitive choice which they explain at the end of the Results\nsection: “...This explains the crucial role of the droplet compressibility in the cluster stability [our\nobservations indicate that increasing surface tension and thus reducing compressibility significantly\ndestabilizes droplet clusters.” This seems like an important result and should be supported with\ndata. If the proposed mechanism is correct, I would expect the gradients in Figures 4b and 4c to\nbecome shallower as surface tension increases. Additionally, this behavior contrasts with MIPS,\nwhich supports stable cluster formation in large, incompressible colloids.\nReply: WethanktheReviewerforraisingthisimportantpoint. Ourchoicetousecompressibledroplets\nwas motivated by the desire to explore the behavior of soft inclusions in active fluids.\nExperimental\nrealizations include micron-scale beads in soap films driven by bacterial activity [29], or lipid vesicles\nsuspended in microtubules suspensions activated with molecular motors – a system which is currently\nunder investigation by some experimental groups. These systems are inherently deformable, and our\naim is to understand how their softness influences collective behavior under active stresses. We are\nalso encouraged by the interest our work has generated among various experimental groups currently\nstudying flexible filaments and soft inclusions in active media.\nThat said, we agree with the Reviewer that the original manuscript did not sufficiently elaborate on\nthe role of surface tension and compressibility. To address this, we have now added a new paragraph in\nthe Section Clustering dynamics, a new figure (Fig. 5 of the revised manuscript), and a Supplementary\nmovie (Movie S6) that systematically explore how increasing surface tension – thus reducing droplet\ncompressibility – affects the stability of clusters. These new results clearly show that deformability\nis an essential ingredient for phase separation: compressible droplets can absorb the active pressure,\nstabilizing clusters, whereas incompressible ones cannot, leading to the destabilization of the separating\nmechanism and melting of the cluster. Furthermore, the Reviewer is correct in anticipating that the\npressure gradients in Fig. 4b becomes shallower with increasing surface tension. This is indeed what\nwe observe and now explicitly show in Fig. 5a of the revised manuscript where we plot the average\npressure gradient for varying surface tension. We further present a phase diagram as a function of\nsurface tension and activity number (Fig. 5b), which demonstrates that increasing activity stabilizes\nactive-DIPS, whereas increased surface tension favors melting and drives the system toward the active\nturbulent regime. The transition boundary approximately follows the scaling relation σ ∼ A1/2. Our\nresults indicate that the observed scaling can be interpreted through the introduction of an active\ncapillary number, defined as\nζaR2\nCa = d, (1)\nact σ2\n√ √\nwhere a denotes the bulk compressibility modulus. The transition line σ ∼ A ∼ ζ corresponds to a\nconstant Ca , highlighting that the phase behavior is governed by the balance between active stresses\nact\nand interfacial forces which resist deformation and regulate droplet compressibility.\nFinally, we appreciate the Reviewer’s observation regarding the comparison with MIPS. In this re-\ngard, we anticipate that although our initial intention was to draw an analogy and refer to the observed\nphenomenon as “inverse MIPS”, we now recognize that this terminology may be misleading. Following\nthe Reviewers’ suggestion, we have renamed this phase activity-driven deformability-induced phase sep-\naration (abbreviated as active-DIPS) to avoid confusion and to more accurately reflect the underlying\nphysical mechanism [see also response to Reviewer Comment 1.3].\nReviewer Comment1.2 — Figure 3g shows the fluid kinetic energy vs activity for different\npacking fractions. In the text the authors comment that “the fluid’s kinetic energy increases\nmonotonically with activity...” for a given packing fraction. This seems true. However, the\nauthors do not comment on the behavior at different packing fractions. Based on the previous\nparagraph (“higher packing fractions reduce the available space for the active nematic instability,\nnecessitatinggreateractivitytogenerateself-sustainedflowsandsustaindropletdiffusivemotion”).\nIwouldexpecttheretobesomefrustratedstateswherethefluidkineticenergyislowbecausethere\nis not enough space for it to undergo the bend instability to generate flows. So, for a given activity\n(active length scale) I would have expected that the fluid energy would decrease with packing\nfraction. If there is any effect, it seems to be the opposite of this intuition (looking at activity\n= 0.7, for example). Can the authors discuss the behavior of the fluid kinetic energy at different\npacking fractions, especially with how it relates to their explanation of droplet MSDs?\nReply: We thank the Reviewer for this thoughtful observation and we agree that the behavior of the\nfluid kinetic energy as a function of packing fraction deserves further discussion.\nFor the lowest packing fraction considered in Fig. 3g (Φ = 0.09, N = 5) and moderate activity\nd\n(A = 432, ζ = 0.03), the system is already in the active turbulent regime. In this case, droplets are\ndispersed throughout the domain, and the active liquid crystal can develop instabilities over a length\n√\nscale comparable to the typical inter-droplet distance, which can be estimated as ∼ 0.5L / N ≈ 70,\ng d\nwhere L is the grid size and N is the number of droplets. This limited length scale restricts the\ng d\ndevelopment of large-scale vortical structures. In contrast, at the same activity but higher packing\nfraction (Φ = 0.35), the system enters the active-DIPS phase (formerly referred to as inverse-MIPS\n[see response to Reviewer Comment 1.3]). In this phase-separated configuration, the liquid crystal has\naccess to the full system size to develop instabilities, leading to stronger flows. As a result, the fluid\nkinetic energy can be higher at larger packing fractions, contrary to the intuition that crowding would\nsuppress flow.\nThis behavior is also evident in Fig. R1 of this referee report, where we show the temporal evolution\nof the active force for a case with Φ = 0.22 (N = 50) and A = 864 (ζ = 0.06). Initially, the system\nd\nis dispersed and exhibits localized regions subject to active energy injection resulting in weak flows. As\nphaseseparationproceedsandtheactivefluidisfreedbypassiveinclusions,large-scaleactiveinstabilities\nand vortical structures emerge, ultimately leading to a significant increase in fluid kinetic energy.\nInterestingly, the behavior of the droplet kinetic energy follows an opposite trend. At low packing\nfractions, the system transitions from a turbulent regime where droplets move independently under\nthe influence of active flows, to a phase-separated state where droplets form a single cluster. In this\nclustered state, droplet motility decreases, as the diffusivity of the larger aggregate becomes smaller\nthan the typical scale of uncorrelated active jets. Conversely, at high packing fractions, the system\ntransitions directly from a gel-like state with minimal droplet motion to a phase-separated state where\nthe cluster is advected by the active flow featuring diffusive dynamics. However, due to its larger size,\nthemotilityoftheclusterremainslowerthanthatofsmallerclustersobservedatlowerpackingfractions,\nand consequently, the overall droplet kinetic energy remains reduced.\nWehavenowaddedtheseconsiderationsintheSectionKineticCharacterizationwhereweextensively\ncomment on the fluid and droplet kinetic energy trends. Furthermore, we have included Fig. R1 in the\nSupplementary Information, where it appears as Fig. S2.\nReviewer Comment1.3 — I found ”inverse MIPS” to be confusing terminology. Naively, I\nwould expect “inverse MIPS” particles to cluster in regions of high activity rather than low;\nor\nperhaps that the motility induces dissolution. Additionally, classic MIPS is a dry theory, so\nit is also confusing that these are droplets in a fluid. After reading the paper, I understand the term,\nbut it was not clear by the end of the abstract. I suggest something like “fluid-driven MIPS” to\nhighlight the major differences from classic MIPS better.\nReply: We thank the Reviewer for sharing their concerns regarding the terminology used in the initial\nversion of our manuscript to describe the clustering behavior observed at high activity. We agree that\nthe term inverse MIPS may be misleading and open to misinterpretation. As the Reviewer rightly points\nout, the term might suggest motility-induced dissolution, which does not accurately reflect the behavior\nwe observe.\nFollowing the valuable feedback from both Reviewer 1 and Reviewer 3 [see also response to Re-\nviewer Comment 3.8], as well as input from colleagues, we have revised the terminology throughout\nthe manuscript. We now refer to this phase as Activity-Driven Deformability-Induced Phase Separation\n(abbreviated as active-DIPS). We believe this new nomenclature more accurately captures the essential\nfeatures of the phenomenon: its non-equilibrium (activity-driven) nature, the central role of droplet de-\nformability, and the emergence of phase separation. We hope this revised terminology improves clarity\nand avoids confusion with standard MIPS scenarios.\nReviewer Comment1.4 — The ordering and labeling of figure panels is often unclear. In\nparticular, the labels in Figures 1d–e, 3g–h, 4f, and 5f are difficult to find. Please revise to ensure\nthat figure panels follow a consistent alphabetical (a–z) order in the figure layout, caption, and\nmain text.\nReply: We have redesigned Figure 1 and Figure 6 (Figure 5 in the initial version of the manuscript)\nto improve their readability and align with the reviewers’ comments. Additionally, we have standardized\nthe placement of panel labels across all figures: labels now consistently appear in the top-left corner of\neach panel, and the alphabetical order has been carefully maintained throughout the manuscript.\nReviewer Comment1.5 — What do the colors represent in Figures 1a, 1b, 1d, and 1e?\nReply: We thank the referee for this question. In the new Figure 1a(i-iv), the colors represent the\nconcentration field ϕ. Specifically, purple indicates the active fluid, while yellow corresponds to the\ninterior of the passive droplets. A threshold has been applied such that regions with ϕ > ϕ0 are shown\n2\nin yellow, representing the droplets, and regions with ϕ < ϕ0 are shown in purple, representing the\n2\nsurrounding active fluid. The positions of topological defects are highlighted in white. We have now\nclarified this in the caption of Figure 1 in the revised manuscript.\nReviewer Comment1.6 — The concept of “defect bonds” was difficult to understand on first\nreading. IrecommendaddingaschematicofanunbondeddropletinFigure2, withthesurrounding\ndefects highlighted in orange, as in Figures 2a and 2b.\nReply: A representation of an isolated droplet surrounded by two −1/2 defects is already shown in\nthe inset of Figure 1a(i). We have now enhanced this figure by highlighting the defects, using the same\nvisualstyleasinFigures2aand2b, toensureconsistencyandimprovereadability. Furthermore, wehave\nrevised the manuscript text to provide a clearer explanation of defect bonds and now include explicit\nreferencestotheinsetinFigure1aandtoFigures2aand2bwherebondedandunbondedconfigurations\nare illustrated.\nBelow we report the changes made to the manuscript. When introducing the defect configuration\nof an isolated droplet, we now write:\nWhen the liquid crystal exhibits either homeotropic or tangential anchoring at the droplet\ninterface, the director field wraps around each droplet in a way that imposes a topological\nchargeof+1. Tosatisfyglobaltopologicalconstraintstheliquidcrystalcompensatesbygen-\neratingtwo−1/2defectsinthesurroundingmedium. Thesedefectsarenotrigidlyattached\n(cid:112)\nto the droplet interface but are energetically stabilized at a distance ℓ = κ /W – the\np F\nso-called penetration length where W denotes the anchoring strength at the droplet surface\n– and, in isolated droplets, they are typically located at antipodal positions (Fig. 1a(i)).\nFurthermore, when discussing the notion of an elastic bond, we now write:\nThe elastic interactions mediated by the defects induced around the suspended droplets be-\ncomesignificantwhenthesystemissufficientlyactivetosustainspontaneousflows, orwhen\nthe packing fraction is high enough to reduce the average inter-droplet distance to the pen-\netration length of the liquid crystal. In these regimes, defect-mediated forces can influence\ndroplet organization and dynamics, contributing to the emergence of collective behavior.\nSpecifically, when two droplets approach each other, anchoring-induced deformations in the\nsurrounding liquid crystal stabilize the formation of defect bonds. These bonds arise when\ndroplets share two (or more) of their surrounding topological defects, causing them to ad-\nhere (see Fig. 2a,b). Such droplet dimers act as aggregation centers, facilitating the growth\nof larger clusters through the absorption of isolated droplets or the merging with other\nclusters. In particular, when the packing fraction becomes sufficiently high (Φ ≳ 0.21), the\ntypicaldistancebetweendropletsbecomessmallenoughforelasticinteractionstodominate.\nThese short-range attractive forces are strong enough to drive a global rearrangement of the\ndroplets, leading the system to evolve into a stationary percolating network (see Fig. 1a(ii)\nand Movie S2), even in the absence of activity.\nReviewer Comment1.7 — Figure 2c and Movie S3 are not discussed in the main text.\nPlease\ninclude commentary explaining their relevance.\nReply: We thank the referee for this point. We have now revised the manuscript to include a\ndiscussion of both Figure 2c and Movie S2 [previously Movie S3] in the main text.\nReviewer Comment1.8 — In Figures 4d and 4e, it is difficult to discern temporal evolution.\nConsider reducing the number of time points shown, or using semi-transparent lines and removing\nthe dots to improve clarity.\nReply: We have followed the advice of the Reviewer and used semi-transparent lines, and reduced the\ndots’ size to improve the graphs readability.\nReviewer Comment1.9 — I see a shift upward in the kinetic energy spectra at long time in\n4e. Is this just noise? Or is it related to the smaller space that the fluid has to mix in?\nReply: The upward shift in the kinetic energy spectrum at long times in Fig. 4e is indeed a meaningful\nfeature and not merely noise. As the Reviewer correctly noticed, this shift reflects two key aspects of\nthe system’s evolution. First, the overall kinetic energy of the active fluid increases, as indicated by the\nupward shift. Second, the dominance of low-wavenumber modes suggests the emergence of coherent\nvortical structures at larger length scales, consistent with the trend observed in Fig. 4f, where the active\nlength scale ℓ increases over time until the fully clustered state is reached. These observations are\nv\nconsistent with the fact that, as droplets cluster, active turbulence develops, strengthening the intensity\nof the flow and fostering the formation of vortical structures in the velocity field.\nThis behavior can be rationalized by considering the space available to the active fluid as droplet\nclustering progresses. Initially, when droplets are randomly scattered, the effective confinement size\nof the active liquid crystal is approximately the average distance between droplets, ℓ ≪ L , where\nd g\nL denotes the linear dimensions of the simulation box. As droplets aggregate into a single cluster,\ng\nthis effectively frees up space for the active fluid to develop instabilities. As a result, the final state\nbecomes more turbulent than the initial one. This interpretation is further supported by the time-lapse\nvisualization of the active force field (see Fig. R1), which shows increasingly large-scale flow structures\nforming over time. We have added Fig. R1 and a short discussion of this point in the Supplementary\nsection Turbulent Flows in active-DIPS, where it appears as Fig. S2.\nReviewer Comment1.10 — The schematic in Figure 5f is difficult to interpret, particularly\nthe angle depictions on the far right. Please consider clarifying or simplifying this panel.\nReply: We appreciate the Reviewer’s feedback. Following this suggestion, we have redesigned Fig. 6\n(previously Fig. 5) with a simplified layout to improve clarity. We believe that the new version provides\na cleaner representation of the activity-quenching results, making the schematic easier to interpret.\nReviewer Comment1.11 — If available, please include citations for colloids embedded in\npassive liquid crystal transitions from quiescence to gelation, if available. I would expect that the\nfirst few rows of the phase diagram in Figure 1c to have been seen previously in both experiment\nand theory? If not seen previously, is it because it the region of phase space has not been explored, or is\nthere something fundamentally different about these simulations -- either the droplet instead of colloid or the small amount of\nactivity.\nReply: Wethanktherefereeforhighlightingthispoint. Thebehavioratloworzeroactivityisconsistent\nwith previous observations in colloids or droplets dispersed in passive nematics, which showed chain-like\nassemblies—linear colloidal chains [38,39], particle aggregates [40,41], and defect-stabilized colloidal\ngels [42]. These works were briefly cited in the Introduction; we have now expanded the literature\nreview to make this connection explicit. We have also added a clarifying remark in the discussion of the\nphase diagram (Fig. 1c) to emphasize the agreement with the passive limit at very low activity.\nReviewer 2\nThe paper describes the dynamical behavior of passive droplets in active nematics. The large\ndataset exploring different packing fraction and activity has been created using Lattice-Boltzmann\nsimulation. The authors show a very rich dynamics with interesting phase changes. The novelty of\nthe paper concerns the effects of the active particles on the passive droplets showing results that\nhave theoretically and possibly applicative impacts. I believe that this is a nice work that could be\naccepted for publication on NATCOMM, however I have a couple of aspects that I believe should\nbe considered.\nReply: Wethankthereviewerfortheirpositiveconsiderationofourwork,foracknowledgingtheimpact\nof our results, and for recommending publication. We appreciate the constructive feedback and have\ncarefully addressed the points raised. Below, we provide detailed responses to each of the Reviewer’s\ncomments.\nReviewer Comment2.1 — The periodic box size used in the simulations can be small for some\ncases. E.g. we see a typical cluster size of the order of 130 with the half-box size at 160. Can the\nauthors perform one simulation for this case with a larger box in order to assess the sensitivity\nto the box size? It should be possible that the clustering dynamics can feel in this case the finite\nperiodic length scale of the simulation.\nReply: We thank the referee for raising this important point. Following the Reviewer’s suggestion, we\nperformed additional simulations with varying box sizes, ranging from L = 160 to 640, while keeping\ng\nthe activity number (A = 864) and packing fraction (Φ = 0.22) fixed to match the specific case shown\nin Fig. 4 of the main text. Figure R2 reports the results for selected grid sizes. These tests show no\nqualitative differences in the final configurations or in the nature of the clustering dynamics, with the\nonly exception of the smallest size considered.\nWe agree with the Reviewer that finite-size effects may become relevant in specific regions of the\nphase diagram, especially at large packing fractions. While our simulations do not show qualitative\nchanges in the clustering behavior for the specific values of activity and packing fraction investigated,\nwe acknowledge that the precise location of the transition line between the gelation regime and the\nactive-DIPS phase may be subject to quantitative shifts in very dense systems due to the limited space\navailable for droplet rearrangement. In such cases, droplets may form elongated clusters that reconnect\nthrough periodic boundaries (as illustrated in the left panel of Fig. R2). This could lead to structure\nfactors and pair correlation functions g(r) exhibiting features of a percolated state, even though the\nunderlying dynamics remain consistent with the active-DIPS regime.\nThis behavior, however, is limited to regions of the phase diagram near the transition between the\nactive-DIPS and gelation phases. We stress that a thorough characterization of this effect would require\nreconstructing the phase diagram at multiple system sizes, which is computationally prohibitive. For\nreference, each simulation used to build the phase diagram in Fig. 1 requires approximately four months\nof machine time, and scaling up to larger systems would scales as L2 due to both increasing number of\ng\nparticles and grid size.\nWe have added Fig. R2 and a discussion of this point in the Supplementary section System Size\nDependence of active-DIPS (see Fig. S5), to clarify the potential limitations of our results.\nReviewer Comment2.2 — One of the most peculiar results of turbulent flows is the multiscale\nnature and the intermittency induced by sparse strong structure. Recently it has been shown\nintermittency using the velocity increments pdf in active turbulence (Yang, Zhu NatCom24). Can\ntheauthorsprovideanintermittencyanalysisonthepassivedropletsbehaviorinactiveturbulence?\nReply: We appreciate the Reviewer’s suggestion, which motivated us to carry out the following ad-\nditional analysis. We have performed an intermittency analysis of the passive droplet behavior in ac-\ntive turbulence, focusing on the droplet velocity component v and the longitudinal velocity difference\nx\n∆v∥(r). Our findings confirm that passive droplets inherit the multiscale and intermittent features of\nthe active turbulent flow, in agreement with the results reported by Yang et al. [74].\nIn particular, we found this analysis particularly insightful, as it can be interpreted as a fingerprint\nof the transition from active turbulence to the active-DIPS regime. Accordingly, we have revised Fig. 4\nof the main text to include this analysis in panel (g) of the updated figure. The newly added Fig. 4g\nillustrates the hierarchical organization of the system, through the transition from active turbulence to\nclustering in active-DIPS. At early times, when droplets move chaotically under the advective action\nof turbulent dynamics, the PDF of the droplet velocity component v is approximately Gaussian, aside\nx\nfrom heavy tails that reflect intermittent behaviors typical of passive particles suspended in active\nturbulent nematics. As clustering advances, the Gaussian behavior progressively disappears, giving way\nto a pronounced peak around the mean velocity. This shift indicates that droplet velocities become\nincreasingly correlated—consistent with the emergence of collective motion, where clustered droplets\nbehave as a single coherent entity. Additionally, Fig. S3 has been included in the Supplementary\nInformation (see Fig. R3 of this reply), presenting the PDFs of the x-component of droplet velocities\nand of longitudinal velocity increments for two separations, R = 2R and R = 12R . This figure\nd d\nillustrates strong intermittency at small scales and its decay at larger separations, consistent with the\nobservations reported by Yang et al. in Ref. [74]. A new section has been added to the Supplementary\nInformation to comment on the intermittency analysis in the active turbulence regime.\nReviewer 3\nTheauthorspresentanumericalstudyofpassivesoftdropletsembeddedinanactivenematicfluid,\nrevealing multiple morphological and dynamical regimes, including a novel inverse motility-induced\nphase separation state. Their phase diagram, kinetic analysis, and demonstration of temporal control\nover morphological transitions are significant contributions. However, the manuscript remains\nlargely descriptive, offering limited physical interpretation of the mechanisms driving the observed\ntransitions. A deeper discussion, either through scaling arguments or through comparison of key\nparameters such as activity, droplet size, and deformability. would significantly enhance the impact\nand generality of the work. For example, the transitions between quiescent, gel, percolated, and\nclustered states appear to depend on balances between active stresses, ‘softness’elasticity of\ndroplets, and packing, but these are not quantitatively explored. The onset of spontaneous flow at\nactivity2×10−4 isnoted, yetitsphysicalsignificanceisnotclear. Typically, inMIPSpapers, phase\ndiagrams are related to the Peclet number to help the reader with intuition. If this is not defined,\nthen the active regimes should be described as to why this is the threshold, and how it sets the\nscale for the rest of the activity regime. Additionally, to relate the novelty here, it would be helpful\nto express activity and relate it to the inclusions in terms of a characteristic time or length, how\nmany times a droplet traverses its own diameter per unit time, and to compare directly to droplet\nsize and defect spacing across the different regimes. Similarly, the role of softness is acknowledged\nbut not characterized. Therefore, several foundational questions remain: Is deformability a source\nof symmetry-breaking essential for defect bonding or clustering? How does it modify the effective\ninteraction between droplets? Is this phase diagram specific to soft inclusions, or might similar\nbehavior arise with rigid particles and have not been studied in these activity regimes? Do the\ndroplets exhibit internal flows that couple to the surrounding nematic? How does the droplet-to-\ndefect size ratio influence the nature and onset of the various transitions? The distinction between\ninverse MIPS and classical MIPS requires clarification. Could the observed behavior be interpreted\nas conventional MIPS, viewed from the co-moving frame of the active medium? Addressing these\npoints more systematically would provide important physical intuition and broaden their relevance\nacross active and soft matter systems.\nReply: We thank the Reviewer for deeming our results a significant contribution and for the thoughtful\nand constructive feedback. We appreciate the Reviewer’s suggestions that allowed us to enhance the\nphysical interpretation of the mechanisms driving the observed transitions, and we have now amended\nthe manuscript to address these in the revised version.\nIn response to the Reviewer’s suggestion, we now present our results using dimensionless quantities\nto improve physical intuition. Specifically, we introduce the dimensionless activity number\nR2 R2|ζ|\nA = d = d , (2)\nℓ2 κ\nact F\n(cid:112)\nwhere R is the nominal droplet radius, ℓ = κ /|ζ| is the active length scale, ζ is the activity\nd act F\ncoefficient, and κ is the Frank elastic constant. Analogous to the P´eclet number commonly used in\nF\nthe context of MIPS, which compares advective to diffusive transport, the activity number A introduced\nhere quantifies the relative importance of active stresses versus elastic relaxation.\nFurthermore, we rescale energy using the characteristic energy of a +1 nematic defect,\n(cid:18) (cid:19)\nL\ng\nϵ = πκ ln , (3)\ndef F\nℓ\nc\nwhere L is the system size and ℓ is the defect core size. Lengths are normalized by the droplet radius\ng c\nR , and time is rescaled by the elastic relaxation time-scale of the liquid crystal,\nd\nγR2\nτ = d, (4)\nκ\nF\nwithγ denotingtherotationalviscosity. ThesedefinitionsarenowclearlystatedintheMethods section.\nFurthermore, we now characterize in details the role of surface tension in cluster stability (see reply to\nReviewer Comment 1.1 and 3.2), which we rationalize through an active capillary number\nζaR2\nCa = d, (5)\nact σ2\nwith σ denoting the droplet surface tension, and a their compressibility modulus.\nFinally, we have expanded our discussion throughout the manuscript and in our reply to address the\nquestions raised by the Reviewer:\n• The role of droplet deformability and its influence on defect bonding and clustering is discussed\nin our response to Reviewer Comment 3.3.\n• Internal flows within droplets and their coupling to the active nematic are addressed in Reviewer\nComments 3.3 and 3.7.\n• The parameters controlling defect size and their interaction with the suspended droplets are\nclarified in Reviewer Comments 3.5 and 3.6.\n• The connection between our observed activity-driven deformability-induced phase separation (or\nshortlyactive-DIPS[previouslyinverse-MIPS])andclassicalMIPSisdiscussedindetailinReviewer\nComment 3.8.\nWe believe that these additions have significantly improved the manuscript, providing a clearer\nphysical interpretation of the observed transitions and enhancing the generality and relevance of our\nfindings.\nReviewer Comment3.1 — WhiletheauthorsfocusontheinverseMIPSregimeathighactivity\nand moderate packing, I would encourage greater attention to the intermediate percolated network\nphase, which seems capable of temporally pinning topological defects and suppressing active tur-\nbulence. This regime could be particularly promising for achieving structural control in active\nmaterials.\nReply: We thank the reviewer for highlighting the importance of the intermediate percolated network\nphase, which was underrepresented in the initial version of the manuscript. Following this insightful\nsuggestion, we have substantially expanded our discussion of the dynamics leading to gelation and the\nrole of this regime in suppressing active turbulence. In particular, we now detail how defect-mediated\ndropletbondinginthepassivesystemcontributestotheemergenceofapercolatednetwork. Asdescribed\nin the revised manuscript:\nThese bonds arise when droplets share two (or more) of their surrounding topological\ndefects, causing them to adhere (see Fig. 2a,b). Such droplet dimers act as aggregation\ncenters, facilitating the growth of larger clusters through the absorption of isolated droplets\nor the merging with other clusters. In particular, when the packing fraction becomes suf-\nficiently high (Φ ≳ 0.21), the typical distance between droplets becomes small enough for\nelastic interactions to dominate. These short-range attractive forces are strong enough to\ndrive a global rearrangement of the droplets, leading the system to evolve into a stationary\npercolating network (see Fig. 1b and Movie S2), even in the absence of activity.\nWe also emphasize how activity influences this regime:\nAthigherpackingfractions, activityacceleratestheformationofapercolatingnetworkeven\nat intermediate activity levels. As a result, the system undergoes a transition from dynamic\nclusteringtoagelatedstateasthepackingfractionincreases(Fig.1b). Specifically,whenthe\npacking fraction is sufficiently large, active jets mobilize droplets into chain-like structures\nthat eventually intersect, forming a passive network of droplets that encapsulates regions\nof active nematic fluid. Notably, the droplet network confines the active liquid crystal into\ndomains smaller than the overall system size, effectively suppressing unstable infrared active\nmodes, leading the system to settle into a stationary state over long timescales (Movie S2).\nTo further characterize this regime, we now explicitly discuss the radial distribution function g(r):\nIn the gelation phase, Fig. 2c displays a sharp, well-defined first peak corresponding to the\nnearest-neighbor droplet distance, followed by a series of secondary peaks, reflecting short-\nrange positional correlations, which are characteristic of an extended, percolated droplet\nnetwork lacking long-range order.\nReviewer Comment3.2 — I would like to emphasize this point again because I find it to be\none of the most interesting and potentially impactful aspects of the study. The choice to model the\ninclusions as soft, deformable droplets is central to the novelty, but the implications are not fully\nexplored. The only mention I could find in the text is: This explains the crucial role of droplet\ncompressibility in cluster sta-bility [our observations indicate that increasing surface tension–and\nthus reducing compressibility–significantly destabilizes droplet clusters]. This should be elaborated\nquantitatively. What key physical phenomena arise specifically due to deformability? Can the\nauthors change the surface tension or viscosity of the droplets? Would similar clustering or inverse\nMIPS behavior be observed with rigid particles, or is the deformation necessary? What is the\ninternal velocity field? A more in-depth discussion of how softness contributes to defect trapping,\norganization,orphaseseparationwouldsignificantlyimprovethemanuscriptandclarifythebroader\nrelevance of the results. Additionally in figure 1-2-3 the deformations are not visible but in figure\n4 they are. Are the conditions the same?\nReply: We thank the reviewer for emphasizing the central role of droplet deformability in our study.\nWe fully agree that this aspect is both novel and impactful, and we have now significantly expanded\nthe discussion in the revised manuscript to better highlight its implications.\nFollowing the Reviewer’s suggestion, we have systematically varied the surface tension σ of the\ndroplets to explore how compressibility affects cluster formation and their stability. These new results\nare presented in a new Fig. 5 of the amended manuscript and discussed in the revised Clustering\nDynamics section, along with a supplementary movie (Movie S6). Quantitatively, we observe that\nincreasing surface tension—thus reducing compressibility—leads to a shallower pressure gradient within\nclusters (Fig. 5a), while the active force field remains largely unchanged, as it originates from the\nsurrounding active fluid. This mismatch between the diminished internal pressure gradient and the\npersistent external active forcing destabilizes the cluster. We comment on this in the manuscript as\nfollows:\nWe next examine the role of surface tension in the stability of active-DIPS, as illustrated in\nFig. 5. In panel (a), we perform a experiment starting from a clustered state and vary the\nsurface tension σ. As surface tension rises, the cluster undergoes melting, which is accom-\npanied by a sharp decrease in the average hydrodynamic pressure gradient ⟨|∇P|⟩R2/κ .\nd F\nThis behavior is further detailed in the inset, showing the temporal evolution of the pressure\ngradient for two representative values of surface tension: low tension maintains a high pres-\nsure gradient which stabilizes the cluster, whereas high surface tension makes the droplets\nstiffer so that they are not able to absorb external pressure eventually leading to melting.\nFurthermore, wenowpresentaphasediagramatvaryingsurfacetensionandactivity(Fig.5b), revealing\nthat high activity stabilizes active-DIPS, while high surface tension promotes melting and drives the\nsystemtowardanactiveturbulentstate,withthetransitionboundaryscalingapproximatelyasσ ∼ A1/2.\nIn the manuscript, we comment as follows:\nThe phase diagram in Fig. 5b further reveals a competition between activity and surface\ntension; high activity stabilizes active-DIPS, while high surface tension promotes melting\nand drives the system toward an active turbulent state, with the transition boundary scaling\napproximately as σ ∼ A1/2. Our data suggest that this scaling can be rationalized by\nintroducing an active capillary number, defined as\nζaR2\nCa = d, (6)\nact σ2\n√\nwhere a is the bulk compressibility modulus. The observed transition line σ ∼ A ∼\n√\nζ corresponds to constant Ca , confirming that the phase behavior is controlled by\nact\nthe competition between active forcing and interfacial stresses that resist deformation and\nregulate droplet compressibility.\nThis new evidence show that deformability is a key ingredient for phase separation: soft, compressible\ndroplets can absorb active stresses, creating local pressure sinks that stabilize clusters. In contrast, rigid\nor incompressible droplets lack this mechanism, leading to cluster destabilization and eventual melting.\nIn response to the Reviewer’s question about whether similar behavior would be observed with rigid\nparticles: our simulations show that clustering is significantly suppressed at high surface tension, and\npre-formed clusters dissolve when surface tension exceeds a critical threshold. This strongly suggests\nthat deformability is necessary for the active-DIPS behavior we report.\nWe also now include in Fig. R4d of this report a snapshot of the velocity field showing that the\ninternalvelocitywithinclustersisrelativelyuniformandsignificantlylowerthaninthesurroundingactive\nfluid, which exhibits strong vortical patterns typical of active turbulence [The figure has been added to\nthe Supplementary Information where it appears as Fig. S1].\nFinally, regarding the visibility of droplet deformations: the simulation parameters are consistent\nacross the whole manuscript, with the only exception of the new Fig. 5 where surface tension was varied\nascommentedabove. Thereasondeformationsarenotapparentinregimesdifferentthanactive-DIPSis\nthat the surface tension is sufficiently high to resist deformation at the single-droplet level. In Figure 4,\nhowever, the collective imbalance of stresses within large clusters leads to visible shape distortions, even\nunder the same nominal surface tension conditions.\nWe thank the reviewer again for encouraging a deeper exploration of this topic. We believe the new\ndata and expanded discussion significantly strengthen the manuscript and clarify the broader relevance\nof our findings to the design of active materials.\nReviewer Comment3.3 — While the authors focus on inverse MIPS, the percolated network\nregime at intermediate activity and high packing fraction is also highly intriguing. It appears\ncapable of temporarily pinning defects and suppressing turbulence. Do the droplets deform? I\nencourage the authors to elaborate on the structure and dynamics of this regime and its potential\nrelevance to materials with embedded topological control.\nReply: In the gelation regime, droplets do not deform appreciably. The external active and elastic\nstresses are not sufficient to overcome the surface tension, which remains high enough to maintain\ndroplet sphericity. This results in a percolating network of nearly spherical droplets connected via\ndefect-mediated bonds.\nFor a more detailed discussion of the structure and dynamics of this regime–including the role\nof defect bonds, the emergence of a stationary percolated network, and its impact on active flow\nsuppression–we kindly refer the Reviewer to our response to Reviewer Comment 3.1.\nReviewer Comment3.4 — Each droplet introduces a +1 charge, leading to two +½ defects.\nThis generates interactions. Please provide more quantitative context: what is the typical distance\nbetween a droplet and its defects? Are defect locations stable over time, or do they reorganize?\nHow do the defects compensate for the charge of the deformed droplet?\nReply: We thank the referee for raising this point. Each droplet indeed introduces a topological charge\nof +1, which is compensated by the emergence of two −1 defects in the surrounding liquid crystal.\ntypical distance between the droplet and these defects is set by the so-called penetration length, which\narises from the competition between the anchoring strength at the droplet interface and the elastic\nresistance of the liquid crystal. This length scale is given by ℓ = (cid:112)κF,\nwhere κ is the Frank elastic\np W F\nconstant and W is the anchoring constant. At this distance, the system minimizes its free energy, and\nthe defects tend to stabilize in these energetically favorable positions. We have clarified this in the main\ntext, as follows:\nWhen the liquid crystal exhibits either homeotropic or tangential anchoring at the droplet\ninterface, the director field wraps around each droplet in a way that imposes a topological\nchargeof+1. Tosatisfyglobaltopologicalconstraintstheliquidcrystalcompensatesbygen-\neratingtwo−1/2defectsinthesurroundingmedium. Thesedefectsarenotrigidlyattached\n(cid:112)\nto the droplet interface but are energetically stabilized at a distance ℓ = κ /W – the\np F\nso-called penetration length where W denotes the anchoring strength at the droplet surface\n– and, in isolated droplets, they are typically located at antipodal positions (Fig. 1a(i)).\nFor completeness, we note that this description pertains to a two-dimensional configuration; in\nthree-dimensional systems, the analogous structure is the so-called Saturn ring disclination line, a\nclosed loop disclination which encircles the droplet to compensate for the topological charge introduced\nby homeotropic anchoring conditions.\nInthelow-activityregime–suchasduringthequiescentorgelationphase–weobservethatthedefects\nremain close to their equilibrium positions. However, as activity increases and flow fields become more\nprominent, the defects can be advected away from the droplet. In the active turbulent regime, the\nsituation becomes even more dynamic: defects are not only displaced but can also be spontaneously\ncreated and annihilated due to active stresses. These processes occur in pairs of +1 and −1 defects,\n2 2\npreserving the overall topological charge of the system, namely −N , with N the number of droplets\nd d\nin the system.\nReviewer Comment3.5 — InFig. 2a,thedefectsappearsignificantlysmallerthanthedroplets.\nIs this scale difference consistent across activity levels? How defect size or separation change with\ndroplet size?\nReply: We thank the referee for this observation. The apparent size difference between defects and\ndropletsinFig. 2aisindeedconsistentacrossactivitylevels. Thetypicalspatialextensionofatopological\n(cid:112)\ncore is governed by the coherence length, defined as ℓ = κ /A , where κ is the elastic constant\nc F 0 F\nand A is the bulk constant of the liquid crystal. This length scale characterizes how rapidly the liquid\n0\ncrystalline order decays near a defect core, effectively setting the defect’s size.\nImportantly, the size of topological defects is determined solely by the intrinsic properties of the\nliquid crystal and does not depend on the size of the droplet that induces the defect. Since we do\nnot vary either κ or A across different simulations, the defect size remains invariant with respect to\nF 0\nactivity level and droplet size.\nIn our simulations, we set ℓ ≈ 0.5, which corresponds to a defect core size of approximately 1–2\nc\ngrid points in diameter. This choice is customary in numerical studies, as the coherence length is often\nused as an ultraviolet cutoff beyond which the hydrodynamic description ceases to be valid. By ensuring\nthat our discretization respects this cutoff, we avoid sampling unphysical length scales. We have added\nthe following sentence in the Method section while discussing the simulation parameters:\n(cid:112)\nWith the chosen parameters, the liquid crystal coherence length is ℓ = κ /A ≈ 0.5.\nc F 0\nThis length scale characterizes how rapidly the orientational order decays near a defect core,\neffectively setting the defect size to approximately one grid point.\nReviewer Comment3.6 — What is the flow in the vesicles? Do the droplets respond purely\nto external flows, or do they exhibit internal velocity or field gradients? This could influence how\nthey deform or interact with defects. I think the paper would benefit from showing how a single\nsoft particle behaves at different activity levels, specifically in terms of its deformation, the flow\nfields inside and outside the particle, and its mean squared displacement.\nReply: We agree with the referee that the internal flow within the droplets is an important aspect\nto consider. In our model, droplets contribute to the overall velocity field through an interfacial stress\nterm, given in Eq. (5) of the manuscript:\nσϕ = (cid:88) N d (cid:20)(cid:18) f −ϕ δF (cid:19) δ − ∂f ∂ ϕ (cid:21) ,\nαβ i δϕ αβ ∂(∂ ϕ ) α i\ni β i\ni=1\nwhere the deviatoric contribution, proportional to k ∇ϕ∇ϕ, generates a flow that promotes relaxation\nϕ\ntoward the equilibrium (spherical) droplet shape. However, in passive systems, the magnitude of these\nflowsistypicallyontheorderofO(10−4)l.u.,whichissignificantlysmallerthantheactiveflowsobserved\nin our simulations, which can reach O(10−2) l.u. or more. As a result, the passive interfacial flows\nare effectively masked by the dominant active dynamics, making their direct visualization and analysis\nchallenging. This means that for the choice of the model parameters considered in this study, the\ninterfacial stress exerted by the droplet is too weak to affect the droplets and defects dynamics.\nRegarding the referee’s suggestion to include a detailed analysis of an isolated droplet at varying\nactivity levels, we note that this scenario has been comprehensively addressed in a recent study by\nSingh and Chaudhuri [64], where the deformation, internal flow, and dynamics of single droplets were\nsystematically explored as a function of activity, droplet size, and surface tension. Our observations are\nin agreement with the findings reported in that work. For this reason, and to avoid redundancy, we have\nchosen not to include a dedicated section on isolated droplets in the present manuscript.\nNonetheless, to provide a broader picture of the flow behavior across different activity regimes, we\nnow include in Fig. R4 a gallery of velocity fields corresponding to the various phases discussed in the\nmanuscript. The figure and a short comment have been added to the Supplementary Information where\nit appears as Fig. S1.\nReviewer Comment3.7 — The inverse MIPS terminology is compelling but potentially con-\nfusing. In classical MIPS, phase separation is driven by self-propelled particle motility. Here, the\ndroplets are passive, and it is the surrounding fluid that is active. The authors should, in my per-\nsonal opinion, name the effect differently to distinguish between the options: a clustering regime\nviewed from a passive-particle perspective or, as the name implies, the breaking of MIPS at high\nactivity. And incorporate the deformability of the particles in the terminology.\nReply: We thank the Reviewer for their comment regarding the use of the term “inverse MIPS.” We\nagree that this terminology may lead to confusion. As the Reviewer points out, classical MIPS refers to\nmotility-induced phase separation in dry active systems, where self-propelled particles slow down at high\ndensity and cluster as a result. In our case, however, the droplets are passive and suspended in an active\nfluid, and the clustering behavior arises from the interplay between activity and droplet deformability,\nrather than from self-propulsion.\nFollowing the suggestions of both Reviewer 1 and Reviewer 3 (see also Reviewer Comment 1.3), we\nhaverevisedtheterminologythroughoutthemanuscript. WenowrefertothisregimeasActivity-Driven\nDeformability-Induced Phase Separation (abbreviated as active-DIPS). As the Reviewer suggests, this\nnew name highlights the key physical ingredients of the phenomenon: the non-equilibrium (activity-\ndriven) nature of the system and the crucial role of droplet deformability in enabling phase separation.\nReviewer Comment3.8 — It would be helpful to include average velocity fields for each phase.\nVisualize how droplets interact with active flows and whether there are correlations between flow\nstructures and droplet clusters. Maybe in Fig. 1 or Fig. 3.\nReply: We thank the Reviewer for their suggestion. We have now included representative velocity field\nconfigurations for each of the phases discussed in the manuscript. These visualizations are provided in\nFig. R4 of this report and have been added to the Supplementary Information as Fig. S1.\nReviewer Comment3.9 — Figure 1 would benefit from clearer labeling and a more intuitive\nlayout. Presenting representative snapshots for all phases in a single row or panel with labels, with\nconsistent activity and packing fraction annotations.\nReply: We thank the Reviewer for the suggestion to improve the clarity and layout of Figure 1.\nFollowing this advice, we have revised the figure to present snapshots of all identified phases in a single\nrow. Eachconfigurationisnowclearlylabeledwithannotationsforactivitynumberandpackingfraction.\nReviewer Comment3.10 — The model is 2D, but many soft biological systems, such as\nvesicles and LLPS condensates, are inherently 3D, and experimental nematics fields with +-1/2\ndefects are 2D. The dimensional mismatch must be discussed. What real systems does this model\nbest represent? For example, might LLPS droplets be in 2D planar confinement?\nReply: We thank the Reviewer for raising the important issue of dimensionality. Several experimental\nsystems relevant to biological and soft matter physics indeed operate under quasi-two-dimensional con-\nfinement. For example, Wu and Libchaber [29] studied passive beads suspended in a freely suspended\nsoap film containing motile E. coli, effectively creating a quasi-2D bacterial bath. Similarly, V´elez-Cer´on\netal. [30]embeddedhydrogelpillars(fixedtothesubstrate)withinaquasi-2Dactivenematiclayercom-\nposed of microtubules and kinesin motors. These systems demonstrate how soft embedded structures\ncan be used to probe active stresses and defect dynamics in confined geometries. One could envision ex-\ntending such setups to include deformable inclusions, thereby realizing a quasi-2D active-passive system\nakin to the one modeled in our study.\nThese experimental realizations illustrate that quasi-2D configurations are not only feasible but also\nhighly relevant for investigating active-passive interactions. Our model is conceptually aligned with such\nsystems and may represent scenarios such as LLPS droplets under planar confinement in microfluidic\nenvironments, where soft inclusions interact with active flows. We have now added a brief discussion in\nthe conclusions, which reads as follows:\nOur model conceptually aligns with quasi-2D experimental systems such as active nematic\nlayersconfinedatawater-oilinterfaceorsoapfilms, wheresoftinclusionslikeLLPSdroplets\nor hydrogel particles can be embedded and interact with active flows. These setups offer a\npromising route to experimentally realize and test the mechanisms uncovered in this study.\nReviewer Comment3.11 — Instead of plotting full MSD curves for every regime (Fig. 3a), a\nmore compact representation showing the short-time and long-time exponents vs activity would be\nmore informative. As done for the late times, the MSDs could be collapsed on the same figure for\neasy comparison with boundaries between 1-2. This would make it easier to pinpoint transitions\nbetween diffusive and turbulent regimes.\nReply: We thank the Reviewer for the suggestion. We have explored several layout options for Fig. 3\nwith the goal of condensing the information as much as possible. After careful consideration, we\nfound that presenting representative MSD curves for each regime offers the clearest and most intuitive\nvisualization. Overlapping all curves in a single panel, while compact, would result in significant visual\nclutter.\nRegarding the suggestion to plot the short-time MSD exponent, we now include a plot of the short-\ntime exponent as a function of activity in the Supplementary Information where it appears as Fig. S4.\nWe hope the Reviewer finds this a reasonable compromise that maintains clarity while addressing the\ncore of their suggestion.\nReviewer 3\nThe revised manuscript is much clearer, and I am happy with the changes and recommend for\npublication.\nReply: We thank the Reviewer for acknowledging the quality of the revised manuscript and the\nimprovements made relative to the initial submission.\nMy only remaining difficulty is:\nReviewer Comment3.1 — Fig. 3. as presented, the MSD curves all appear to lie within a\nsimilar range, and although the panels are divided into two activity regions, there is no clear visible\ndependence on the morphological regime. It is therefore not obvious what additional insight is\ngained by presenting so many separate panels. I see this primarily as an issue of presentation; the\nunderlying analysis does not need to be changed.\nReply: We appreciate the Reviewer’s observation regarding the graphical presentation of the MSD\ncurves. We have previously explored the possibility of collapsing these data into a single panel; however,\nwe found that doing so makes it nearly impossible to distinguish the five distinct physical behaviors\ndue to the overcrowding. By maintaining separate panels, we are able to highlight subtle but physically\nsignificant features that would otherwise be obscured, such as the varying sharpness and timing of the\ntransitionstolate-timebehavior. Thispresentationallowsforaclearercomparisonbetweenregimes—for\nexample, distinguishingthesub-ballisticexponentsinclusteringfromtheballisticnatureofactiveturbu-\nlence—and specifically illustrates the continuous transition unique to the active-DIPS case. We believe\nthat this graphical approach is the most effective way to ensure the reader can accurately identify these\nsignatures across the different morphological states.\nReviewer Comment3.2 — Error estimates. For both Fig. 3 and Fig. 5a, it is unclear whether\nthe data represent single simulations for each condition or averages over multiple realizations. The\ncurves fluctuate and are not very smooth, which is reasonable for simulation data. However, I do\nnot see any error bars or other indications of variability, and the captions do not discuss errors\nor averaging. If multiple simulations were performed, the authors should indicate how many, how\nthe data were averaged, and include an appropriate measure of variability ( error bars or shaded\nconfidence intervals). If only single realizations are shown, this should be stated explicitly in the\ncaptions, along with a brief justification that this level of statistics is consistent with the data.\nReply: We thank the Reviewer for highlighting this point. Regarding Figure 3, the caption states\nthat the long-time MSD exponents are obtained by averaging over n = 5 independent simulations.\nFurthermore, the standard error of the mean for the long-time MSD exponents is reported in panel f.\nPanels g and h display the average of the fluid and droplet kinetic energy at late simulation times once\nthe non-equilibrium steady state was reached. The calculated standard eeror of the mean is consistently\nsmaller than the marker size, meaning the error bars fall entirely within the data points and are not\nvisually distinguishable.\nAlso in Figure 5a, the reported data represent averages calculated over long simulation times once\nthe non-equilibrium steady state was achieved. We have now included error bars in panel a (denoting\nstandarddeviation)torepresentthetemporalvariabilityofthesemeasurements. Thisadditionquantifies\nthe fluctuations in the active-DIPS regime, which are inherently more pronounced than the stable\nbehavior observed in the melted state."}]

Decision letter