Spontaneous self-propulsion and nonequilibrium shape fluctuations of a droplet enclosing active particles

Active particles, such as swimming bacteria or self-propelled colloids, spontaneously assemble into large-scale dynamic structures. Geometric boundaries often enforce different spatio-temporal patterns compared to unconfined environment and thus provide a platform to control the behavior of active matter. Here, we report collective dynamics of active particles enclosed by soft, deformable boundary, that is responsive to the particles’ activity. We reveal that a quasi two-dimensional fluid droplet enclosing motile colloids powered by the Quincke effect (Quincke rollers) exhibits strong shape fluctuations with a power spectrum consistent with active fluctuations driven by particle-interface collisions. A broken detailed balance confirms the nonequilibrium nature of the shape dynamics. We further find that rollers self-organize into a single drop-spanning vortex, which can undergo a spontaneous symmetry breaking and vortex splitting. The droplet acquires motility while the vortex doublet exists. Our findings provide insights into the complex collective behavior of active colloidal suspensions in soft confinement. The presence of geometric boundaries is known to affect the collective behaviour of active particles. Here, the authors unravel exotic patterns of self-propulsion and non-equilibrium shape fluctuations in a system of active particles enclosed in a droplet and interacting with its soft boundaries.

A ctive (self-driven) particles such as motile colloids present novel opportunities for the engineering of smart materials that can self-heal or change properties on demand 1,2 . The reconfigurability of active materials stems from the propensity of active particles to self-organize into large-scale dynamic structures that can be modulated by external cues such as electric or magnetic fields 3-10 , light 11,12 , or chemical reactions [13][14][15] . The spatio-temporal patterns can be also manipulated by geometric boundaries. For example, while unconfined suspensions of bacteria exhibit turbulent-like flow [16][17][18][19] , confinement into a long and narrow macroscopic "racetrack" geometry stabilizes bacterial motion into a steady unidirectional circulation 20,21 , and when subjected to 2D arrays of vertical pillars arranged in a square pattern bacterial suspensions transforms into a lattice of hydrodynamically bound vortices with a long-range antiferromagnetic order 22 . Inside of a droplet, the bacteria form a macro-scale bidirectional vortex 23,24 . Similar behaviors are also observed with synthetic microswimmers. For example, the Quincke rollers, which are powered by a spontaneous electro-rotation of a dielectric sphere exposed to a static electric field, self-organize into a long, polar band and undergo directional motion in the racetrack microfluidic channel 25 , while in strong confinement (rectangular geometries), the band state is replaced by a single macroscopic vortex 26 .
A confinement that is responsive to the particles activity adds another degree of freedom that can unlock novel collective states and dynamics. Active matter confined in drops can exhibit spontaneous symmetry breaking 27,28 leading to droplet motility 29,30 . Active droplets also mimic cell behaviors such as growth, division, and reshaping 31,32 . Bacteria and self-propelled colloids encapsulated in a droplet or a vesicle (bilayer membrane sac) drive strong shape deformations [32][33][34][35] and can cause net motion 33,[36][37][38][39][40] .
Here, we explore the relation between the particle's activity, deformations and motility of the soft confining container. As active particles we employ the Quincke rollers since their speed and locomotion pattern can be easily manipulated 41,42 . We use a soft container comprised of a liquid droplet sandwiched between two electrodes, which creates quasi two-dimensional geometry. We find that at low activity (quantified by the speed of the rollers) the droplet contour fluctuates, while the droplet stays nearly stationary, and the rollers self-organize into a vortex spanning the whole system. Increasing the activity leads to a growth of the shape fluctuations that exhibit a power spectrum consistent with active fluctuations driven by particle-interface collisions. Coupling of activity and soft boundary fluctuations often results in bursts of the droplet translation in a randomly selected direction. The net propulsion is driven by a spontaneous formation of a vortex doublet composed of two counter-rotating vortices.

