Geometric control of motility-induced phase separation

Abstract

Curvature fundamentally alters the collective properties of soft, active, and biological materials. Here we study motility-induced phase separation (MIPS), a canonical non-equilibrium transition, and demonstrate that even weak and slowly varying curvature provides robust geometric control over the dense MIPS phase. This includes dictating both the location and morphology of the MIPS cluster, even in regimes where curvature has minimal effect on the overall phase boundaries. Focusing on active Brownian particles confined to the surface of a torus, we show that varying the aspect ratio drives a structural transition of the dense cluster from a disk localized at the outer equator to a band wrapping the minor circumference. We then discuss how the curved geometry provides a platform for comparing different theoretical frameworks for the MIPS phase: by analyzing the geometries of the cluster boundaries, we show that the dense phase shape is more consistent with a boundary-length-minimizing, thermodynamic picture than with the simplest kinetic picture in the large-particle-number limit. Our results establish curved space not only as a tool to shape and guide non-equilibrium dynamics, but as a uniquely sensitive arena for probing the fundamental mechanisms of active matter.

---

1 Introduction

Curvature fundamentally alters the properties of materials and is an inherent feature of many biological environments. In equilibrium soft matter, curvature dictates the geometry and dynamics of phase separation.^1–9 In living systems, this geometric influence scales from the localized distribution of molecules in cell membranes¹⁰ up to the curvotactic guidance of migrating cell collectives.^11–14 Motivated by the link between geometry and biological function, a growing body of work has sought to understand how synthetic active materials respond to curved environments. This has revealed that introducing spatial curvature to active systems yields novel phenomena, including modified flocking behavior,^15–17 the trapping of semi-flexible filaments,¹⁸ coordinated continuous rotation,^19–22 and curvature-dependent rigidity in epithelial models.^23,24

Motility-induced phase separation (MIPS), in which self-propelled particles phase separate into coexisting dense and gaseous phases, is one of the most well-studied phenomena to emerge in scalar active matter systems.^25–27 This phenomenon is observed in computational studies of minimal active models,^28–31 in experimental studies of self-propelled colloids,^32,33 and in collective motion of biological systems such as bacteria^34,35 and fire ants.³⁶ The introduction of curvature to these models, whether through boundary confinement³⁷ or continuous surface geometry,^38,39 fundamentally alters MIPS phase behavior. Simulations of active Brownian particles (ABPs) have shown that surface curvature can alter the non-equilibrium phase diagram, with large curvature scales able to substantially alter the parameter regime in which phase separation occurs.^38,39 Schönhöfer and Glotzer, in particular, showed that positive curvature can promote clustering while negative curvature inhibits it, with the strength of this effect depending on the relative magnitude of the Gaussian curvature K to the size scale of the individual particles σ.³⁹ At lower curvature scales, for which the quantity |K|σ² ≪ 1, the phase boundaries approach the Euclidean space limit, as would be expected.

However, whether slowly varying curvature in this regime can exert control over the morphology of the dense phase remains an open question. Here we show that it indeed has important consequences, even in parameter regimes where |K|σ² is small and its effect on the overall phase boundary is minimal. By studying MIPS on the surface of a torus, which contains regions of positive and negative Gaussian curvature, we demonstrate that curvature provides robust geometric control over both the shape and location of the dense cluster. Reminiscent of how equilibrium phases seek optimal perimeter-to-stress balances on curved surfaces,⁹ at a low torus aspect ratio ξ, the dense active phase forms a disk that preferentially sits near the outer equator in the region of positive curvature. At high ξ the cluster instead forms a band that wraps the minor circumference of the torus.

Our curved space simulations also provide a platform for breaking the geometric degeneracy between different mechanisms that can drive the formation and shape of the dense cluster. Many features of MIPS map directly onto expectations for equilibrium interfaces,^25,40,41 and theoretical models suggest that cluster interfaces follow surface-tension-driven area-minimizing principles.^42,43 An alternative interpretation treats the cluster as a kinetic structure that is shaped by the balance of particle deposition and escape.^29,44–46 In Euclidean space, both of these interpretations lead to the averaged dense phase cluster shape being a disk or a band, depending on the area fraction covered by the dense phase and the aspect ratio of the system.³¹ However, in regions of varying curvature, these interpretations predict distinguishable shapes for the cluster boundaries. By leveraging the geometry of the torus, we directly compare the dense phase boundaries against the structures predicted by thermodynamic perimeter-minimization and kinetic frameworks. Our results highlight how curvature can act as a powerful tool for revealing the physical principles that govern active processes.

2 Methods

We investigate active Brownian particles (ABPs) with soft repulsive interactions that move and interact entirely within a curved surface. The motion of the particles is governed by the coupled overdamped Langevin equationswhere r_i is the position of particle i, θ_i is the orientation of its heading measured relative to a reference axis, and dots denote time derivatives. The translational motion of the particles has contributions from the self-propulsion v_i and the net force F_i on a particle due to neighbor interactions.

The self-propulsion is given by v_i = v₀(cos θ_iê₁ + sin θ_iê₂), where v₀ is the self-propulsion speed and ê₁ and ê₂ are orthogonal coordinates in the local tangent plane at a given point on the surface. We set v₀ and the time-scale τ such that a free particle self-propels one particle diameter σ in time τ, giving v₀ = σ/τ, and we use dt = 0.01τ as our simulation timestep. Particle orientations are subject to Gaussian white noise obeying 〈η^r_i〉 = 0 and 〈η^r_i(t)η^r_j(t′)〉 = 2D_rδ_ijδ(t − t′). Here D_r is the rotational diffusion constant, which has a corresponding timescale τ_r = 1/D_r.

Particles interact with each other to give a net force on a particle due to its neighbors. Here we use a soft harmonic repulsion

where r̂_ij is a unit vector tangent to the geodesic connecting points i and j at point i, and k is the stiffness. We set k = 30 and the mobility μ = 1 such that the elastic relaxation time scale τ_e = (μk)⁻¹ is small.

Due to challenges posed by curvature, many studies of particulate systems on curved surfaces calculate interparticle forces and particle motion using a Euclidean metric with additional forces imposed to ensure the particles remain within the surface.^{8,9,15,18,22,39} Here, we instead use a fully curved-space approach where all particle motion and interactions occur within the surface using geodesic distances to calculate forces and in-surface vector transport for particle translation, allowing us to truly probe the effects of inherent curvature in this system. We conduct all simulations using curvedSpaceSim,⁴⁷ which performs these calculations on arbitrary connected non-self-intersecting embedded surfaces represented by triangulated mesh objects. The equations of motion are integrated using a geodesic-based Euler–Maruyama scheme, with rotational noise applied in the tangent plane of the particles.

