Emergence of Bimodal Motility in Active Droplets

Artificial model swimmers offer a platform to explore the physical principles enabling biological complexity, for example, multigait motility: a strategy employed by many biomicroswimmers to explore and react to changes in their environment. Here, we report bimodal motility in autophoretic droplet swimmers, driven by characteristic interfacial flow patterns for each propulsive mode. We demonstrate a dynamical transition from quasiballistic to bimodal chaotic propulsion by controlling the viscosity of the environment. To elucidate the physical mechanism of this transition, we simultaneously visualize hydrodynamic and chemical fields and interpret these observations by quantitative comparison to established advection-diffusion models. We show that, with increasing viscosity, higher hydrodynamic modes become excitable and the droplet recurrently switches between two dominant modes due to interactions with the self-generated chemical gradients. This type of self-interaction promotes self-avoiding walks mimicking examples of efficient spatial exploration strategies observed in nature. DOI: 10.1103/PhysRevX.11.011043


I. INTRODUCTION
In response to physical constraints in nature, microorganisms have adapted and developed various locomotion strategies. Depending on cues from the environment, these strategies range from the more commonplace helical swimming [1,2], run-and-tumble, and switch-and-flick motility [3] to more sophisticated transient behaviors, e.g., peritrichous bacteria switching poles in response to a steric stress [4], octoflagellate microalgae exhibiting run-stop-shock motility with enhanced mechanosensitivity [5], and starfish larvae maximizing fluid mixing, and thereby nutrition uptake, through rapid changes of ciliary beating patterns [6]. Such intricate gait-switching dynamics [7,8] enable organisms to navigate in external flows [9,10], to follow gradients [11], or to efficiently explore their environment [12,13]. Recent efforts in the development of synthetic swimmers have led to the synthesis of systems that are capable of mimicking some of the aforementioned features of their natural counterparts such as rheotaxis [14,15], chemotaxis [16,17], and gravitaxis [18]. However, dynamic multimodal motility in the absence of external actuation has not been explored before in artificial swimmers, and the mechanisms underlying unsteady behavior in self-actuating systems are not well understood, particularly with respect to distinguishing states with increased random fluctuation from ones featuring true multimodal behavior [19,20].
Paradigms for biomimetic artificial swimmers include autophoretic microswimmers, powered by chemical activity at their interface, which are able to generate long-living chemical gradients in the environment [17]. In this regard, droplet microswimmers driven by micellar solubilization [21] provide a sophisticated experimental realization. Unlike most synthetic swimmers which are inherently asymmetric, active droplets are isotropic. Interfacial activity spontaneously breaks the symmetry, allowing for the emergence of different flow patterns depending on the environmental parameters. Here, we use such active droplets as model systems to demonstrate the physical principles guiding the emergence of multimodal motility in response to changes in environmental conditions.
We show that active droplets adapt to an increase in the viscosity of the swimming medium by exhibiting increasingly chaotic motion-a counterintuitive response given that increasing viscous stress generally tends to stabilize noninertial dynamics. Using time-resolved in situ visualization of the chemical and the hydrodynamic fields around the droplet interface, we find that the emergence of the chaotic dynamics correlates with the onset of higher hydrodynamic modes at increasing Péclet number Pe. Once these higher modes prevail, the droplet exhibits an unsteady bimodal exploration of space triggered by its interaction with a self-generated, slowly decaying chemical gradient. The conditions for the onset of this dynamical transition are quantitatively predicted by an advection-diffusion model for the transport of the chemical species, which takes into account the nonlinear coupling between the hydrodynamic and chemical fields. The visualization technique and the findings presented here lay the groundwork for future investigations of emergent dynamics in active phoretic matter.