Results
Particle dynamics inside the droplet. The experimental system consists of 40 μm polystyrene spheres dispersed in hexadecane/ AOT medium (see "Methods" section for details). A small volume of the solution (approx. 5 μL) is sandwiched between two indium tin oxide (ITO) coated glass electrodes spaced 150 μm apart to form a liquid bridge with a high aspect ratio that produces a quasi-two-dimensional droplet (Figs. 1a and 2a). The particles are allowed to sediment on the bottom electrode before the application of a uniform DC electric field E between the electrodes. Above a threshold magnitude E Q particles develop a steady rotation due to the Quincke effect 25,43 (see Supplementary note I for an overview of the phenomenon) and start to roll on the bottom surface.
Hydrodynamic and electrostatic interactions between the rollers promote alignment of their translational directions and result in a formation of multiple transient flocks and vortices of particles (Fig. 1b, c). Eventually, multiple vortices and flocks merge to form a macroscopic global vortex inside the droplet, see Fig. 2b, c. Similarly to vortices formed by magnetic rollers 44,45 , the vortex spontaneously selects its handedness (clockwise or counterclockwise) that changes from experiment to experiment. The particles velocity fields in the droplet reveal a dramatic change in its appearance compared to rollers confined by a solid boundary 26 , where the rollers accumulate near the confining interface. In the droplet system, the rollers are distributed throughout the droplet interior being more packed in the center of the droplet and less densely packed towards the drop edge. This is most clearly seen in Fig. 2c as one bright blur in the center of the droplet-a big crowd of rollers rotating as one vortex-and several bright blurs around-flocks orbiting in the direction of the rotation. Equal direction of rotation for the whole system is confirmed by a single central peak in the vorticity field of The equilibrium droplet contour in the absence of activity (below the onset of Quincke rotation) is a circle (Fig. 2a). Once the rollers become motile, the interface begins to fluctuate and during the process of vortex formation the droplet shape can become very non-circular (see Fig. 2b and Supplementary Movie 1). Even when the vortex is formed the shape continues to fluctuate (see Fig. 2c, Supplementary Movie 2, and Supplementary Fig. S1). We quantify the droplet deformations by the asphericity parameter Δ (see "Methods" section), with Δ = 0 corresponding to a perfect circle. Figure 2d shows that once the system is energized, the droplet undergoes pronounced shape deformations until a macroscopic vortex is formed and Δ decreases back to nearly pre-activation values (nevertheless shape fluctuations are still present). The gradual growth of the macroscopic vortex is illustrated by a correlation length, r corr , defined as the first zero crossing of the spatial velocity correlation function, C norm (r) (Fig. 2d inset). r corr exhibits a monotonic growth while the individual rollers organize into flocks and transient vortices (the velocity fluctuations phase), and it reaches a plateau when the globally correlated state (vortex) is formed.
Nonequilibrium shape fluctuations. The shape of the active droplet undergoes strong fluctuations with amplitude reaching 10-15% of the initial droplet radius, see Fig. 3a. The deviation from the equilibrium circular shape, h(s, t) = r s (s, t) − R drop , where r s is the droplet interface position at arclength s, is fitted with Fourier series, hðs; tÞ ¼ ∑h q ðtÞ expðiqsÞ (see "Methods" section). The power spectrum of the shape fluctuations exhibits a power-law dependence on wavenumber q, 〈|h q | 2 〉~q −4 (Fig. 3b). Thermal fluctuations of the droplet interface would exhibit a power spectrum with q −2 capillary scaling for surface-tensioncontrolled shape relaxation, suggestive that the origin of the fluctuations is non-equilibrium.
The temporal behavior of the shape fluctuations also indicates out-of-equilibrium dynamics. In the active state, the relaxation time obtained from the time autocorrelation function, , displays a power-law decay with q with exponent close to −3/2, see Fig. 3c. In equilibrium systems, thermal shape fluctuations with mean-squared amplitudẽ q −4 are exhibited by bilayer membranes and semiflexible polymers and are due to bending rigidity. The bending forces drive relaxation with rates~q 3 for membranes and~q 4 for polymers [46][47][48] . However, in non-equilibrium systems the excitation and relaxation need not to have same driving forces. The q −4 -law of the power spectrum is observed for vesicle shape fluctuations due to particle-interface collisions 32,35 and active ring polymer 49 . In the latter system, the relaxation times are predicted to be~q −2 and q −4 for tension and bending controlled modes, respectively. In our system, the fluctuations are driven by the strong flows generated by the Quincke rollers. The deformations are opposed by the surface tension and the relaxation rates are set by dissipation due to the motion of the contact line 50,51 . An estimate based on the balance of surface tension γ and viscous dissipation by viscosity μ gives for the relaxation time τ q = μ/γθ 3 q, where θ is the contact angle 50 . This 1/q dependence of relaxation time is weaker than the experimentally observed one, suggesting that the relaxation is not purely driven by contact line elasticity and likely impacted by the activity.
To quantify the non-equilibrium nature of the fluctuations, we test for broken detailed balance in the transitions between microscopic configurations 52 (see "Methods" section). The configurations correspond to the shapes defined by different Fourier modes. In equilibrium, it is equally likely for the forward and backward transition to occur between any two different Fig. 1 Quincke rollers in a droplet. a A sketch of the experimental system: A small amount of weakly conductive liquid (hexadecane) surrounded by air forms a bridge between two planar electrodes. Inside of the quasi-two-dimensional drop are polystyrene spheres that start to roll upon application of a uniform DC electric field, E. Quincke rollers initially move in random directions along the surface of the bottom electrode. Blue and green arrows indicate the direction of translation and rotation of the rollers, respectively. The image shows the top view of the drop and enclosed rollers. b The monolayer of rolling particles generates an essentially two dimensional flow. Streamlines reveal the formation of transient vortices throughout the droplet. c Snapshot of the velocity field inside the droplet as determined by particle image velocimetry (PIV). Experimental parameters: packing fraction of the Quincke rollers φ = 0.58 ± 0.11, droplet area is A drop = 11.05 ± 0.03 mm 2 , the driving electric field E = 3.193 ± 0.007 MV m −1 . . Rollers velocity and vorticity fields in the droplet. Both clockwise (blue) and counter-clockwise (red) vorticity is simultaneously present indicative of transient vortices. Even though the time-variation of the asphericity Δ appears periodic, this behavior is specific for this particular experiment and it is not universally observed. c Globally correlated state (vortex phase). Velocity and vorticity fields indicate the presence of a single macroscopic vortex. d Evolution of the droplet shape characterized by the asphericity, Δ (cyan circles, left scale). The right scale presents the time evolution of the correlation length r corr on the same dataset, as defined in the inset from C norm (r), the normalized spatial correlation of velocity field, first zero crossing, normalized by the drop radius R drop . The inset shows C norm (r) computed from the data in the velocity fluctuations (#) and vortex (*) phases. The dashed lines delineate the regimes of velocity fluctuations, transition to vortex, and developed vortex. The error bars are the standard error of the mean and the purple line is a running average of the data points to guide the eye. When error bars are not visible they are smaller than the symbols. Experimental parameters: rollers packing fraction φ = 0.18 ± 0.04, droplet area A drop = 13.56 ± 0.03 mm 2 and electric field strength E = 5.104 ± 0.008 MV m −1 .
Fourier modes. A non-equilibrium system, however, would display a probability flux in the phase space of shapes. Figure 3d (Fig. 3f). The broken detailed balance analysis also reveals that the modes gradually go out of equilibrium, starting with the short wavelengths. The longest wavelength modes are in equilibrium early in the velocity fluctuation phase (Fig. 3d, f) and become nonequilibrium in the vortex phase (Fig. 3e, f). This reflects the structure evolution: rollers cluster into aggregates with growing size that eventually become the droplet-spanning vortex.
Activity enhancement of the shape fluctuations. Thermallydriven droplet shape fluctuations are negligible in our system due to strong surface tension. For an oil/air interface the surface tension is γ = 10 mN m −1 , which corresponds to interfacial energy far exceeding the thermal energy,~10 12 k B T, for a droplet`with radius of 1 mm; the amplitude of the fluctuations calculated from 〈|h q | 2 〉 = k B T/γq 2 L 2 , where the contour length is L = 2πR drop , is below a nanometer even for the lowest wavemode. However, in the active state the rollers generate flows that can be strong enough to overcome surface tension and deform the interface. An individual roller with radius a moves with speed a _ G, where the generated flow strain rate is _ (see SI for details). In our system, the Maxwell-Wagner time is τ mw~1 ms.
The flow can exert stress on the interface in the order of μ _ G $ 1 Pa, which is comparable to the capillary stress γq. Since the flow is created by the rollers, the active shape fluctuations are expected to increase if either the rollers velocity or number increases.
To quantify the effects of the activity, we define an effective energy of the system from velocity fluctuations 53 where the index N runs over the individual rollers, and u N = v N − v is the difference between the individual roller velocity, v N , and the macroscopic flow, v. Experimentally, however, we have access neither to the individual particle trajectories nor the detailed hydrodynamic flow, therefore we consider v to be the instantaneous velocity of the droplet center and v N the velocity field pixel from PIV velocity fields of the particles. PIV procedure is based on image contrast correlations and detects only the movement of the particles in the droplet making the velocity field a good approximation for the actual particle velocities as long as the packing fraction φ is large. Assuming that all particles have equal mass m N = m p and replacing the summation by the temporal average value of the square of the velocity fluctuations hu 2 N i yields for the energy density where ρ is the buoyant density of particles. In our system the velocity of individual particles depends on the driving electric field  The hu 2 N i linearly increases with ðE=E Q Þ 2 À 1 (see Fig. 4a inset) and, as a result, the energy density e is proportional to ðE=E Q Þ 2 À 1 (Fig. 4a). The energy density also follows a linear dependence on the packing fraction φ as predicted by Eq. (2) (Fig. 4b), demonstrating that addition of active agents is directly reflected in increased energy injection into the droplet.
Activity is known to enhance fluctuations of ring polymers 49 and vesicles [55][56][57][58][59] . However in these systems the active particles are located at the interface, while in our system the active particles are in the bulk, and the interface deforms in response to both direct particle collisions and flow due to the collective motion of the rollers. The magnitude of the active shape fluctuations in our system scales as 〈|h q | 2 〉 = αq −4 (Fig. 3b), in agreement with the scaling predicted for fluctuations due to particle collisions with a fluid interface governed by surface tension 35 . The active pressure associated with the particle collisions is a volumetric energy density, thus the theory suggests that α should scale quadratically with the energy density, which is in agreement with our experimental results (Fig. 4c). A quadratic increase of the fluctuations magnitude with activity, quantified by the Peclet number, is also predicted for the ring polymer system 49 . It is intriguing that different systems (e.g., our system is dense and particle motions are strongly correlated, while the model of ref. 35 considers uncorrelated collisions) display qualitatively similar behavior.
Spontaneous droplet self-propulsion. The shape fluctuations are accompanied by a motion of the droplet as a whole. The behavior of the mean square displacement (MSD) 〈Δr 2 〉 for the center of the droplet exhibits a typical diffusive behavior at the long time scales (see Fig. 5a), and it is in agreement with the behavior observed with droplets enclosing bacteria 38 or active nematics 30 , and predicted by simulations for microswimmers in a vesicle 33 . As activity increases, the droplet shape fluctuations grow and the system occasionally undergoes a symmetry-breaking instability. The global vortex spontaneously splits into two counter-rotating entities that drive a significant elongation of the droplet (as illustrated by the changes in the droplet perimeter shown in Fig. 5b  to the direction of the droplet travel (corresponding angle to β ≈ π/2), see the inset of Fig. 5b) in contrast to the regular droplet excursions that do not show correlations between β and velocity of the droplet (see Supplementary Fig. S1). The droplet translation persist for about 2-3 s after which the vortex pair recombines into the global vortex and the droplet restores its nearly-circular shape. The events of spontaneous splitting of the self-assembled vortex into two entities leading to a droplet elongation and subsequent bursts of self-propulsion are probabilistic and statistically rare compared to the regular behavior but the phenomenon is robust.

Conclusions
In this work, we experimentally explored the dynamics of motile colloids in soft confinement. We employ Quincke rollers in a droplet as a model system, with the rollers speed easily controlled by the applied field strength. We find that the interplay between deformable confinement and activity-driven flow gives rise to several previously unobserved phenomena. At low activity, the rollers form a vortex spanning the whole droplet, in contrast to rollers in solid-wall confinement, where the particles accumulate near the boundary. The droplet contour fluctuates about a circular shape and the fluctuations power spectrum is consistent with active fluctuations driven by particle-interface collisions. The non-equilibrium nature of the shape fluctuations is revealed by a broken detailed balance of the shape dynamics. While the interface deformation is driven by the particle-induced flow, the relaxation appears mainly controlled by surface tension as evidenced by the time correlations of the shape fluctuations. Shape fluctuations grow with the activity and can result in a spontaneous extension of the droplet in one direction driven by a formation of the vortex doublet. The spontaneous droplet elongations are accompanied by bursts of persistent self-propulsion in a direction perpendicular to the extension axis. The vortex splitting and recovery lasts for few seconds during which time the droplet can travel a distance of several droplet radii. The timing of the excursion and the direction of motion are randomly chosen. Their control requires understanding of the symmetry-breaking mechanism that leads to the vortex doublet formation. Our results provide insights into the complex dynamic behavior of active colloidal suspensions confined by a deformable boundary and provide new directions for future research in the engineering of selfpropelled micromachines. Inset: Linear dependence of the square average of velocity fluctuations hu 2 N i on ðE=E Q Þ 2 À 1. In both cases the line is a least squares fit to the experimental points. All experiments were performed on the same droplet with packing fraction φ = 0.58 ± 0.11 and droplet area A drop = 11.05 ± 0.03 mm 2 . b As predicted by Eq. (2), the energy density e is linearly dependent on φ. The line is a least squares fit to the experimental points. These experiments were performed at the same E = 4.25 ± 0.01 MV m −1 but different droplets with A drop in the range from 8.97 mm 2 to 11.25 mm 2 . c A graph of the power spectrum slope α, defined as 〈|h q | 2 〉 = αq −4 , versus e demonstrates a quadratic dependence (red line). The experimental points were obtained from several realizations of the active droplet: E from 3.19 to 4.52 MV m −1 ; φ from 0.13 to 0.58; and A drop from 8.97 to 11.25 mm 2 . Error bars are the standard error of the mean.

Methods
Experimental system. Spherical polystyrene particles (d = 2a = 40 μm diameter, Phosphorex) and density 1.06 g cm 3 in hexadecane (ρ hexadecane = 0.77 g cm 3 ) are dispersed in hexadecane with 0.1 mol L −1 dioctyl sulfosuccinate sodium (AOT) salt (Sigma Aldrich). A small amount of the solution with the particles (approx. 5 μL) is sandwiched between two indium tin oxide (ITO) coated glass slides to form a liquid bridge with a high aspect ratio. The distance between the electrodes is set by glass spheres with diameter 150 μm (Novum Glass) embedded in vacuum grease. Particles are allowed to sediment on the bottom electrode. We explore several packing fractions φ, ranging from 0.13 to 0.58, of Quincke rollers inside droplets. The packing fraction is determined by a thresholding method described in details in Supplementary Note II.
Imaging and droplet shape analysis. The recordings are made with a fast CMOS camera (IDT) at 300 fps and 500 fps mounted on a stereoscope (Leica). Velocity, vorticity fields and streamlines reflect the motion of the Quincke rollers and were obtained by a particle image velocimetry (PIV) package MatPIV for Matlab. Velocity fields together with the movement of the entire liquid bridge served as the input to calculate the energy per unit area E tot /A drop of the system.
To characterize the droplet contour fluctuations each image was automatically thresholded by Otsu's method in Matlab to obtain the border outline. The center of the droplet R c ! was determined as the mean coordinates of the droplet area pixels, A drop , and equivalent radius was calculated The perimeter of the equivalent circle with R drop around R c ! served as the baseline for the coordinate s and the radial deviations form the equivalent circle h(s) were determine from the images (Fig. 3a inset). We used the square of the fast Fourier transform algorithm in Matlab to compute the power spectrum of the fluctuations and averaged it over time. The number of frames to calculate the temporal averages of all the measured quantities was >2000, error bars represent the standard error of the mean value. The droplet radius of gyration R g and asphericity Δ were determined by the where index n runs over all area pixels of the liquid bridge drop (N in total). Q has eigenvalues λ 1 and λ 2 which define R g ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi λ 1 þ λ 2 p and Δ ¼ ðλ 1 À λ 2 Þ 2 =ðλ 1 þ λ 2 Þ 2 .
Detailed balance analysis of the shape fluctuations. The methodology follows ref. 52 to analyze the transitions between microscopic configurations defined as the shapes corresponding to different Fourier modes. We compute the trajectory in the phase space spanned by the two modes from the time series of their amplitudes. Then the phase space is discretized into equally sized, rectangular boxes each of which represents a discrete state. The probability is defined as the ratio of the time spent at a given state and the total trajectory time. The arrows indicate the currents across box boundaries determined by counting transitions between boxes. The transitions between neighboring discrete states occur when the system trajectory crosses box boundaries. Computing the contour integral of the probability current, Ω ¼ , shows a non-zero circulation for a system is out of equilibrium. Details of the methodology can be found in Supplementary Note III.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
All relevant code used in this study is available from the corresponding author upon reasonable request.
Received: 21 October 2021; Accepted: 22 March 2022;  Supplementary Fig. S2. Inset: A closer inspection reveals that the direction of motion is perpendicular to the maximal caliper diameter defined with an angle β (see inset snapshot, the green line is the obtained maximal caliper (Feret) diameter and the red line is the trajectory up to that moment, Scale bar is 1 mm). Experimental parameters: roller packing fraction φ = 0.28 ± 0.02, droplet area A drop = 10.07 ± 0.03 mm 2 , electric field strength E = 4.464 ± 0.008 MV m −1 .