In this system there are three relevant length scales—the particle diameter σ, the persistence length l_p = v₀/D_r = 100σ, and a length scale l set by the surface. In Euclidean space, the phase behavior of this system can be characterized in terms of the rotational Péclet number Pe_r = l_p/σ, which characterizes how far a particle moves before it changes orientation, and the total area fraction ϕ of particles on the surface.^29,30 Schönhöfer and Glotzer showed that the effect of curvature on the phase behavior can be characterized through the quantity Kσ², where K is the Gaussian curvature, while the collective properties of the system are independent of l_p for l_p > 50σ.³⁹ Here we highlight the effect of curvature on properties of the dense phase in regimes where it has minimal effect on the overall phase separation. We therefore set D_r = 0.01τ⁻¹ to give Pe_r = 100, and investigate ϕ = 0.4, which sits well within the region of phase space where MIPS is expected to occur.^29,30,39 We treat particles as circular disks with the diameter σ determined such that a system of N particles obeys Nπ(σ/2)² = ϕA, where A is the total surface area of the system. The approximation that the area of a particle is its Euclidean area is reasonable for sufficiently large N, as used here. We verified that the phase behavior of our system is consistent with prior expectations by investigating interaction strengths between k = 10 and k = 100 and area fractions between 0.2 and 0.8 for a torus with aspect ratio ξ = 1.5.

Having established this baseline, we focus on how the specific geometric parameters of the torus influence MIPS. Coordinates in the surface are parameterized by

where the major axis R is the distance from the center of the torus to the center of the tube, and the minor axis r is the radius of the tube. The toroidal angle θ is measured around the major circumference, while angles around the minor circumference are measured using the poloidal angle ψ (see Fig. 1(a)). The torus aspect ratio is ξ = R/r. We set the unit of length by fixing the surface area of the torus A = 4π²rR to A = 4π²l², and set l = 1. The major and minor axis lengths of the torus are therefore R = ξ^1/2l and r = ξ^−1/2l respectively. We vary the aspect ratio in the range 1.1 ≤ ξ ≤ 3 and study the resulting dense phase behavior.

Fig. 1 (a) Snapshot of a steady-state configuration for a torus with aspect ratio ξ = 1.5 and N = 20 000 particles. Here the dense phase forms a disk. (b) Steady-state snapshot for ξ = 3.0 and N = 5000. Here the dense phase forms a band that wraps around the minor circumference of the torus. (c) and (d) Corresponding time-averaged density profiles for particles in the dense phase. Cluster locations are centered in the θ direction before averaging. Dashed lines show the contour enclosing the average dense phase area fraction . The outer equator has ψ = 0.

To ensure we are probing regimes where the curvature has minimal effect on the phase behavior, we need the curvature scale to be sufficiently small relative to σ. On the torus, we have

which attains a maximum positive value at the outer equator (ψ = 0) and a minimum negative value at the inner equator (ψ = π). Ref. 39 identified that the low-curvature regime is bounded by |K|σ² ≲ 0.003, which sets a ξ-dependent lower bound on our choice of N. At the outer equator, this varies from N ≈ 3500 for ξ = 1.1 up to N ≈ 5000 for ξ = 3.0. At the inner equator, the magnitude of the curvature is larger, meaning we would require N ≈ 74 000 for ξ = 1.1, which comes down to N ≈ 10 000 for ξ = 3. Given the computational intensity of our simulations, here we use N = 5000 for the majority of our simulations exploring the behavior of the dense phase as a function of ξ. To verify that our results are robust at larger sizes, we have performed simulations with N = 20 000 for select values of ξ that span the different regimes we observe.

In all simulations, we initialize the system in a random configuration and allow the dynamics to reach steady state. We then record snapshots spaced at intervals of at least 200τ for a duration of at least 40 000τ. Full details of the data collected for each ξ and N are listed in Table 1.

Table 1 Summary of data used in analysis of different state points on the torus. Here N_s is the number of independent samples, N_f is the minimum number of frames recorded from a given sample after reaching steady state, and τ_f is the time interval between recorded frames

N	ξ	Ns	Nf	τf/τ
a 10 samples have 200 frames, 10 samples have 450 frames.
5000	1.5	20	200/450^a	200
	1.6, 1.7, 2.1	10	450	200
	1.2, 2.0, 2.2	5	450	200
	3.0	10	540	400
	1.1, 1.4, 1.8, 1.9, 2.4, 2.6, 2.8	1	250	200
20 000	1.2	1	450	300
	1.5, 2.0, 2.2	10	200	300
	1.8, 2.4	5	300	300

---

3 Results

3.1 Morphological transition and localization on the torus

We begin by investigating general features of the dense phase as the aspect ratio ξ of the torus is varied. In the low ξ regime, the dense phase is typically disk-shaped, while at high ξ the dense cluster instead spans the torus minor circumference to form a band. Typical snapshots in each of these regimes are shown in Fig. 1(a) and (b) for ξ = 1.5 and 3 respectively. To investigate the typical cluster shape, we plot the averaged density profiles for these ξ values in Fig. 1(c) and (d). These are measured after centering the clusters in each snapshot in the θ coordinate direction, and highlight the distinct disk and band shapes in the low- and high-ξ regimes.

To characterize features of the dense phase structure, we calculate its curved-space area and perimeter. Details of the calculations are given in the Appendix. Fig. 2(a) shows the resulting distributions of the area fraction ϕ_d of the surface covered by dense phase as ξ is varied. We observe that the overall shapes of the distributions are approximately Gaussian for all ξ values, and that the mean ϕ_d across all of the N = 5000 data is . There is slight shift in the distributions as ξ is varied, which we further investigate by plotting the mean area fraction of the dense phase for each ξ in Fig. 2(b). This indicates that ϕ_d initially decreases as ξ increases, before increasing again around ξ = 2 to reach a plateau at large ξ. This total change in the area fraction is small compared to the distribution width, and represents a change of ∼3% of the total number of particles in the dense phase. The effect also becomes even smaller at larger N, suggesting the slight variation could be due to the effect of curvature being greater at low N near the inner equator, causing a slight reduction in cluster formation at intermediate ξ. Overall this result shows that the dense phase cluster area is largely independent of the torus aspect ratio at state points well within the phase separated regime.