II. DROPLETS PROPELLED BY MICELLAR SOLUBILIZATION
Our experiments use a specific subclass of active droplets: oil droplets that are slowly dissolving in supramicellar aqueous solutions of ionic surfactants. The droplets spontaneously develop self-sustaining gradients in interfacial surfactant coverage, resulting in Marangoni stresses which lead to self-propulsion [22]. This interfacial instability may be understood as follows [Figs. 1(a) and 1(b)]: During the solubilization of the droplet, oil molecules migrate into surfactant micelles in a boundary layer around the droplet interface, causing the micelles to swell and take up additional surfactant monomers from the aqueous phase, therefore reducing the local density of monomers c below its equilibrium, the critical micelle concentration (CMC). Unless there are empty micelles present to restore the CMC by disintegration, this local mismatch reduces the interfacial surfactant coverage, such that the interfacial tension increases with the local ratio of filled to empty micelles. Following an advective perturbation in the vicinity of the droplet, the initially radially isotropic cloud of filled micelles is distorted; the resulting fore-aft asymmetry generates a surface tension gradient toward the trailing oil-filled micelles which drives the droplet forward. Because of this self-sustaining gradient, the droplet propels continuously while leaving behind a trail of swollen micelles [ Fig. 1 As proposed by hydrodynamic theory models [23][24][25][26][27], such spontaneous self-propulsion stemming from an advection-diffusion-driven interfacial instability arises only if the Péclet number Pe, which characterizes the ratio of advective to diffusive transport, exceeds a critical threshold. In a simplified description, the surfactant dynamics are approximated by treating the droplet interface as a sink for surfactant monomers [23][24][25]. In this framework, on which we base our subsequent mode stability analysis, Pe can be shown to be a monotonically increasing function of the swimming medium (outer) viscosity μ o , here nondimensionalized as μ ¼ μ o =μ i using the constant inner viscosity μ i [see Appendix B 2 for a step-by-step derivation of Eq. (1)]: where V t is the theoretical terminal droplet velocity in a surfactant gradient [25,28] D ¼ ðk B T=6πr s μ o Þ the diffusion coefficient for the surfactant monomer (length scale r s ∼ 10 −10 m), q s the isotropic interfacial surfactant consumption rate per area, and ζ ∼ 10 nm the characteristic length scale over which the surfactants interact with the droplet [24,28]. Increasing μ o corresponds to an increase in Pe, besides the increase in viscous stresses. Henceforth, we represent an increase in μ o by this corresponding increase in Pe, as tabulated by the color map in Fig. 2. We note that, in view of the necessary simplifications in the derivation of Eq. (1), all experimental Pe values should be regarded as approximate.
In experiments, we control μ o via water-glycerol mixtures as the swimming medium (viscosity values in Appendix A, Fig. 7), varying the glycerol content from 0 to 70 vol % and adding the surfactant tetradecyltrimethylammonium bromide (TTAB) at 5 wt % to generate activity. Monodisperse CB15 oil droplets of radius R d ¼ 30 μm are studied in quasi-2D reservoirs 60 μm in depth for 5-6 min, a time over which the droplet radius should not shrink by more than approximately 7%. Therefore, we do not consider any droplet size effects on Pe (see also the experimental materials and methods section in Appendix A).

III. SIMULTANEOUS VISUALIZATION OF CHEMICAL AND HYDRODYNAMIC FIELDS
To visualize the chemical and hydrodynamic fields involved in the droplet activity, we directly image the chemical field of swollen micelles by adding the hydrophobic dye Nile Red to the oil phase [Figs. 1(c) and 1(d); see also Appendix A 5 and Supplemental Video S1 in [29] ]. The dye comigrates with the oil molecules into the filled micelles, which fluoresce when illuminated. We seed the surrounding medium, a supramicellar aqueous surfactant solution, with green fluorescent tracer colloids and measure the flow field using particle image velocimetry (PIV). The emission spectra of dye and colloids are  [29]. The color bar relating experimental Pe estimates to the viscosity ratio μ ¼ μ o =μ i applies to all subsequent figures. sufficiently nonoverlapping to be separately detected in dual-channel fluorescence microscopy. Consequently, both fields can be simultaneously observed and analyzed; we provide an example micrograph with an overlay of the extracted droplet trajectory in Fig. 1(d). Because of the large size (approximately 5 nm) of the filled micelles, the timescale of their diffusive relaxation exceeds that of the droplet motion; thus, there is a persistent fluorescent trail in the wake of the droplet.

IV. DESTABILIZED MOTION WITH INCREASING PÉCLET NUMBER
We begin, however, with an overview of the droplet dynamics using trajectory plots and statistical analyses of speed and orientational persistence taken from brightfield microscopy (Fig. 2). With increasing Pe, the droplet propulsion changes from uniform speeds and persistent motion to unsteady motion with abrupt reorientations [Figs. 2(a)-2(d)]. We define PðjδθðtÞjÞ as the distribution of the reorientation angle δθ of the 2D droplet velocity VðtÞ during a fixed time step δt [30]: PðjδθðtÞjÞ broadens significantly, corresponding to more frequent and sharper reorientation events [ Fig. 2(e)]. The faster decay of the angular velocity autocorrelation function,

V. SIGNATURES OF UNSTEADY DYNAMICS IN THE TIME EVOLUTION OF CHEMICAL AND HYDRODYNAMIC FIELDS
To investigate the origin of this unsteady behavior, we study the evolution of chemical and hydrodynamic fields around the droplet. We extract the tangential flow velocity u θ ðθÞ and the red fluorescence intensity IðθÞ of the chemical field close to the interface [ Fig. 3(d) and Appendix A 6] and map them in kymographs Iðθ; tÞ and u θ ðθ; tÞ, respectively.
For low Pe ≈ 4, at persistent propulsion, Iðθ; tÞ shows a single fixed-orientation band marking the origin of the filled micelle trail at the rear stagnation point of the droplet [ Fig. 3(a) and Supplemental Video S6 in [29] ]. The two bands in u θ ðθ; tÞ correspond to a steady flow field with dipolar symmetry that is consistent with the Iðθ; tÞ profile. On the right side of Fig. 3(a), we superimpose the streamlines of this dipolar flow field on the corresponding chemical micrograph at the time marked by I in the Iðθ; tÞ kymograph.
For intermediate Pe ≈ 36 [ Fig. 3(b) and Supplemental Video S7 in [29] ], Iðθ; tÞ shows secondary branches forming at the anterior stagnation point of the droplet and subsequently merging with the main filled micelle trail. This observation coincides with a transient second hydrodynamic mode with quadrupolar symmetry [ Fig. 3(b), II], causing the accumulation of an additional aggregate of filled micelles at the droplet anterior (see also Appendix C, Fig. 11 for additional flow field examples).
The ratio of the diffusive ðR 2 d =D fm Þ to advective ðR d =VÞ timescales for the migration of filled micelles is ðVR d =D fm Þ ≫ 1 for all experiments, assuming a diffusion coefficient D fm ¼ k B T=6πμ o r fm , with a micellar radius of Oðr fm Þ ∼ 2.5 nm. Therefore, the aggregate is unlikely to dissipate by diffusion and will continue to grow as long as the quadrupolar mode exists. However, this mode is not stable. Eventually, the dipolar mode dominates and advects the secondary aggregate toward the main trail [ Fig. 3(b), III]. The transport of the aggregate along one side of the droplet locally disturbs the interfacial flow, leading to an abrupt reorientation of the swimming direction [ Fig. 3(a), I-III]. As shown in the trajectories in Figs. 2(b) and 2(c), these reorientation events become more frequent with increasing Pe; accordingly, u θ in Fig. 3(b) exhibits quasiperiodic reorientation patterns.
For high Pe ≈ 293 [ Fig. 3(c) and Supplemental Video S8 in [29] ], the quadrupolar mode eventually prevails, resulting in a predominantly symmetric extensile flow around the droplet [ Fig. 3(c), I], as shown by a pronounced fourfold pattern in the additional kymograph u r ðθ; tÞ of the radial velocity. Because of the nonpropelling quadrupolar mode, the droplet is trapped in place. The gradual accumulation of filled micelles at the two stagnation points with radially outward flow manifests in two stable branches in the chemical kymograph [marked by I in Fig. 3(c)]. The growth of the two micellar aggregates locally generates a lateral chemical gradient, which eventually pushes the droplet out of its self-made trap. Concomitantly, the two points of filled micelle emission move along the droplet interface and merge on the new rear side of the droplet into a single filled micelle trail [ Fig. 3(c), II and III]. The chemorepulsion from the local field micelle gradient induces an apparent dipolar mode which gradually decays as the droplet leaves the self-made trap. Now, the quadrupolar mode resaturates, with an aggregate growing at the droplet anterior, until the droplet is trapped again and a new bimodal "stop-and-go" cycle begins. Since the escape direction is always lateral, consecutive runs are approximately perpendicular, resulting in the sharp reorientation events apparent in the trajectories in Figs. 3(c) and 2(d), as well as the broadening jδθj distribution in Fig. 2(e).

VI. DEPENDENCE OF HYDRODYNAMIC MODES ON THE PÉCLET NUMBER
In order to understand the dependence of the onset of bimodal motility on Pe, we analyze the underlying advection-diffusion problem for the active droplet within the framework of an axisymmetric Stokes flow as established in Refs. [23,25,26,32] (see Fig. 4 and Appendix B). At the smallest value of μ, Pe is approximately equal to the critical value of 4 necessary for the onset of the first hydrodynamic mode (n ¼ 1), i.e., the mode with dipolar flow symmetry [23,25,26]. With increasing μ, Pe [markers in Fig. 4(a)] eventually exceeds the critical values necessary for the onset of the higher hydrodynamic modes [lines in Fig. 4(a)], specifically the second hydrodynamic mode (n ¼ 2), i.e., the mode with quadrupolar symmetry. A linear stability analysis around an isotropic, quiescent base state (see Appendix B 3 and Refs. [23,26]), which is the idealized starting point for each experiment, shows that, for small to moderate Pe, the nondimensionalized instability growth rate λ for n ¼ 1 exceeds that for n ¼ 2 Each frame corresponds to the point in time indicated on the kymographs by I, II, or III. (d) defines the mapping of the profiles of red light intensity I (filled micelle concentration) and tangential velocity u θ around the droplet circumference onto the y axis of the kymographs in the middle column. All u θ profiles are in the translational droplet reference frame but with θ ¼ 0 fixed at the laboratory x direction to visualize the reorientation dynamics. In (c), the third kymograph corresponds to the radial velocity u r in the laboratory reference frame to better depict the quadrupolar symmetry of the flow field. The second hydrodynamic mode starts to appear at intermediate Pe and dominates the dynamics for high Pe. See also Appendix C, Fig. 11 Fig. 3(c)], it appears that the n ¼ 2 mode can also evolve from a nonquiescent state and prevail in a similar Péclet regime as derived from the performed stability analysis. Note that we restrict our analysis to the first two hydrodynamic modes, since these two are solely responsible for the droplet propulsion and the associated far-field hydrodynamic disturbance.

