Hyperuniformity and phase enrichment in vortex and rotor assemblies

Ensembles of particles rotating in a two-dimensional fluid can exhibit chaotic dynamics yet develop signatures of hidden order. Such"rotors"are found in the natural world spanning vastly disparate length scales - from the rotor proteins in cellular membranes to models of atmospheric dynamics. Here we show that an initially random distribution of either ideal vortices in an inviscid fluid, or driven rotors in a viscous membrane, spontaneously self assembles. Despite arising from drastically different physics, these systems share a Hamiltonian structure that sets geometrical conservation laws resulting in distinct structural states. We find that the rotationally invariant interactions isotropically suppress long wavelength fluctuations - a hallmark of a disordered hyperuniform material. With increasing area fraction, the system orders into a hexagonal lattice. In mixtures of two co-rotating populations, the stronger population will gain order from the other and both will become phase enriched. Finally, we show that classical 2D point vortex systems arise as exact limits of the experimentally accessible microscopic membrane rotors, yielding a new system through which to study topological defects.

Two-dimensional (or nearly so) fluid flows show rich and complex vortical dynamics. These can arise from flow interactions with boundaries (1,2), the inverse cascades of 2D turbulence (3)(4)(5), from Coriolis force dominated atmospheric flows (6), and from quantization effects in super fluid He-II (7,8). Point vortices have long been staples for the modeling of such inertially dominated inviscid flows. Kirchoff (9) was the first to describe point vortices using a Hamiltonian framework and his work was extended by many others [e.g. (10)(11)(12)], notably, Onsager (13) in his statistical mechanics treatment of 2D turbulence as clouds of point vortices.
Remarkably, structurally identical Hamiltonian and moment constraints can arise in the microscopic viscously-dominated realm from a strict balance of dissipation with drive on immersed rotating objects. These objects include models of interacting transmembrane ATP-synthase "rotor-proteins" (14)(15)(16), and the planar interactions of rotors -microscopic particles driven to rotate by an external torque (17,18). We refer to such systems as BDD systems, as in balanced drive and dissipation. In modeling rotational BDD systems other physical effects may also come into play, such as steric interactions, that can yield interesting complexities (16). Interacting assemblies of driven-to-rotate particles has become an area of intensifying interest in the active matter community (17)(18)(19)(20)(21)(22)(23).
Here we study both point vortices and a BDD rotor system of rotationally-driven microscopic particlesmembrane rotors -immersed in a flat membrane. We show that in both systems, their Hamiltonian conservation laws lead to distinct structural states -hyperunifor-mity, phase enrichment and crystallization (see Fig. 1), not yet observed for either system. We use the Hamiltonian to derive a bound for spatial correlations requiring hyperuniformity. We demonstrate numerically that rotational dynamics robustly self-assembles particles into a disordered hyperuniform 2D material; This self-assembly is insensitive to the details of the hydrodynamic interactions, steric repulsion, or the presence of impurities in the form of different rotation rates. At steady state, the long wavelength configuration is characterized by an isotropically vanishing structure factor, S(q → 0) → 0 (where q is the wavevector), leading to an isotropic band-gap (24)(25)(26).
In classical mechanics, symmetries of the Hamiltonian H restrict the phase-space of the conjugate variables, position and momentum. However, in 2D point vortex or BDD point rotor systems, the conjugate variables are the actual spatial coordinates of the ensemble {x i } and {y i }. The conservation laws are therefore geometrical in nature, bounding the proximity and distribution of the particles. For both point vortices and membrane rotors, as well as for a myriad of other 2D rotating systems (17-20, 22, 27), the dynamics are dictated by Hamilton's equations, where ∂ ⊥ i = (∂y i , −∂x i ), v i is the velocity of rotor i, and Γ i is the circulation (proportional to the magnitude of the torque for rotors). Our finding, as we will show, is that the spatial arrangements of point vortices, as measured by S(q), are dictated by the Hamiltonian, To derive Eq. 2 and to find the Hamiltonian of N particles, we first describe the flow due to a single vortex in arXiv:2103.00296v1 [cond-mat.soft] 27 Feb 2021 Three different structural states of 2D vortices/rotors -hyperuniformity for Euler point vortices (A) and QG rotors/surface rotors (B), (C) phase enrichment induced by circulation differences where green (black) represents vortices of high (low) circulation, and (D) crystallization arising from hydrosteric interactions. The insets of (A), (B) and (C) show the structure factor, S(q). In (A) and (B) S(q) decays to zero at small q, indicating that the distribution is hyperuniform. In (C) the structure factor shows the six distinct peaks of a hexagonal lattice.
an ideal Euler fluid and show its equivalence to a point rotor in a viscous membrane. We then use the linearity of the equations to extend the result to the many-body case. An ideal point vortex is given by a singular vorticity, ω = ∇ × v = δ(r). A 2D incompressible fluid can be described using a stream function Ψ such that the velocity, v, is given by v = ∂ ⊥ Ψ. This equation, combined with the equation above gives, Ψ = − 1 2π log r (12). The flow, v(r), therefore, scales as 1/r, where r = |r|.
We switch now to a point rotor in a viscous membrane, driven by an external torque τ . Following Saffman and Delbrück's seminal work (28), and many others that followed (29)(30)(31), we assume that the membrane is incompressible (∇·v = 0), and that inertia is negligible. Under these assumptions, the Stokes momentum conservation equation for the membrane reads, where v is the 2D velocity in the plane of the membrane, u ± is the 3D flow in the outer fluids, η 2D is the 2D viscosity, and η 3D is the viscosity of the outer fluids. The second term on the right hand side is the surface shear stress of the outer fluids, and the third term is the force due to a rotating point object. There is no pressure contribution when the motion is purely rotational. This equation is coupled to the equations of the outer fluids. It is easy to solve the above equations using a 2D Fourier Transform where Γ = τ /η 2D , and λ = η 2D /2η 3D is the Saffman Delbrück length. At small distances (r λ) momentum travels in the plane of the membrane. At large distances (r λ) momentum travels through the outer fluid as well (32,33). In real space Ψ(r) = 1/4(H 0 (r/λ) − Y 0 (r/λ)), where H 0 and Y 0 are zeroth order Struve function and Bessel function of the second kind respectively.
In the limit of small distances, r λ, the stream function is, Ψ ≈ − 1 2π log r, i.e. exactly the same as for an ideal point vortex. In the opposite limit, r λ, the stream function becomes Ψ = 1 2πr as in quasigeostrophic (QG) flows -atmospheric or oceanic flows coming from gradients in pressure coupled to the Coriolis force (34), or driven rotors on the surface of a fluid (21). A membrane rotor, therefore, transitions from a point vortex for Euler at small distances to that of QG flow at large distances. The velocity is given by derivatives of Ψ and is thus proportional to 1/r (1/r 2 ) in the limit of small (large) distances (see Fig. 2B). For simplicity, we work primarily in the limit of small distances, r λ, since in this limit the dynamics in a membrane converge with those of point vortices (many results still apply to the more general case). In what follows, we will use "point vortices" when there are only hydrodynamic interactions and "rotors" when the particles have steric interactions in addition to hydrodynamic ones.
The dynamics of N point vortices follows from the The Hamiltonian depends on the conjugate variables r i = (x i , y i ), [normalized by the circulation |Γ i | sgn(Γ i )], i.e. the positions of the vortices (12). The symmetries of the Hamiltonian correspond to conservation laws (36). In this case, we have symmetries with respect to translation in time, space, and rotation, corresponding to conservation of the Hamiltonian itself, and of the first and second moments of the distribution, L = i Γ i r i (= 0 wlog), and M = i,j Γ i r 2 i . Thus, the initial area cannot change dramatically, particles cannot drift to infinity since the second moment is fixed, nor can they collapse to a point since the Hamiltonian is conserved. These properties are readily observed in simulations. Figure 2D shows typical trajectories of 200 membrane rotors. The initial distribution is random in a predefined finite area, and the dynamics are chaotic (37). The final configuration occupies nearly the same region of space as the initial configuration does, and the conservation laws hold The velocity field due to a membrane rotor (solid line) which scales as a point vortex v ∼ 1/r at small distances (dotted), r/λ 1, transitioning to a QG behavior at large distances v ∼ 1/r 2 (dashed). (C) Contour dynamics of an ellipse with radii ratios r l /rs ≤ 3, where r l (rs) is the major (minor) axis. Starting from the same contour, the dynamics differ according to the radius relative to the SD length. Blue is in the limit r l λ. In this limit the ellipse is rotating as a rigid body, as predicted by Kelvin (35) for an elliptic patch in an Euler fluid. Black is in the limit r l λ, no longer conserving its shape since the large distance flow is in the quasigeostrophic regime. (D) 200 point membrane rotors, blue is the initial random configuration, black is the final configuration. Solid line shows typical trajectory of an individual vortex. Note that the area did not change considerably since the system of vortices is self-bounding. (E) the relative error in H and M over a few cycle times, tc.
to high precision in our simulations, as shown in Fig. 2E. This self confining property of vortex dynamics has further consequences, as we now show.
Hyperuniformity. Hyperuniformity is the suppression of density-density fluctuations at small wavenumbers (or correspondingly, at large distances) (38)(39)(40). Disordered hyperuniformity can emerge due to short ranged interactions such as those that arise in sheared suspensions (27,41,42), jammed materials (43), and for spinning particles (44). Here we will show hyperuniformity emerging from long ranged interactions, similar to its emergence in sedimentation of irregular objects (45). A good way to characterize hyperuniformity is the structure factor, defined as S(q) = N −1 | ρ(q)| 2 , where ρ(r) = i δ(r − r i ) is the coarse grained density. In a hyperuniform material, S(q) goes to zero as a power law at small wavevnumbers. We argue that point vortices must be hyperuniform due to the conservation of the Hamiltonian. For a density of rotors, the Hamiltonian is given by H[ρ(r)] ∼ Γ 2 2 dr dr ρ(r)ρ(r )ψ(|r − r |). Using the convolution theorem, we find a general relation between the Hamiltonian and the structure factor In the case of point vortices, Ψ(q) = 1/q 2 , which gives Eq. 2. For the integral of Eq. 2 to converge in 2D, S(q) ∼ q α near the origin, and we must have α > 0. In other words, an ensemble of point vortices is hyperuniform (a similar calculation in the QG limit, where Ψ = λ/q, yields α > −1). Figures 3B and 4C, show an apparent α ∼ 1.3 scaling for point vortices, consistent with the above argument.
Using simulations we show that a set of N vortices, uniformly distributed within a radius R, evolves to a disordered steady-state with a hidden order visible to the naked eye (compare Figures 3A left and right). We quantitatively characterize the system in steady-state in three ways: (1) The structure factor. At steady-state S(q) shows a distinct cavity, at q ≈ 0, S(q) → 0, for both points vortices (Fig. 3A) and rotors (Fig. 3C). All simulations produce a hyperuniform arrangement. (2) Perturbations. We demonstrate that hyperuniformity is robust under different perturbations, be it in the form of numerical errors, repulsive interactions, or impurities (in the next section). For point vortices, the steady state appears later and later as the timestep is decreased, suggesting that perturbations are necessary for convergence, here very small but persistent timestepping errors (46). Adding steric interactions, hyperuniformity appears on a timescale that is independent of the timestep. Moreover, with steric interactions, as the area fraction φ of the particles is increased, the system transitions from disordered hyperuniform, to an ordered hyperuniform hexagonal lattice at φ ∼ 0.5, as can be seen in Fig. 3C. The inset of Fig. 3B shows the averaged structure factor where at intermediate area fractions we see Percus-Yevick type features for the structure factor of disks (47). (3) The returnity. We observe that at late times the ensemble of point vortices rotates almost as a rigid body and each particle goes back to its position at the previous cycle. We measure particle deviations by what we term the "returnity" (see Fig. 3D for details). The system may seem to have reached an absorbing state, but the motion of vortices over many cycles is still chaotic.
Rotation induced phase enrichment. We now show that for mixed populations of fast and slow rotating particles, there is phase enrichment of both populations and hyperuniformity of the fast ones. Consider a mixture of two equally numbered populations (ρ l = ρ h at t = 0) initially placed within the same radius R. ρ l rotates slowly with Γ l Γ h , where Γ h is the circulation of the second population. Figure 4A  a disk of only slightly smaller size than their initial area (Fig. 4B). The slow particle distribution shows a significant expansion. In addition, there is a striking difference when comparing the independently computed structure factors of these two populations, the fast vortices are hyperuniform with the same scaling as before, S(q) ∼ q 1.3 , whereas the slow ones show no signs of hyperuniformity (Fig. 4C). This difference is dramatic enough to be visible in a cursory examination of the separate distributions; see Fig. 4A.
Using a heuristic model, we show that the conservation laws allow two solutions at steady-state. In one solution, the two populations remain confined to a circle of the same radius. In the second solution, the radius of the slower population expands, while the radius of the faster population contracts. We then show that the segregated solution is the one that maximizes the number of states in the system. For simplicity, we assume that the final steady states are uniform (not true for the slow vortices as is clear from Fig. 4B). There are two possible solutions where H and M are conserved -in the first, the initial radius, R, does not change; in the second, the radius of the fast vortices slightly decreases to R h , allowing the slow vortices to expand to a larger radius R l given by Fig. 4D). Linearly expanding in 1/γ, we find that R h R(1−β/γ) for the high circulation vortices, where β is a positive prefactor of order 1. The slow vortices asymptote to R l R √ 1 + 2β +O(1/γ). The simulation results indicate that the outer radius indeed asymptotes to a larger valued constant as γ increases and does not increase indefinitely (see Fig. 4D).
A solution with two different radii is therefore possible and is indeed observed at large circulation ratios. Such a solution is favored entropically since it maximizes the available states. Asymptotically at large γ, the main entropical contribution is volumetric, ∆S volume = 2N log(R final /R initial ). Since the high circulation vortices hardly change radius, R h γ→∞ − −−− → R, the change in entropy is coming mainly from the expansion of the low circulation vortices and is given by ∆S total ∼ N log(1+2β) > 0. Coupling the two populations allows one population to expand where before it was bounded (48). The situation is analogous to depletion interactions, where the net entropy of a system increases by condensing the large particles allowing for the small particles to explore a larger volume (49).
A simple way to estimate the entropy in a system is by using LOSSLESS compression, as suggested by Refs. (50,51). Compressing plots of particle positions in a system of 10,000 point vortices with circulation ratio Γ h /Γ l = 128 shows an increase in file size for ρ l and a decrease for ρ h , while the combined system is increasing, see Fig. 4E.
Discussion. We have shown that driven particles in a membrane or a soap film, as well as point vortices in an ideal 2D fluid, have geometrical conservation laws which limit their distribution. These conservation laws dictate different possible structural states -namely hyperuniformity and phase enrichment. We have shown that hy-peruniformity is robust to several forms of perturbations whether arising due to numerical errors, steric interactions, or impurities in the form of low circulation vortices. For rotors with steric interactions, the unbounded ensemble crystallizes into a hexagonal lattice when the area fraction φ 0.5 (see also (16)). We have limited the discussion to membrane rotors and vortices, but the results hold for other settings in which mass is conserved in the 2D plane, e.g. particles at the surface of a fluid.
What is especially interesting about our particular BDD system is its potential for experimental realizability, its moment and Hamiltonian structure, and that its near-field interactions (i.e. below the Saffman-Delbruck length) are identical to those of Euler point vortices. Further, the far-field interactions of membrane rotors are identical to those of point vortices of the semiquasigeostrophic equations (34,52,53) used to model atmospheric flows. Thus, to observe the interesting dynamical features we describe, one does not need to go to the atmospheric scale, or cool a fluid to near-zero temperature. In principle, one can simply observe microscopic particles on a soap film, in smectic films, a membrane, or even at the surface of a fluid (18,21,54,55). Methods.
Simulations. Simulations were performed in Python. Random initial configurations within the unit disk were found by rejection sampling (points in the unit rectangle were sampled uniformly, transformed to the rectangle [−1, 1] 2 , and those with r > 1 were discarded). The initial Hamiltonian H 0 is computed at t = 0, and the relative error (t) = |Ht − H 0 |/H 0 is monitored as a measure of fidelity. For simulations of rotors (i.e. with steric repulsion), a 5th order explicit Runge-Kutta method based on the Dormand-Prince scheme (56) with a fixed timestep size of δt = 10 −7 was used. Long integration times were required for simulations of point vortices, and for these simulations an exlpicit eighth-order adaptive method based on the Dormand-Prince scheme (57, 58) was used, with both relative and absolute tolerances set to 10 −6 . The specific implementation of the scheme used was the DOP853 method of scipy.integrate (59). For simulations of 10,000 point vortices with Γ = 2π, (t) < 1.6 × 10 −3 up to t ≈ 16, 000 cycles, while for simulations with 5,000 vortices with Γ = 2π and 5,000 vortices with Γ = 256π, (t) < 5 · 10 −3 up to t ≈ 10 cycles. Time is normalized by the average cycle time, tc ≈ 4π 2 R 2 / i Γ i , where R is the initial radius.
Steric interactions were taken as the repulsive part of a harmonic potential, i.e. for two particles whose centers are distance r i apart, F = −kr ij if r ij < 2a and zero otherwise. The use of a harmonic potential, rather than a sharp step function for hard core particles, provided improved numerical stability and convergence. A large k value was chosen to ensure no overlap between particles, k = 1·10 6 , for particles of size a = 0.01.
Structure factor. To accurately compute the structure factor S(q) we use a type-1 non-uniform fast-Fourier transform (60). Explicitly, points are restricted to a windowing region which is confined entirely within the unit disk. The frequencies ρ(q) are computed for the first 512 modes in each direction, and the average value (i.e. ρ(0)) is set to 0. This results in structure factors in the plane, such as those shown in Fig. 3. Except in those cases where crystallization occurs, these structure factors are azimuthally isotropic. To summarize this information, the angular average over the structure factor was calculated by slicing the result to 1000 equal bins between q min and qmax and taking the mean of the results that fell within each slice.
Compression. A plot of the positions of the point vortices was compressed using PNG with AGG backend. Each vortex was plotted by a single pixel. The total size of the plots was kept fixed in time. The figure size was chosen to minimize overlap between neighboring vortices but maintaining a computationally accessible file size.