Fig. 2 (a) Probability density distribution of the area fraction ϕ_d of the torus covered by the dense phase for different aspect ratios ξ. The color scale is shown in panel (c). (b) Variation in the mean area fraction with aspect ratio for different numbers of particles N. (c) Probability density distribution of the perimeter P of the dense phase, scaled by the perimeter of an ideal band P_b. (d) Peak value of the perimeter distribution for different ξ and N. (e) Fraction of snapshots f_b that form a band as ξ is varied. The solid lines are fits to for each N, which give (ξ_c, h) = (2.02, 0.3) for N = 5000 and (2.07, 0.2) for N = 20 000. The ξ_c values are indicated by the dotted lines while the turquoise and purple dashes respectively show predictions for ξ_c for disks of constant geodesic curvature (ξ_c = 1.68) and radius (ξ_c = 2.10). (f) Probability density distribution of the ψ coordinate center of the dense phase for ξ = 1.5 for different N. The outer equator of the torus has ψ = 0. Distributions in panels (a), (c), and (f) are histograms constructed with a Gaussian smoothing kernel. The color scale in panel (b) also applies to panels (d)–(f).

By contrast, the cluster perimeters, calculated as the total geodesic length of the cluster boundaries, show a much stronger ξ dependence. The distributions, shown in Fig. 2(c) with the perimeter P scaled by the perimeter of an idealized banded configuration (P_b = 2C_r, where is the minor circumference of the torus), have a clear peak and a long tail. At larger ξ, where the dense phase forms a band, the peaks of the distributions appear to converge to a value of 1.4P_b, indicating that in this regime the typical measured perimeter length relative to the minor circumference is constant. Examining the position of the distribution peaks P_Max in Fig. 2(d), we find that the typical perimeter remains constant in the low-ξ disk regime. In the banded regime, P_Max decreases with ξ such that the ratio P_Max/P_b plateaus.

We characterize the transition from a disk-shaped dense phase to a band by determining the fraction of banded snapshots f_b at a given ξ, details of which are given in the Appendix. Fig. 2(e) shows how f_b increases from zero at small ξ, and saturates at one for larger ξ. By fitting a curve of the form

we estimate that the critical aspect ratio for the transition is ξ_c = 2.02 and the width of the transition is h = 0.3 for N = 5000. At larger N, our select data points show that ξ_c shifts to a larger value (ξ_c = 2.07, h = 0.2 for N = 20 000), indicating that our estimate for ξ_c should be taken as a lower bound on the true transition value.

Returning to the cluster density profile for low ξ (Fig. 1(c)), we observe that the resulting averaged cluster shape suggests that the cluster center is typically close to the outer equator of the torus, in the region of positive Gaussian curvature. This suggests that the cluster does not uniformly explore the entire surface, as it would in Euclidean space. To investigate this further, we plot the distribution of the cluster ψ coordinate center ψ_c in Fig. 2(f). This confirms that the cluster centers are peaked around the outer equator, with very few clusters centered in the inner region of the torus (|ψ| > π/2). At larger N, the distribution becomes narrower, highlighting that the localization of the dense phase to the outer equator is robust to larger N and is not merely a consequence of the potential reduced clustering in the higher curvature regime at the inner equator for N = 5000.

3.2 Thermodynamic and kinetic driving mechanisms

What causes the cluster to preferentially sit at the outer equator, and its eventual transition to a band? To answer this, we must examine the physical mechanisms governing the dense phase interface. In Euclidean space, the two dominant theoretical frameworks for MIPS yield identical predictions for the macroscopic cluster shape. Thermodynamic models, based on capillary tension and interfacial dynamics, suggest that the cluster minimizes its boundary length for a given area, analogous to an equilibrium fluid droplet.^41,42 Alternatively, kinetic theories posit that the cluster is shaped by a local balance of particle deposition and escape rates, resulting in an interface that on average is equidistant from a central core.^29,46 In flat space, both mechanisms produce a circular disk or a straight-edged band.

Curved space breaks this geometric degeneracy. The thermodynamic framework of an emergent surface tension maps onto the mathematical isoperimetric problem; this is in general very challenging to solve, but the allowed solutions for a broad class of surfaces of revolution, including the torus, are known. On the torus, isoperimetric domains are bounded by curves of constant geodesic curvature κ_g,^48,49 which quantifies the curvature of a curve γ relative to the underlying surface—geodesics have κ_g = 0, while nonzero κ_g gives information about how γ curves within the surface.⁵⁰ These domains are bands bounded by meridians (lines of constant θ, which are geodesics and therefore have κ_g = 0), and disks with constant κ_g that are symmetric with respect to the equators.⁴⁹ Due to the Gaussian curvature profile of the torus, a constant κ_g disk centered at the outer equator will always have a shorter perimeter than one enclosing the same area at the inner equator. This directly aligns with our observation in Fig. 2(f) that the low-ξ dense phase preferentially sits at the outer equator.

Conversely, a naive application of the kinetic framework predicts that a stable cluster maintains an equal geodesic distance from its center. In curved space, this corresponds to a curve of constant geodesic radius, r_g measured relative to a central point or line. In regions of varying curvature, curves of constant κ_g and constant r_g take on distinct morphological profiles. By analyzing the dense phase on the torus, we can therefore directly probe which geometric bound the MIPS cluster actually obeys.

To first determine the expected disk-to-band transition point, ξ_c, under the isoperimetric assumption, we investigate how the perimeter P of a constant κ_g disk centered at the outer equator varies with ξ. Details of how the constant κ_g disks are constructed are given in the Appendix. The results, plotted in Fig. 3(a), show that for small ϕ_d the ratio P/P_b is always less than one, indicating that the disk solution is always optimal. At larger ϕ_d, there is a critical aspect ratio above which the enclosed region can reduce its total perimeter by transitioning from a disk into a band.

Fig. 3 (a) Variation in the perimeter P of a constant geodesic curvature κ_g disk centered at the outer equator, scaled by the perimeter P_b of an ideal banded solution as the aspect ratio ξ is varied. Curves correspond to different values of ϕ_d, the area fraction covered by the disk. The solid turquoise line indicates , the average dense phase area fraction in the simulations, and the corresponding vertical dashed line shows the predicted transition value ξ_c = 1.68, above which a banded solution has a lower perimeter than a disk. The torus inset shows the shape of the disk solution at the predicted ξ_c. (b) Variation in the minimized perimeter P_ψ0 of a disk with fixed area fraction ϕ_d = 0.1 as the center position ψ₀ of the disk is varied for different values of ξ. Perimeter values are scaled by the perimeter P₀ of a disk centered at ψ₀ = 0. Solid markers show solutions obtained via numerical optimization while open markers indicate values for which the solution is a curve of constant κ_g. Note that for ξ = 1.1, the numerical approach fails to find a valid solution near the inner equator. At the inner equator, the disk solution self-intersects slightly, indicating that the solution would instead be a cylinder spanning the inner equator. Inset tori show solutions for ξ = 2.1 for .