VII. INTERACTIONS WITH SELF-GENERATED CHEMICAL GRADIENTS CAUSE SPEED BURSTS
It remains to explain the broadening of PðVÞ with increasing Pe [ Fig. 2(e)], particularly the remarkable bursts in speed for high Pe. While the dipolar mode is propulsive, the quadrupolar mode is not. Hence, the growth and decay of the respective modes affect the droplet speed. As shown in Fig. 3, recurrent transitions between the two hydrodynamic modes lead to abrupt reorientation events; we therefore investigate the correlation between changes in speed and reorientation angle jδθj.
In a typical trajectory for intermediate Pe ≈ 36, each sharp turn is preceded by a deceleration and followed by an acceleration, as shown in the plot of the positional data color coded by speed in Fig. 5(c). Signatures of these correlations in the droplet dynamics appear in the conditional averages of jδθj, V, and tangential acceleration a t as quantities X for all sharp reorientation events i in the trajectory, centered ; the events are identified by choosing a threshold value of jδθj > 0.2 (see Appendix C, Fig. 12).
We can now directly compare these dynamics to the higher-resolution fluorescence data taken at Pe ≈ 36 presented in the kymographs in Fig. 3(b). Figure 5(b) shows a series of micrographs of the chemical field, with arrows marking the droplet velocity vector (black) and the position of the secondary filled micelle aggregate (white). The aggregate accumulates, is then entrained, and finally merges with the posterior trail, corresponding to the creation and merging of a secondary chemical branch in the kymograph.
For t < 0, the droplet decelerates while the secondary aggregate is accumulating. t ¼ 0 marks the point in time where V is minimal and the aggregate is on the cusp of leaving the anterior stagnation point. For t > 0, the aggregate is advected to the droplet posterior, and the droplet accelerates due to the resaturation of the dipolar mode. V peaks once the aggregate has merged with the main trail-creating an amplified fore-aft gradient-at t ≈ 1 s, which is comparable to the advective timescale R d =V ≈ 1 s. In the wide-field data analysis in Fig. 5(a), this is the time τ 1 it takes the droplet to reach maximum speed after a reorientation.
We now use the correlation function between V and jδθj, C jδθj;V ðΔtÞ ¼ hjδθðtÞj · Vðt þ ΔtÞi t , plotted in Fig. 5(d), to estimate the growth times of the second mode from our data for Pe > 10. Since V is minimal at maximum jδθðtÞj approximately corresponds to the timescale for the growth and resaturation of the n ¼ 2 mode during the bimodal motility (i.e., starting from a nonquiescent base state). Nevertheless, we compare this experimentally obtained τ 2 − τ 1 with the theoretical growth times for the n ¼ 2 mode starting from the isotropic base state, λ −1 n¼2 R d =V t [ Fig. 4(c)], for different values of Pe. Figure 5(e) shows that these two timescales, which are strictly speaking different, still are of the same order of magnitude and show a similar decreasing trend with increasing Pe. We note that the growth time of the dipolar flow above Pe ≈ 100 cannot be used for comparison to λ n¼1 , since this flow is imposed by the lateral chemical gradient. However, we can assume that this gradient increases with Pe, resulting in faster acceleration, markedly higher swimming speeds, and, hence, reduced τ 1 , as observed experimentally [ Fig. 5(d)].