Using from our simulation data leads to a predicted value of ξ_c = 1.68 above which the banded solution becomes optimal. This value is notably smaller than ξ_c = 2.0–2.1 that we observed in our simulation data. One reason for this apparent disagreement could be that even above the predicted ξ_c the constant κ_g disks are still local solutions to the isoperimetric problem, even though the banded solution is the global optimum. Close to the transition, the boundary of the constant κ_g disk is mostly in the regions with positive or zero Gaussian curvature, as shown in the inset in Fig. 3(a). As such, there is a finite energy barrier for the system to transition between the disk and band-shaped structures. As ξ increases, the boundaries of the disk become closer to the inner equator, meaning that it is easier for fluctuations to allow the system to find the globally optimal banded state.

Although the optimal solution for low ξ is a disk centered at the outer equator, our simulation data indicates that the cluster center does deviate from this line (see Fig. 2(f)). We therefore explore how moving the disk away from the equator changes the resulting boundary perimeter. Given that the isoperimetric solutions are centered at the outer equator, these deviations give an indication of the energy penalty of the center straying from the outer equator. Any disk centered away from the equator is, at best, a local solution to the isoperimetric problem.

To find these local solutions as the disk center location is varied, we use the shape optimization software Morpho⁵¹ (see the Appendix for details). We consider a disk with ϕ_d = 0.1 both for technical reasons in the numerical implementation and so that the disk is a valid solution close to the inner equator without self-intersecting. Fig. 3(b) shows the disk perimeters P_ψ0 as the center location ψ₀ is varied, plotted relative to the perimeter P₀ of the disk at the outer equator (ψ₀ = 0). Note that the P₀ values at the outer equator are almost identical across ξ values, with less than 2% variation between the minimum and maximum value. We observe that, as expected, the perimeter of the disk increases as the cluster center approaches the inner equator, with the ξ-dependent maximum perimeter being 1.3–1.65P₀. Most of the increase in perimeter occurs in the region ψ₀ > π/2, where the Gaussian curvature of the surface is negative. Up to ψ₀ = π/2, the perimeter increase is only around 10%, meaning that deviations of the cluster center within the region of positive Gaussian curvature (ψ < π/2) give only a small penalty on the minimum possible length of the boundary compared to the optimal solution.

To determine how the observed dense phase shape compares with theoretical expectations for disk-shaped clusters, we return to the averaged density profiles for ξ = 1.5, shown in Fig. 1(c) for N = 20 000. From these profiles, we extract contours of constant ϕ_d, and compare their shape to area-matched regions bounded by curves with constant κ_g and r_g centered at the outer equator. Fig. 4(a) shows such a comparison for a contour from the N = 20 000 data enclosing an area fraction . Qualitatively, the cluster contour is slightly closer in shape to the constant κ_g solution than the constant r_g one, although its shape sits between the two.

Fig. 4 (a) and (b) Comparison between a contour of constant density, taken from the density plot in Fig. 1(c) and area matched curves of constant geodesic curvature κ_g and geodesic radius r_g centered at the outer equator, shown (a) on the torus surface, and (b) in torus coordinate space. The area of the disks is taken as the average area fraction of the dense phase cluster . The color scale for the constant κ_g and r_g curves indicated in panel (b) applies to all figure panels. (c) Variation in the perimeter P of these different curves as the area fraction ϕ_d is varied. Results are shown for contours taken from the N = 5000 samples (black) and N = 20 000 samples (gray). The dashed vertical line indicates . (d) Histogram of the perimeter deviation w between individual disk-shaped clusters with centers within ψ_Max = π/8 of the outer equator and area-matched curves of constant κ_g and r_g centered at the outer equator. Results are shown scaled by the particle diameter σ. Solid lines are for ξ = 1.5, N = 5000, dashed lines are for ξ = 3.0, N = 5000 and lighter shades are for ξ = 1.5, N = 20 000. (inset) Mean deviation w̄/σ for the ξ = 1.5 distributions as ψ_Max is varied. Solid markers are for N = 5000, open markers are for N = 20 000. Error bars, which are typically smaller than the marker size, indicate the standard error across independent trajectories.

To quantify the differences between these shapes, in Fig. 4(b) we plot the perimeters of each for a range of ϕ_d values. At smaller ϕ_d, the constant κ_g and r_g solutions are almost indistinguishable, which is to be expected for small regions close to the outer equator where the Gaussian curvature is approximately constant. As ϕ_d increases and the constant κ_g and r_g solutions deviate, the perimeter of the constant r_g solution becomes increasingly larger than the constant κ_g solution. Interestingly, we observe that the N = 5000 contour perimeters track the constant r_g perimeters very closely, while for N = 20 000 the perimeters instead track the constant κ_g perimeters for . This therefore indicates that as the system approaches the continuum limit, the dense phase cluster shape appears consistent with the boundary-length-minimizing isoperimetric solution.

In addition to comparing the shapes of the averaged clusters, we can also compare the typical shapes of individual clusters to the two different types of geodesic disk. We do this by measuring the averaged “width” w of cluster interfaces, which quantifies the deviation between the cluster shape and area-matched disks with constant κ_g and r_g. We calculate w using⁴¹

where s is a measure of arc length around the curve, δh(s) is the distance between the reference curve and the interface, and P_d is the total perimeter of the reference curve. Further details of the calculation are given in the Appendix.

Here, we again focus on ξ = 1.5, and consider only disk-shaped clusters that have their centers within an angle ψ_Max = π/8 of the outer equator. Fig. 4(d) shows the resulting histograms for w, measured in units of σ. For N = 5000, we observe that the constant r_g histogram peaks at a smaller value than the constant κ_g histogram, consistent with the observation that the averaged cluster data at this N more closely tracks the constant r_g data. This is also reflected in the inset to Fig. 4(d), which shows how the mean width w̄ relative to the constant r_g disks is consistently smaller as ψ_Max is varied. For N = 20 000, we instead observe that the two histograms are very similar. This is also refelcted in the w̄ values, for which the constant κ_g values are slightly smaller but not within statistical significance. Note that σ for N = 5000 is twice that for N = 20 000, meaning that the actual w̄ values are smaller for the higher N in both cases. These results are consistent with the averaged cluster shape analysis, again showing that the shape of the lower-N data is closer to that of a constant r_g disk, but at higher N the shape becomes closer to the constant κ_g disk.

In the banded regime, we find that for ξ = 3 with N = 5000 particles, the cluster shape is much more consistent with the isoperimetric solution. Using the averaged cluster density profile in Fig. 1(d), we calculate the total perimeters of contours enclosing fixed ϕ_d. We observe that the contour perimeters remain roughly constant as ϕ_d is varied and that the contour perimeter P_c enclosing area fraction is 1.02P_b, where P_b is the perimeter of the constant κ_g banded solution. By contrast, an area-matched region bounded by bands with constant r_g has a perimeter 1.2P_c. Calculating w/σ for these clusters (Fig. 4(d)) also shows that the typical value for constant κ_g is smaller than for constant r_g, with mean values 4.2σ and 7.3σ respectively.

3.3 Geometric trapping on complex surfaces

Having established that variable curvature on the torus dictates the steady-state shape of the dense phase while creating kinetic barriers to structural transitions, we next investigate whether these competing geometric principles can be leveraged to intentionally trap active clusters. To test this, we consider an hourglass surface consisting of two spheres connected by a narrow neck, as shown in Fig. 5(a). This surface is constructed such that a dense phase covering the typical area fraction observed on the torus could fit entirely on the upper sphere, as indicated by the shaded region in Fig. 5(a). Details of the surface and simulations are given in the Appendix. On the hourglass, the isoperimetric solutions include regions bounded by meridians (lines of constant z), as well as all spherical caps with boundaries entirely within one of the two spheres. The globally optimal location for a disk covering an area is therefore to form a spherical cap in the smaller, upper sphere, as this yields the shortest total boundary length. Alternative solutions, such as a spherical cap on the larger lower sphere, possess perimeters that are at least twice that of the optimal solution. Note that spherical caps on either sphere also have constant r_g.

Fig. 5 (a) Snapshot of the system on the hourglass surface with the coordinate directions marked. The shaded region represents the optimal isoperimetric solution on the surface for the expected dense phase area. (b) Time-dependence of the z-coordinate center of the dense phase cluster for two representative trajectories. In one trajectory, the dense cluster crosses between the lower and upper spheres; in the second, it remains centered on the lower sphere for the entire time sampled. Horizontal dashed lines indicate the boundaries of the neck region. (c) Probability distributions for the dense phase area fraction for clusters centered in different regions of the surface.

To determine whether the dense phase localizes to this thermodynamically predicted location, we track the z-coordinate of the cluster center across trajectories (see Fig. 5(b) for two representative trajectories). Contrary to purely thermodynamic predictions, the dense phase does not preferentially localize to the upper sphere. Instead, we observe that it spends the majority of its time in the lower sphere: in four of the ten recorded trajectories, the dense phase remains in the lower sphere for the entire 200 000τ simulation window, and in the remaining trajectories the cluster transitions between the spheres. In total, roughly 80% of our snapshots have the cluster center below z = 0, despite this lower region accounting for only 65% of the total surface area. The narrow, negative-curvature neck appears to act as a substantial bottleneck for the cluster crossing between the two spheres. When the cluster does cross between spheres, the transition through the neck region occurs over relatively short timescales compared to the time the cluster spends localized within either of the two spheres. These trajectories indicate that the surface effectively creates a two-level system, where the cluster is almost always centered on one sphere or the other for a sustained duration, separated by an effective energy barrier to crossing. We also note that while trajectory data indicates that transitions between spheres are relatively quick, statistical sampling reveals the cluster center occasionally lingers in the neck region, further highlighting how this narrow constriction acts as a kinetic barrier to the cluster transferring between the two spheres.

To understand why the dense phase dynamics deviate from the purely shape-based isoperimetric predictions, we examine the area fraction of the dense phase on the upper and lower spheres. We include only clusters centered above or below the respective sphere equators to avoid the immediate influence of the neck region. As shown in Fig. 5(c), we find that the area distribution for the lower sphere clusters is consistent with observations on the torus, with a mean value ϕ_d = 0.32. On the upper sphere, however, the distribution is much broader, and the mean shifts to ϕ_d = 0.28. We also observe this distribution is sensitive to the choice of center cutoff, with the smaller peak disappearing and the distribution becoming narrower when the cutoff is set further up the upper sphere. This suggests that the dense phase is unable to easily maintain its preferred size on the upper sphere, likely because the nearby high-negative-curvature neck region prevents enough particles from remaining in the cluster. These observations demonstrate how highly localized curvature regimes, in conjunction with global shape-based considerations, create complex energetic landscapes that can trap and localize the dense phase.

4 Conclusions

Our results demonstrate that spatial curvature exerts robust geometric control over the dense phase in motility-induced phase separation, even in parameter regimes where the scale of the Gaussian curvature relative to the particle size (|K|σ²) is sufficiently small that the overall phase boundaries remain largely unaffected.³⁹ While the total area of the dense phase remains independent of the underlying geometry, its morphology and spatial localization are dictated by the surface. By leveraging the distinct geometric properties of curved space, we explored the structural degeneracy present in flat-space MIPS and differing theoretical frameworks governing the dense phase.

Our numerical results and analysis reveal an interesting tension. The steady-state preferred locations of the dense clusters on low-aspect-ratio tori align very well with thermodynamic, surface-tension-driven models.^41,42 The shape of these dense-phase disks is more ambiguous: simulations of larger numbers of particles conform more closely with the thermodynamic picture, whereas those with fewer particles align better with the kinetic framework (see Fig. 4(c) and (d)). By contrast, in the banded regime the cluster shape is consistent with the thermodynamic picture even at the lower number of particles. Our results therefore point towards the cluster shape being ultimately consistent with expectations for thermodynamic, surface-tension-driven shapes in the limit of vanishing curvature (which our larger-N simulations tend toward).

However, while the steady-state shape may be controlled by effective surface tensions, our data suggests that the transition between morphologies is kinetically limited. The shift from a localized disk to a continuous band occurs at an aspect ratio (ξ_c ≈ 2.0–2.1) that is significantly larger than the purely isoperimetric prediction (ξ_c = 1.68). Notably, our simulations with increasing N indicate a shift of the transition point even farther from the isoperimetric prediction. We hypothesize that this discrepancy arises from a substantial geometric energy barrier. At the predicted ξ_c, transitioning from the optimal disk to a band would require massive interfacial fluctuations. Interestingly, calculating the aspect ratio at which a purely kinetic, constant r_g disk of area fraction self-intersects yields ξ_c = 2.10, closely mirroring our observed transition point. This implies that while the cluster strives for an isoperimetric shape, the physical transition to a banded state is dictated by deposition-escape kinetics and the necessity for the fluctuating cluster boundary to self-intersect. We note that while we clearly observe a structural transition on the torus, to our knowledge the equivalent critical aspect ratio for a disk-to-band transition in a flat, periodic system is currently unknown. Characterizing this flat-space counterpart will be a valuable next step to fully isolate the influence of the curved manifold from the effects of the periodic boundary conditions.