VIII. CONSEQUENCES FOR SPATIAL EXPLORATION
Reminiscent of gait-switching dynamics in biological locomotion, we demonstrate the emergence of complex swimming behavior in a minimal active droplet system by tuning the Péclet number. We find a transition from persistent swimming at low Pe to chaotic bimodal swimming at high Pe-the latter results from the excitation of higher hydrodynamic modes beyond critical Pe values, while the continuous switching between them is caused by the self-generated chemical gradient in the environment.
This gradient sensitivity causes trail avoidance [16], which, in turn, affects the way these droplet swimmers explore their environment. With increasing reorientation frequency, we find a transition from quasiballistic propulsion to a 2D self-avoiding walk (2D SAW). This effect is illustrated by the trajectories in Figs. 2(a)-2(d) and also by the fact that C VV in Fig. 2(e) does not decay to zero. For a statistical analysis, we plot mean squared displacements for selected Pe values in Fig. 6(a), which reproduce the expected scaling with t 2 (ballistic) for Pe ≈ 4 and a transition to t 3=2 (2D SAW [33]) for Pe ≳ 36, with the crossover time decreasing with increasing Pe. While transitions to random walks governed by run-and-tumble gait switching are common in bioswimmers [34], selfavoidance requires chemical self-interaction [35].
Examples of anomalous diffusion driven by repulsive biochemical signaling are found in the spreading of slime FIG. 5. Interactions with self-generated chemical gradients cause speed bursts at reorientation events. (a) Conditional averaging of tangential acceleration a t , speed V, and reorientation angle jδθj, for abrupt reorientation events at Pe ≈ 36 (see Appendix C, Fig. 12, for an illustration of the identification criteria). The dotted line marks the maximum speed at t ¼ τ 1 after reorientation. (b) Video stills of the chemical field for one such event with t ¼ 0 s set to the point of minimum speed; white arrows track the accumulation of the secondary filled micelle aggregate at the anterior stagnation point and its advection along the interface, and black arrows correspond to the droplet velocity vector. The droplet speed is maximal when the secondary aggregate and the trail merge at t ¼ 0.93 s. See also Supplemental Videos S9 and S10 in [29]. (c) An example trajectory for Pe ≈ 36. Any reorientation event (curved arrows) is preceded by a deceleration and followed by an acceleration. The lowest speed occurs at the point with the highest curvature. (d) Correlation function between reorientation angle and speed, C jδθj;V ðΔtÞ for increasing Pe. Times τ 1 and τ 2 (next reorientation event) are identified by the respective peak and dip in C jδθj;V . (e) Timescale for the growth of the n ¼ 2 mode versus corresponding Pe: experimentally obtained, τ 2 − τ 1 (∘), compared to values from stability analysis, λ −1 n¼2 R d =V t (□).