While this analysis contrasts the thermodynamic isoperimetric shapes against a baseline kinetic model that assumes a constant geodesic radius, a more complete kinetic theory of MIPS in curved space will likely be more complex. For instance, effective diffusion coefficients are known to depend on the local Gaussian curvature,⁵² which could theoretically introduce spatial variations in the particle influx and escape rates along a cluster boundary. We believe these effects to be subdominant here due to the low curvature scale and the nearly ballistic motion of our Pe_r = 100 active particles. Specifically, escaping particles tend to travel along geodesics for long persistence lengths rather than diffusing back into the cluster, and the geometric variation in available area naturally offsets curvature-induced changes in effective diffusion to maintain a uniform gas influx. While developing an extended kinetic picture—and, indeed, a more complete thermodynamic picture that treats local curvature as a potential source of difference in local energy density^53,54—remains an important theoretical target, we believe the constant-r_g shape serves as a valuable null model for the low-curvature regime studied here. Ultimately, the large-N results strongly favor an interpretation in which the system approaches the isoperimetric solution.

However, while the system ultimately seeks an isoperimetric optimum, our results on the hourglass manifold demonstrate that geometry can easily frustrate this thermodynamic preference. On this surface, both the global isoperimetric solution and simple curvature-dependent diffusivity arguments lead to the expectation that the cluster will be found on the smaller upper sphere; yet, we observe that the dense phase is predominantly trapped on the larger lower sphere. The narrow, negative-curvature neck acts as a kinetic bottleneck that disrupts the local deposition balance and prevents the cluster from reaching its global geometric minimum. This likely stems from the neck region falling outside the low-curvature regime identified in ref. 39, meaning the local phase boundary for MIPS is expected to be shifted in this region. Consequently, this constriction yields distinct cluster trajectories compared to those reported in ref. 39 for a similar surface with a wider neck connecting identical spheres. Whereas their cluster frequently crosses between spheres to undergo curvature-induced size fluctuations, here we observe that the cluster remains kinetically trapped onto one of the two spheres with only infrequent transitions across the neck. These observations suggest that further simulations at larger system sizes may be needed to determine whether these features prevail in lower curvature regimes.

Together, these findings highlight how geometry can serve as a powerful tool for probing the physics of non-equilibrium systems, in this case by breaking degeneracies that exist in flat space. Looking forward, the coupling between curvature-dependent phase behavior, isoperimetric shape preferences, and kinetic bottlenecks opens the door to another flavor of programmable active matter. By designing complex surfaces with specific curvature gradients and topological bottlenecks, it should be possible to engineer stable traps, guide optimal trajectories, and dictate the macro-scale morphology of active clusters. Investigating whether similar geometric control mechanisms are exploited in biological active systems, such as in the curvotaxis of cell collectives or the formation of curved biofilms, remains an exciting future direction.

Conflicts of interest

There are no conflicts to declare.

Data availability

Simulations of active Brownian particles on curved surfaces were performed using the open source software curvedSpaceSim, available at https://github.com/sussmanLab/curvedSpaceSim/blob/SPPBranch/README.md. Trajectory data, sample notebooks including all of the data in the plots, and Morpho code for numerical minimization of the shapes of disks on a torus is available at DOI: 10.5281/zenodo.20169111.

Appendices

A. Analysis on the torus

A.1 Identifying clusters and bands. We identify the dense phase cluster by constructing a nearest neighbor graph for the particles in [Doublestruck R]³ using a Euclidean distance cutoff of 1.1σ. The value of σ means that the Euclidean separation is a very good approximation for the geodesic separation. To determine whether the cluster is a disk or a band, we first map the edges in the largest cluster graph into toroidal coordinate space (θ, ψ) and construct the eight closest periodic images of the cluster edges. If the largest connected component of this combined graph is the same size as the original cluster, the dense phase is a disk; if the size increases, it is a band. Here our focus is primarily on banding around the ψ direction (minor circumference), although bands that span the θ direction at the inner equator sometimes emerge at small ξ. We therefore record separately the direction in which the bands form.

To compare the averaged cluster shape across snapshots, we center each of the clusters in the θ direction. In curved space, the center of an object can be interpreted in different ways. Here we treat the center as the center calculated in coordinate space.

A.2 Dense phase area and perimeter. The area of a cluster on the torus is given by

To calculate this in our system, we construct a concave hull mesh for the centered cluster in torus coordinate space. For disk-shaped clusters, the integration is performed over the region defined by the mesh. If the cluster is a band, the mesh is constructed to include image clusters to eliminate artifacts at the periodic boundaries. The region of integration is then the intersection of the mesh with the primary cell.

The cluster perimeter is calculated by identifying edges on the mesh boundary and calculating the total geodesic lengths of all the edges. For bands, only the edges within the primary cell are included. For small separations, as is the case here for neighboring particles on the edge of the mesh, we use the approximate geodesic separation Δs to speed up calculations. For two points with coordinates (θ₁, ψ₁) and (θ₂, ψ₂), this is given by

where Δθ = θ₂ − θ₁, Δψ = ψ₂ − ψ₁ and ψ_av = (ψ₁ + ψ₂)/2.

We note that in this approach, the exact numerical values of area and perimeter depend on choices of the distance cutoff threshold for connecting points in the concave hull mesh, with higher thresholds leading to a cluster with smoother edges and a slightly larger enclosed area. In torus coordinate space, particles at the inner equator of the torus have a larger separation than those at the outer equator. We therefore choose the cutoff threshold for the mesh to be as small as possible while ensuring all of the particles in the cluster are included within the mesh. At smaller ξ, where this effect is most pronounced, this may lead to a slight over-estimation of the cluster area and underestimation of the perimeter as the resulting mesh is slightly smoothed to incorporate points near the inner equator. However, we believe that the overall trends reported here are robust to reasonable small changes in this threshold.

A.3 Using Morpho for numerical optimization. We use the open-source shape optimization software Morpho⁵¹ to find the isoperimetric solution for a disk centered at an initial location (θ₀, ψ₀) on the torus. In the problem setup, the disk is a circular 2D mesh with a prescribed area A_d = ϕ_dA_T, where ϕ_d = 0.1 is the area fraction and A_T is the final surface area of the torus. The optimization minimizes the perimeter of the disk subject to the constraint A_d is fixed. The disk is constrained to the torus surface through a level-set constraint that imposes an energy penalty for deviation from the surface.

In the initial setup, the disk is placed at a chosen starting location with its center in contact with the torus surface and the disk in the local tangent plane to the torus. To allow the disk to wrap onto the torus properly, we initially use a torus with minor radius r = 2r_d, where r_d is the initial radius of the disk. During optimization, we iteratively allow the disk to find its preferred shape before reducing the torus size at fixed aspect ratio. Here we perform this optimization procedure five times.

Using a relatively small ϕ_d means that the disk remains centered around its initial location, and we are able to more carefully track how the local curvature variation changes the resulting perimeter. Note that for small ξ, the optimal solution near the inner equator is topologically a cylinder, not a disk. We therefore exclude points where the disk overlapped from the plot in Fig. 3.

A.4 Curves of constant geodesic curvature. Consider a curve on the torus γ(s) = (θ(s), ψ(s)), parameterized by arc length s. Using the parameterization of the torus in eqn (3), we can define a set of orthonormal basis vectors ê_θ, ê_ψ that live in the tangent plane to the surface at each point and are parallel to the two coordinate directions. The surface normal ê_n = ê_θ × ê_ψ gives us a complete basis. The tangent vector to γ(s) can then be written as

T̂ = cos β(s)ê_ψ + sin β(s)ê_θ, (9)

where we have introduced β(s) as the angle between the curve and the ê_ψ direction. This obeys

The geodesic curvature of the curve is given by, see e.g. ref. 50

κ_g(s) = T′(s) [ê_n(s) × T̂(s)], (12)

where , N(s) is the unit normal to γ(s) and κ(s) is the curvature of γ(s). For our system, this gives

κ_g(s) = β′(s) − sin ψ(s)θ′(s) (13)

Curves with a constant κ_g can therefore be constructed by numerically solving eqn (14) with κ_g prescribed.

Here, our goal is to determine the value of κ_g that encloses a fixed area. Starting from an initial guess for κ_g, we use Stokes' theorem to calculate the area of the resulting disk. We then use interval bisection to refine this estimate of κ_g until the resulting disk area is within a prescribed tolerance of the target area. The perimeter of the disk is determined as the point where the numerical solution self-intersects.

A.5 Curves of constant geodesic radius r_g. To construct curves of constant geodesic radius, we first find a geodesic γ(t) that starts from the chosen center (θ₀, ψ₀) at an initial angle β₀, measured relative to the ê_ψ direction. The resulting geodesic equations obeyed by γ(t) = (θ(t), ψ(t)) are

We set the initial conditions

ψ(0) = ψ₀, θ(0) = θ₀ (17)

such that γ(t) has unit speed parameterization.

We numerically solve these equations and record the end position after traveling the prescribed geodesic radius r_g along γ(t). To construct the full constant r_g circle, we repeat this procedure for a range of β₀ values spanning the entire circle (typically 60 or 120 values). The perimeter of the circle is calculated as the sum of the geodesic separations between neighboring pairs of sampled end points, and the area is numerically integrated. To find a constant r_g disk with a prescribed area, we start with an initial guess for r_g and calculate its area. We then use interval bisection to refine the estimate until it is within a prescribed tolerance of the target area.

A.6 Cluster interface width. We calculate the mean interface width w of individual clusters using eqn (6). For disk-shaped clusters (ξ = 1.5), we focus on clusters that are centered close to the outer equator, and consider all clusters for which the coordinate center ψ_c obeys ψ_c ≤ ψ_Max, with ψ_Max = π/8. This corresponds to 52% of snapshots for N = 5000 and 94% for N = 20 000.

To find w for a given snapshot, we calculate the area of the dense phase cluster, then construct the corresponding constant κ_g and r_g solutions with their θ-coordinate centers aligned. We then sample points along the perimeter of each of the constant κ_g and r_g curves, and calculate the geodesic distance between each curve and the edge of the cluster. This is calculated as the geodesic distance traveled along a geodesic ray emitted perpendicular to the reference curve at a given point. For constant κ_g, points are sampled evenly along the arc length of the curve. For constant r_g disks, points are sampled at equally spaced angles measured from the center of the disk, while for the constant r_g bands points are sampled evenly around the ψ-direction. In each case, at least 60 points around the reference curves are used to calculate w.

B. Hourglass surface

The hourglass surface consists of two spheres with different radii and offset centers that are connected using a cubic spline that ensures the surface is smooth at the point of connection. In cylindrical coordinates (r, θ, z), the radius of the surface at a given z is given bywhere R_U = 3 and R_L = 2.2 are the radii of the upper and lower spheres, c_U = hR_U and c_L = hR_L with h = 1.4 are the center locations of the two spheres along the z-axis, and the polynomial coefficients a ≈ 0.002, b ≈ 0.45, c ≈ 0.06, and d ≈ 0.46 are numerically determined to ensure the surface is smooth.

This construction means that the upper portion of the surface (z > 0) accounts for 35% of the surface area, meaning that a dense phase cluster at the typical size observed in the torus would fit entirely in this region. Here, we use N = 5000 particles with ϕ = 0.4. The curvature scales on the lower and upper spheres are Kσ² = 0.002 and Kσ² = 0.0038 respectively. The lower sphere is therefore in the “low curvature” regime identified in ref. 39, meaning curvature has minimal effect on phase behavior, while the upper sphere is at the lower end of the “intermediate curvature” regime, meaning the positive curvature may slightly promote clustering. The neck region has negative Gaussian curvature with a minimum value of Kσ² = −0.036, suggesting this region would reduce clustering.³⁹

The z-coordinate center of the clusters is calculated by projecting the [Doublestruck R]³ coordinate center of the cluster onto the closest point on the surface. For a spherical cap, this approach exactly aligns with the coordinate center of the cluster, which is equidistant from all points on the cluster boundary. For our clusters, this approach provides a good approximation for the center, especially when the center is on one of the two spheres and away from the neck region. In the neck region, alternative methods for determining the cluster center could yield slightly different coordinate values. However, we find this approach gives a fast method for calculating the center that provides reasonable values across the entire surface.

Acknowledgements

We thank Tim Atherton and Philipp Schönhöfer for useful discussions. HSA acknowledges funding from the Tarbutton Postdoctoral Fellowship. This material is based upon work supported by the National Science Foundation under Grant No. DMR-2143815.

Notes and references