IX. CONCLUSION
In this work, we demonstrate that the manner in which hydrodynamic and self-generated chemical fields are coupled determines the nonlinear dynamics of autophoretic microswimmers. The fluorescence-based visualization technique used to simultaneously probe this coupling can provide insight into many recent autophoretic models [20,21,23,24,[38][39][40]. For example, extensive theoretical studies [41][42][43][44] demonstrate the importance of quantifying far-field and near-field contributions, coupling to chemical fields, and the effects of confinement to understand how swimmers approach each other or form bound states, which is vital to nutrient entrainment, food uptake, and mating in bioswimmers.
While many microswimmer models incorporate unsteady dynamics via stochastic fluctuations, we show that the interplay of nonlinear dynamics and interaction with the history of motion also allows for the emergence of memory-driven chaotic behavior. An appealing example from a different field are droplet walkers on a vibrated bath [45], which show a transition from persistent to a bimodal, stop-and-go motion based on an effective "system memory" parameter [46,47]. The corresponding theoretical framework [46] is general enough to also apply to bimodal chaotic motion in droplet swimmers.

Materials and characterization
Our samples consist of droplets of (S)-4-cyano-4'-(2methylbutyl)biphenyl (CB15) doped with the fluorescent dye Nile Red in an aqueous solution of the cationic surfactant TTAB corresponding to 5 wt % (50 mg in 1 ml of solution) in pure water, with a critical micelle concentration of CMC ¼ 0.13 wt %. We purchase CB15, TTAB, and Nile Red from commercial suppliers (Synthon Chemicals and Sigma-Aldrich) and use them as is. We control the viscosity of the swimming medium μ o by adding glycerol to the aqueous TTAB solution.
We use an Anton Paar MCR 502 rotational rheometer to characterize the shear viscosity of water-glycerol-surfactant solutions (Fig. 7). Experiments are carried out using a coneplate geometry, to find shear-rate versus shear-stress curves at a fixed temperature and viscosity versus temperature at a fixed shear rate. To limit effects of solution evaporation, the coneplate geometry is surrounded by a water bath and covered by a Peltier hood. Over the shear rate range 0.01 s −1 < _ γ < 100 s −1 , viscosity is found to be constant, such that our solutions are well described as Newtonian, as should be expected: Water-glycerol mixtures are used as Newtonian standard media throughout the existing literature.
To estimate the surfactant consumption rate q s in Eq. (1), we extract the droplet shrinking rate dR d =dt from the bright-field microvideography data presented in Fig. 2. We find a moderate dependence on the glycerol fraction (Fig. 8), which we include as a first-order approximation, via linear regression (blue line), to evaluate q s in the Pe estimates in the main manuscript.

PDMS soft lithography for droplet generation
For the production of monodisperse oil droplets, we fabricate microfluidic channels in house, using standard soft lithography techniques. First, 2D photomasks are designed in AutoCad and then printed onto an emulsion film in high-resolution (128 000 dpi) by a commercial supplier (JD Photo-Tools). Next, the photoresist SU-8 3025 (MicroChem) is spin coated onto a 4-inch-diameter silicon wafer (Si-Mat), where spin speed and duration are adjusted to give a controllable uniform thickness. A negative mold is cured in the SU-8 through the photomask by UV light exposure. After further chemical treatment with photoresist developer, uncured SU-8 is removed, leaving behind cured SU-8 microstructures on the silicon wafer.
We then pour a poly(dimethyl siloxane) (PDMS, Sylgard 184, Dow Corning) mixture of 10∶1 volumetric ratio of base to cross-linker over the wafer and bake for 2 h at 80°C, producing a solid PDMS layer with microstructured indentations. We peel the indented PDMS from the wafer and punch holes through it to create liquid inlets and outlets at opposing ends of the channels. The structured PDMS surface, as well as a glass cover slip, are cleaned and treated with partial pressure air plasma (Pico P100-8; Diener Electronic GmbH+Co. KG) for 30 s and then pressed together, bonding the two surfaces. Figure 9 shows a micrograph of such a PDMS chip during droplet production.
The walls of these microfluidic chips are selectively treated to hydrophilize the channels where surfactant solution will flow. This treatment prevents oil from wetting the walls during droplet production. We follow the technique of Petit et al. [48]: First, the channel walls are oxidized by a 1∶1 mixture of hydrogen peroxide solution (H2O2 at 30 wt %, Sigma-Aldrich) and hydrochloric acid (HCl at 37 wt %, Sigma-Aldrich). This mixture is flushed through the channels for approximately 2 min by using a vacuum pump system. After the oxidation, the channel is rinsed by flushing double distilled water for 30 s. Next, a 5 wt % solution of the positive polyelectrolyte poly(diallyldimethylammonium chloride) (PDADMAC, Sigma-Aldrich) is flushed for 2 min through the oxidized channel of the device. The PDADMAC binds to the activated channel walls by ionic interactions. Finally, a 2 wt % solution of the negative polyelectrolyte poly(sodium 4-styrenesulfonate) (PSS, Sigma-Aldrich) is flushed for 2 min.

Droplet generation
Once the chips are treated, we mount syringes of oil and 0.1 wt % aqueous TTAB solution to a microprecision syringe pump (NEM-B101-02B; Cetoni GmbH), connect these to the two inlets of the microfluidic chip via Teflon tubing (39241; Novodirect GmbH), and tune the flow speed through the chip until the desired droplet size is reached. Once droplet production is monodisperse (after approximately 5 min) and at a steady state, these droplets are collected in a bath of 0.1 wt % TTAB solution. This solution is of a high enough concentration to stabilize the droplets against coalescence but not high enough to induce solubilization.

Fabrication of the observation Hele-Shaw cell
The swimming behavior of the droplets is observed in a quasi-2D Hele-Shaw reservoir, which we fabricate directly from SU-8 photoresist without PDMS casting. To fabricate the reservoirs, we therefore use a photomask with inverted polarity. We spin coat the photoresist directly onto a glass slide (50 × 75 mm 2 ) and follow the same procedure for photolithography as outlined in Appendix A 2. This process results in a layer of cross-linked SU-8 (thickness approximately 60 μm) with reservoirs of the dimensions 8 × 13 mm. These reservoirs are filled with the samples, sealed with a glass cover slip, and put under a microscope.

Double-channel fluorescent microscopy technique
We use double-channel fluorescent microscopy for simultaneous imaging of the chemical and hydrodynamic fields. A schematic of the setup is shown in Fig. 10 diameter), which visualize the fluid flow around the droplet, are excited with a 488 nm laser and emit light at a maximum of approximately 510 nm. The emitted light is separated using a beam splitter and appropriate filters for each emission maximum. We also use a spatial pinhole (confocal microscopy) to enhance image quality. Examples of snapshots recorded on each channel are shown in Figs. 10(b) and 10(c).