1. P. Fonda, M. Rinaldin, D. J. Kraft and L. Giomi, Phys. Rev. E, 2018, 98, 032801.
2. D. Marenduzzo and E. Orlandini, Soft Matter, 2013, 9, 1178–1187.
3. S. Busuioc, H. Kusumaatmaja and V. E. Ambruş, J. Fluid Mech., 2020, 901, A9.
4. P. Fonda, M. Rinaldin, D. J. Kraft and L. Giomi, Phys. Rev. E, 2019, 100, 032604.
5. G. Vandin, D. Marenduzzo, A. B. Goryachev and E. Orlandini, Soft Matter, 2016, 12, 3888–3896.
6. E. M. Horsley, M. O. Lavrentovich and R. D. Kamien, J. Chem. Phys., 2018, 148, 234701.
7. L. R. Gómez, N. A. García, V. Vitelli, J. Lorenzana and D. A. Vega, Nat. Commun., 2015, 6, 6856.
8. J. O. Law, A. G. Wong, H. Kusumaatmaja and M. A. Miller, Mol. Phys., 2018, 116, 3008–3019.
9. J. O. Law, J. M. Dean, M. A. Miller and H. Kusumaatmaja, Soft Matter, 2020, 16, 8069–8077.
10. C. Arnold, I. Tahmaz, M.-L. Chapon, H. Maayouf, V. Luchnikov, J.-L. Milan, F. Borghi and L. Pieuchot, BMC Biol., 2025, 23, 296.
11. H. G. Yevick, G. Duclos, I. Bonnet and P. Silberzan, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 5944–5949.
12. W. Xi, S. Sonam, T. Beng Saw, B. Ladoux and C. Teck Lim, Nat. Commun., 2017, 8, 1517.
13. W. Tang, A. Das, A. F. Pegoraro, Y. L. Han, J. Huang, D. A. Roberts, H. Yang, J. J. Fredberg, D. N. Kotton, D. Bi and M. Guo, Nat. Phys., 2022, 18, 1371–1378.
14. L. Pieuchot, J. Marteau, A. Guignandon, T. Dos Santos, I. Brigaud, P.-F. Chauvy, T. Cloatre, A. Ponche, T. Petithory, P. Rougerie, M. Vassaux, J.-L. Milan, N. Tusamda Wakhloo, A. Spangenberg, M. Bigerelle and K. Anselme, Nat. Commun., 2018, 9, 3995.
15. R. Sknepnek and S. Henkes, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2015, 91, 022306.
16. S. Shankar, M. J. Bowick and M. C. Marchetti, Phys. Rev. X, 2017, 7, 031039.
17. C. L. Hueschen, A. R. Dunn and R. Phillips, Phys. Rev. E, 2023, 108, 024610.
18. G. Janzen, E. D. Mackay, R. Sknepnek and D. A. Matoz-Fernandez, Soft Matter, 2026, 22, 474–484.
19. S. Praetorius, A. Voigt, R. Wittkowski and H. Löwen, Phys. Rev. E, 2018, 97, 052615.
20. L. Happel, D. Wenzel and A. Voigt, Europhys. Lett., 2022, 138, 67002.
21. L. Happel and A. Voigt, Phys. Rev. Lett., 2024, 132, 078401.
22. B.-Q. Ai, B.-Y. Zhou and X.-M. Zhang, Soft Matter, 2020, 16, 4710–4717.
23. D. M. Sussman, Phys. Rev. Res., 2020, 2, 023417.
24. M. De Marzio, A. Das, J. J. Fredberg and D. Bi, Phys. Rev. Lett., 2025, 134, 138402.
25. M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys., 2015, 6, 219–244.
26. J. O'Byrne, A. Solon, J. Tailleur and Y. Zhao, Out-of-equilibrium Soft Matter, The Royal Society of Chemistry, 2023.
27. M. E. Cates and C. Nardini, Rep. Prog. Phys., 2025, 88, 056601.
28. Y. Fily and M. C. Marchetti, Phys. Rev. Lett., 2012, 108, 235702.
29. G. S. Redner, M. F. Hagan and A. Baskaran, Phys. Rev. Lett., 2013, 110, 055701.
30. Y. Fily, S. Henkes and M. C. Marchetti, Soft Matter, 2014, 10, 2132–2140.
31. J. Bialké, J. T. Siebert, H. Löwen and T. Speck, Phys. Rev. Lett., 2015, 115, 098301.
32. I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger and T. Speck, Phys. Rev. Lett., 2013, 110, 238301.
33. M. N. van der Linden, L. C. Alexander, D. G. A. L. Aarts and O. Dauchot, Phys. Rev. Lett., 2019, 123, 098001.
34. G. Liu, A. Patch, F. Bahar, D. Yllanes, R. D. Welch, M. C. Marchetti, S. Thutupalli and J. W. Shaevitz, Phys. Rev. Lett., 2019, 122, 248102.
35. W. J. Ridgway, M. P. Dalwadi, P. Pearce and S. J. Chapman, Phys. Rev. Lett., 2023, 131, 228302.
36. C. Anderson and A. Fernandez-Nieves, Nat. Commun., 2022, 13, 6710.
37. T. Knippenberg, A. Jayaram, T. Speck and C. Bechinger, Phys. Rev. Lett., 2024, 133, 048301.
38. P. Iyer, R. G. Winkler, D. A. Fedosov and G. Gompper, Phys. Rev. Res., 2023, 5, 033054.
39. P. W. Schönhöfer and S. C. Glotzer, Soft Matter, 2022, 18, 8561–8571.
40. A. P. Solon, J. Stenhammar, R. Wittkowski, M. Kardar, Y. Kafri, M. E. Cates and J. Tailleur, Phys. Rev. Lett., 2015, 114, 198301.
41. A. Patch, D. M. Sussman, D. Yllanes and M. C. Marchetti, Soft Matter, 2018, 14, 7435–7445.
42. L. Langford and A. K. Omar, Phys. Rev. E, 2024, 110, 054604.
43. L. Li, Z. Sun and M. Yang, Proc. Natl. Acad. Sci. U. S. A., 2025, 122, e2505010122.
44. J. Tailleur and M. E. Cates, Phys. Rev. Lett., 2008, 100, 218103.
45. R. Soto and R. Golestanian, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2014, 89, 012706.
46. R. Soto, M. Pinto and R. Brito, Phys. Rev. Lett., 2024, 132, 208301.
47. T. H. Webb and D. M. Sussman, Comput. Phys. Commun., 2025, 311, 109545.
48. M. Ritoré, Commun. Anal. Geom., 2001, 9, 1093–1138.
49. A. Cañete, Indiana Univ. Math. J., 2007, 56, 1629–1660.
50. R. D. Kamien, Rev. Mod. Phys., 2002, 74, 953.
51. C. Joshi, D. Hellstein, C. Wennerholm, E. Downey, E. Hamilton, S. Hocking, A. S. Andrei, J. H. Adler and T. J. Atherton, Nat. Comput. Sci., 2025, 5, 170–183.
52. J. Faraudo, J. Chem. Phys., 2002, 116, 5831–5841.
53. L. Braverman and D. R. Nelson, Phys. Rev. E, 2025, 112, 044221.
54. G. Tarjus, F. Sausset and P. Viot, Adv. Chem. Phys., 2011, 148, 251–310.

---

Footnote

† These authors contributed equally to this work.

---