Image processing and data analysis
To observe the long-time statistical behavior of the active droplets, as in Fig. 2, we observe their motion in a glassbounded Hele-Shaw cell (quasi-two-dimensional reservoir, 13 × 8 mm and height h ≈ 60 μm) under a bright-field microscope (Leica DM4000 B) at low magnification (5×) compared to the double-channel fluorescence microscopy setup. Videos are recorded at a frame rate of ten frames per second using a Canon (EOS 600d) digital camera (1920 × 1080 px). The droplet coordinates in each frame are extracted from video frames using the common Python libraries numpy, PIL, and openCV (scripts available on request). Steps include background correction, binarization, blob detection by contour analysis, and minimum enclosing circle fits. Swimming trajectories are obtained using a frame-by-frame nearest-neighbor analysis.
To acquire the kymographs of the chemical field and tangential and radial velocities around the droplet interface, we observe the droplet behavior by double-channel fluorescent microscopy as described in Appendix A 5. We use a 512 × 512 pixel camera at a frame rate of 14 frames per second connected to a 20× objective. First, we split the red (Nile Red, filled micelles) and green (tracer particles) channels. Then, the red frames are used to extract the droplet coordinates via the blob detection algorithm described above. We use a MATLAB script that centers the droplet and records the red light intensity value along the interface at a distance 15.6 μm for Pe ≈ 4 and 36 and 20.4 μm for Pe ≈ 293. We note that it is not possible to record the intensity closer to the interface, because the strong fluorescence from the large amounts of dye inside the droplet create a very bright region extending several micrometers beyond the actual interface. We plot the extracted profiles versus time to generate spatiotemporal kymographs.
For a quantitative analysis of the flow field around the droplet, we perform particle image velocimetry (PIV) on the tracer particles images (green channel) using the MATLAB-based PIVlab interface [49]. The objective is focused on the midplane of the Hele-Shaw cell. We define a moving mask for the area covered by the droplet. We perform the analysis in 16 × 16 pixel interrogation windows with 75% overlap. The spatial resolution is 1.2 μm=px. After obtaining the velocity vector field, we center the droplet and read the velocity vectors at a certain distance from the droplet interface (3.6 μm for Pe ≈ 4 and 36 and 8.4 μm for Pe ≈ 293). The tangential (u θ , in the droplet reference frame) and radial (u r , only for Pe ≈ 293, in the lab reference frame) velocity components are then calculated and plotted in the kymographs. Because of the impermeability boundary condition, the radial component of the velocity directly at the interface is supposed to be zero; however, since we read the values at a certain distance from the interface, there is an inward and outward radial contribution to the flow. We use this observation, in particular, at Pe ≈ 293 to show the quadrupolar symmetry of the flow field at the stopping moment.
In Fig. 1(a) and Supplemental Video S1 in [29], we track the droplet and center it in the image. To obtain the path lines of the tracer particles in the video, we use FlowTrace [50] to convolve a superposition of ten frames for each image. For Fig. 1(a), we superimpose 30 frames. To visualize the motion of the tracer particles in Fig. 4(b), IV, and the Supplemental Videos S6-S9 in [29], we process the green channel of the input video (8-bit RGB) as follows: For each pixel coordinate, the intensity is replaced by its standard deviation within a 20-frame window around the current frame. Each frame is subsequently contrast maximized within a [0, 255] intensity range. The red and blue channels are not modified. This procedure is inspired by ImageJ's Z projection algorithm; the respective Python code is available on request.

APPENDIX B: VISCOSITY DEPENDENCE OF HYDRODYNAMIC MODES
In this appendix, we describe the mathematical framework for the coupled hydrodynamic and advectiondiffusion problems pertaining to the active droplet system. Note that we follow the solution methodology of Refs. [23,25,26] and rework each step of the analysis for the present system. The appendix shows the origins of all expressions and equations (including the scaling analyses necessary for simplifications) needed to understand the theoretical framework and, importantly, the origin of Fig. 4. We especially show each step of the linear stability analysis so that the derivation of the equations governing the instability growth rates for the hydrodynamic modes is clear.

Governing equations and boundary conditions
for the active droplet system Considering an axisymmetric Stokes flow (Reynolds number for the swimming of the active droplet Re ∼ 10 −4 ) and the impermeability of the droplet interface, the flow field around and inside the spherical active droplet (capillary number Ca ≪ 1) can be expressed in terms of the nondimensional stream function ψ, in ðr; θÞ coordinate system, as [25,26,32] Here, and in the subsequent discussions, superscripts o and i refer to quantities outside and inside the active droplet, respectively, r is the radial coordinate nondimensionalized by droplet radius R d , η ¼ cos θ, and P n ðηÞ is the Legendre polynomial of degree n with the prime denoting its derivative; n here physically represents the nth hydrodynamic mode. The nondimensional radial and tangential flow velocity components around and inside the droplet are related to ψ as u r ¼ −ð1=r 2 Þð∂ψ=∂ηÞ and u θ ¼ −½1=rð1 − η 2 Þ 1=2 ð∂ψ=∂rÞ. The coefficients a n and b n in Eqs. (B1) and (B2) are constrained by the following boundary conditions [25,32].
The coefficients on the right-hand side of Eqs. (B3) and (B4) essentially stem from the nondimensionalization of the classical boundary conditions. Note that the flow velocity is nondimensionalized using V t ¼ ½q s ðγ c R d þ 3μ i MÞ= Dð2μ o þ 3μ i Þ, which is a theoretical estimate for the terminal velocity of the active droplet considering the contributions of both the Marangoni and the diffusiophoretic effects [25,28]. Furthermore, μ ¼ μ o =μ i is the ratio of the swimming medium viscosity μ o to the droplet viscosity μ i , and the nondimensional parameter m represents the relative strengths of diffusiophoretic to Marangoni effects [25]. Essentially, m can be considered as a ratio of the diffusiophoretic velocity scale to the viscocapillary velocity scale representing the Marangoni effect. Accordingly, is the diffusiophoretic mobility [24,28], γ c ≈ k B Tζ is the leading-order change in the interfacial surface tension γ with surfactant concentration c (alternatively, γ c ¼ ðdγ=dcÞ can be considered to be a measure of the change in γ with c assuming a linear variation) [24,25], and ζ ∼ 10 nm is the characteristic length scale over which the surfactants interact with the droplet in the interfacial region. For the active droplet system, OðmÞ ∼ 10 −3 − 10 −2 for the entire range of experiments; hence, for the present physical problem, the diffusiophoretic effect is much weaker as compared to the Marangoni effect. However, the former is considered in the analysis here for the sake of generality. In the definition of V t , q s is an isotropic and constant interfacial surfactant consumption rate per unit area necessary for the droplet activity, and D ¼ ðk B T=6πr s μ o Þ is the diffusion coefficient for the surfactant monomer (length scale for surfactant monomer r s ∼ 10 −10 m). q s can be approximately estimated by assuming that the total number of surfactant monomers necessary per unit time to account for the volumetric reduction rate of the droplet due to the formation of the filled micelles is equal to the total interfacial surfactant consumption rate. Hence, jdV d =dtjN s =v fm ≈ q s 4πR 2 d , which implies that q s ≈ ð3N s jdR d =dtjÞ=ð4πr 3 fm Þ. Here, OðN s Þ ∼ 25 is the number of surfactant monomers per filled micelle, v fm ¼ 4=3πr 3 fm is the filled micelle volume with a micellar radius of Oðr fm Þ ∼ 2.5 nm, and jdR d =dtj is the droplet solubilization rate as given in Fig. 8.
Equations (B3) and (B4) delineate the dependence of the swimming hydrodynamics on the distribution of the nondimensional surfactant concentration c in the vicinity of the droplet. Naturally, c is governed by an advection-diffusion relation [25,26,32]: The distribution of c is subject to the following boundary conditions: (i) isotropic and constant surfactant consumption at the droplet interface (r ¼ 1). ∂c ∂r (ii) the bulk condition.
Note that Eq. (B6) addresses the depletion of the interfacial surfactant monomers due to the creation of the filled micelles by considering the isotropic and constant interfacial surfactant adsorption rate per unit area of q s , corresponding to a flux with unit of number per area per time (in dimensional form, D∇c Ã ·n ¼ q s ; this estimate gives a scale for the surfactant concentration as ∼ðq s R d =DÞ) [25,26]. Pe in Eq. (B5) is the system Péclet number-the details of which are discussed in the following subsection. The above system of equations [Eqs. (B1)-(B7)] can be solved for ψ (therefore, u r and u θ ) and c using the singular perturbation technique for certain limiting cases [25,26]. The solvability condition clearly shows that the actuations of different hydrodynamic modes depend on certain threshold values of Pe [ Fig. 4(a) in the main text] [25]. Furthermore, the asymptotic analysis also provides a physical understanding of the hydrodynamic and surfactant concentration fields corresponding to the different modes, specifically n ¼ 1 and n ¼ 2 [ Fig. 4(b) in the main text].

The system Péclet number
The important thing to understand now is the dependence of Pe on μ. Classically, Pe can be written as Pe ¼ ðV t R d =DÞ, where V t ¼½q s ðγ c R d þ3μ i MÞ=Dð2μ o þ3μ i Þ is the theoretical estimate for the terminal velocity of the active droplet considering the contributions of both the Marangoni and diffusiophoretic effects, as mentioned in the preceding subsection [25,28]. Utilizing the aforementioned definition of V t and following some simple algebraic manipulations, Pe can be expressed in terms of system constants and the parameter μ as In the last step of Eq. (B8), the approximate expressions for M and m (see Appendix B 1) and the definition of D (see Appendix B 1) are utilized to derive the final expression for Pe. Equation (B8) expresses Pe as a monotonically increasing function of the viscosity ratio μ [markers in Fig. 4(a) in the main text]. Note that q s is approximately estimated by relating the dissolution rate of the active droplet to the isotropic and constant surfactant consumption at the droplet interface [24]; the dissolution rate of the active droplet is dependent on the glycerol concentration (Fig. 8), which effectively makes q s dependent on μ o . We further note that the second term in the numerator within parentheses Oðζ=R d Þ ∼ 10 −4 ; this small magnitude further substantiates the fact that the diffusiophoretic effect is much weaker compared to the Marangoni effect for the present system.

Linear stability analysis about a motionless (isotropic) base state
For the linear stability analysis (also see Refs. [23,26]), the time-dependent form of the advection-diffusion equation [Eq. (B5)] is used: Next, the desired quantities are expressed in terms of the unsteady (instability) modes-ψ ¼ e λt P nψ n ðrÞP n ðηÞ and c ¼ −ð1=rÞ þ e λt P ncn ðrÞP n ðηÞ, where λð>0Þ is the nondimensional growth rate for the instability modes. Using the aforementioned expressions for ψ and c and linearizing Eq. (B9), the governing equations for the first two modes can be obtained as d dr d dr where λ s ¼ ffiffiffiffiffiffiffi λPe p and a 1 and a 2 are the coefficients of the first and second modes, respectively, of the outer stream function [as in Eq. (B1)]. Equations (B10) and (B11) are solved to evaluatec 1 andc 2 , respectively: Here, x ¼ rλ s is a rescaled spatial variable, ChiðxÞ and ShiðxÞ are the hyperbolic cosine integral and hyperbolic sine integral functions, respectively, and α 1 and α 2 are the constants of integration. Note that Eqs. (B12)  Using Eqs. (B12) and (B14), α 1 can be evaluated as Similarly, using Eqs. (B13) and (B14), α 2 can be evaluated as Considering the hydrodynamic boundary conditions [Eqs. (B3) and (B4)] and using the orthogonality condition for Legendre polynomials, a set of two simple algebraic equations for the coefficients a n and b n for each of the first two modes can be written as (i) first mode (n ¼ 1): (ii) second mode (n ¼ 2): Note thatc n in the above equations is explicitly dependent on a n [see Eqs. Equations (B21) and (B22) are solved numerically to evaluate the variations of the nondimensional growth rates ½λ ¼ ðλ 2 s =PeÞ with Pe for the first and second instability modes, respectively [Fig. 4(c) in the main text]. Note that Eq. (B21) is identical to that derived for the spontaneous motion of an autophoretic isotropic particle [23]. Furthermore, it is important to note here that the inverse of the timescale used for nondimensionalizing the growth rate is V t =R d , which is consistent with the entire analysis. chemical field kymographs plotted for a longer period of 60 s. Supporting Videos S6-S8, respectively, correspond to the kymographs in Figs. 11(a)-11(c).
In Fig. 12, we plot the long-time tangential acceleration, speed, and reorientation angle for Pe ¼ 36. This dataset is used to identify the abrupt reorientation events. We identify these events based on a cutoff criterion for the reorientation between video frames jδθj¼0.2 rad [Figs. 12(c) and 12(d)], aligned and overlaid the profiles of all events with the turning point (jδθ max j) set as t ¼ 0, and calculated the timedependent average (hi represents ensemble averaging over all events).
In Fig. 13, we plot the long-time acceleration signal for Pe ¼ 293 to demonstrate signatures of bimodal swimming. Such events can be identified by intermittent strong fluctuations in the acceleration profile. The enlarged view further demonstrates the difference between stopping (n ¼ 2) and swimming modes (n ¼ 1). Constant transitions between these modes result in the anomalous diffusive behavior shown in Fig. 6 in the main text.