Anomalous dynamics of a passive droplet in active turbulence

Motion of a passive deformable object in an active environment serves as a representative of both in-vivo systems such as intracellular particle motion in Acanthamoeba castellanii, or in-vitro systems such as suspension of beads inside dense swarms of Escherichia coli. Theoretical modeling of such systems is challenging due to the requirement of well resolved hydrodynamics which can explore the spatiotemporal correlations around the suspended passive object in the active fluid. We address this critical lack of understanding using coupled hydrodynamic equations for nematic liquid crystals with finite active stress to model the active bath, and a suspended nematic droplet with zero activity. The droplet undergoes deformation fluctuations and its movement shows periods of “runs” and “stays”. At relatively low interfacial tension, the droplet begins to break and mix with the outer active bath. We establish that the motion of the droplet is influenced by the interplay of spatial correlations of the flow and the size of the droplet. The mean square displacement shows a transition from ballistic to normal diffusion which depends on the droplet size. We discuss this transition in relation to spatiotemporal scales associated with velocity correlations of the active bath and the droplet.

The random motion of a passive tracer particle moving in a fluid at equilibrium is the starting point of describing fluctuating dynamics [1].Microrheology suggests that the local and bulk mechanical properties of a complex fluid can be extracted from the dynamics of such tracer particles [2,3].A similar strategy can be used for an active fluid, where particle movement is coupled to nonequilibrium degrees of freedom.An active fluid consists of a collection of particles or cells suspended in a fluid, which can convert chemical energy into mechanical work by generating stresses at the microscale [4][5][6][7].For instance, a flock in a fluid can be considered as a liquid crystal which is made of the swimming particles and an oriented state of the system possesses an 'active stress' which is proportional to the orientation [8].Active fluids have been shown to exhibit unusual transport properties such as enhanced diffusion, accumulation near boundaries, and rectification [9,10].Experiments using passive particles in active fluids have revealed superdiffusive behavior of a passive particle at short times and a dramatically enhanced translational diffusion at long times [11][12][13][14][15][16][17][18][19][20].
The collective dynamics of the suspended swimmers in the fluid such as microbial suspensions, cytoskeletal suspensions, self-propelled colloids, and cell tissues, generates spontaneous flows.These flows are characterized by chaotic spatiotemporal patterns and are referred to as active turbulence [21,22].Experimental studies of active turbulence look at the velocity fields and their statistics to characterize long range hydrodynamic flows [23] and show the creation, transportation and annihilation of topological defects by active stresses [24].Theoretical studies of passive tracers in active fluids use the generalized Langevin equation with different properties of the friction and fluctuating force due to the fluid or use particle based simulations of passive and active objects [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43].To the best of our knowledge, numerical studies of passive particle dynamics in active fluids have almost always neglected hydrodynamic interactions of the swimmers because of the complexity involved [44] or have focused on the drag on the passive particle in the presence of applied force [35].Recent experiments of mobile rigid inclusions in an active nematic system have shown the importance of the complex interplay and feedback between the motion of such an inclusion and the active nematic [45].In this letter we computationally study the dynamics of a mobile deformable inclusion, which we call a droplet, in an active turbulent fluid.
One of the theoretical approaches to study active turbulence is active liquid crystal theory where the continuum equations of motion are built from symmetry arguments.A particular class is active apolar nematic systems characterized by head-tail symmetry of the constituent particles -tending to align with each other.These systems exhibit characteristics such as instability of the isotropic nematic state and intriguing pattern formation [2,6,21,[49][50][51], universal scalings of the energy spectra [52][53][54], topological chaos [55], and in general far from equilibrium dynamics [56].The homogeneous isotropic state is unstable and goes into a sequence of instabilities [57].Once the system is in a developed state, the ensemble averaged mean kinetic energy in the system approches a statistically steady state [3] [Fig. 1 (e)].We suspend a passive nematic droplet at this stage of the active turbulence and study its dynamics [Fig. 1 (ad)].Inverse to this setup are systems consisting of active droplets immersed in passive media which demonstrate complex behaviors, for example, spontaneous onset of motility and division of active nematic droplets [59], selforganization and division in active liquid droplets [60], dynamic defect structures in active nematic shells [61], and emulsification in mixtures of active and passive component media [62].Average kinetic energy per unit mass with time.Gray symbols for individual realizations and the red line is ensemble average.Time from t = 0 to t ≈ 300, when the active turbulence is still developing, is neglected and a time window of 300 ≤ t ≤ 800 is considered for all the statistical results in this paper, including the above subfigures (a-d) which are taken at t = 500.(f) Time averaged energy spectrum with spectral exponent close to −4, a behavior also observed by [46,47].(g) Individual realizations of trajectories of the geometric center r(t) of passive nematic droplet of radius R = 15 and R = 45, and (h) droplet interface during typical realizations for R = 15 and R = 45.(i) Position of the droplet center relative to initial position, |r(t)| with time, exhibiting "run" and "stay" chracteristics along the droplet trajectories.(j) Ratio of the horizontal to vertical span of the droplet for R = 15, 30 and 45.Dynamics of translations as well as the aspect ratio slows down as the radius is increased.(k) Probability distribution function of the translational steps taken by the geometric center of the droplet.PDF(∆r) narrows down upon increasing the radius of the droplet.(l) PDF of discrete changes in the direction of motion of the droplet.PDF(∆θ) scales nearly exponentially for intermediate range of angles (∆θ ≈ 25 o to 65 o ) with heavy tailed rare events for > 90 o or close to complete reversal.
We computationally characterize a system where the environment or the bath is an active liquid crystal and a passive nematic soft object is the suspended phase.It is known that an isotropic liquid droplet in a passive nematic liquid crystal under applied shear undergo different modes of movement (oscillatory, breakup, or motile) when the anchoring conditions at the surface of the droplet are changed [63].However inside an active turbulent medium, the forcing is far from simple shear, or random.The forces on the droplet are correlated over certain spatial and temporal scales and are an important aspect of its dynamics.
Apolar active systems in the continuum limit can be modeled using the passive nematic liquid crystal theory with added active stresses, combined with hydrodynamics.The orientation of nematics is described by a director field n, with n ≡ −n for apolar nematics.The symmetric traceless second rank tensor Q = q d d−1 (nn − I/d) characterizes the nematic order, where q is the strength of the nematic order, I is the identity tensor and d denotes the spatial dimension.In our simulations, d = 2.The evolution of Q is described by the nematodynamic equation [2,50] where and Ω = [∇u−(∇u) T ]/2 being the strain rate and vorticity tensors respectively.The tensor H describes relaxation of Q towards the minimum of a free energy and can be written as where K is the elastic constant and C is a material constant [2].The Q-tensor field is coupled to velocity field u via the momentum equation where the passive elastic stress tensor . The stress tensor σ A = −ζQ gives rise to activity or motility in the system [64].The activity parameter ζ ≥ 0 (≤ 0) for a contractile (extensile) active nematic.The term −µu mimics substrate friction or drag, and p is the pressure.The fluid is considered incompressible so that ∇ • u = 0. We suspend a droplet as a second phase inside the active medium.The force f I is due to the interfacial tension and acts only at the interface between passive droplet and the outer active medium [1,66,67].The coupled set of nematodynamic, momentum, and incompressibility equations are integrated computationally, utilizing a finite-volume based pressure projection algorithm, utilizing the Eulerian one-fluid approach [1,67,68].In this, the properties and parameters such as density, viscosity, activity, alignment, friction, etc. are also considered as fields and take different values in the two phases (see supplemental for parameter values in individual phases).The spatial distributions of the properties is updated with time using the front-tracking algorithm [1,66,67] in which the interface location is first updated using the velocity field solution and then the property fields and forces due to interfacial tension are reconstructed from the known interface location [1,66,67].One feature of the present formulation is that the anchoring conditions at the surface of the suspended droplet emerge from the solution itself, instead of implementing them as an input.The advective fluxes in the continuum equations are reconstructed using a sixth order weighted essentially non-oscillatory scheme which is quite volume conservative [69,70].
The system of equations is unstable around the state Q = 0, u = 0; the system passes through states consisting of vortices and then proceeds towards spatiotemporal chaos.Eventually active stresses, on an average, are balanced by the passive, dissipative and frictional stresses and the system reaches a statistically steady state [Fig. 1 (a-d)].To characterize the active turbulent state, we look at the kinetic energy spectrum where k is the magnitude of the wavevector or the wavenumber, û(k, t) = u(r, t) exp(−ik.r)dr,and denotes time average.We find a power law scaling E(k) ∼ k −4 [Fig. 1 (f)] which has also been observed by [46,47].
Droplets of radii in the range 15 ≤ R ≤ 45 are suspended in the active turbulent state, and for a given radius the simulations are repeated at least 15 times to obtain an ensemble of stochastic trajectories.A typical set of trajectories of the geometric center of the droplet is shown in Fig. 1 (g-h) for radius R = 15 and R = 45.To be precise, the position of the center of the droplet is defined as r(t) = i x i (t)/N where the perimeter of the droplet is divided into N points having positions x i -consistent with the front-tracking algorithm [1,66,67].As the droplet interacts and is forced by the surrounding active bath, we measure its distance relative to the initial position with time |r(t)|; typical cases are depicted in Fig. 1 (i).Along the trajectories, the droplets undergo periods of "runs" and "stays" as marked on the plot.The "stays" are prolonged as the radius of the droplet is increased.We also note that the time scale of undulations of droplet deformation parameter A(t), defined as the ratio of horizontal to vertical span of the droplet, increases as we increase the radius [Fig. 1 (j)].With increasing size the droplet dynamics slows down and we transit from a "fast" process towards a "slow" process.
The trajectories consist of discrete translational steps taken by the droplet center ∆r ≡ |r(t + ∆t) − r(t)| in time step ∆t.In addition to translational steps ∆r, we define angular steps ∆θ = cos πA/(4τ α 3/2 ) with discretization time step τ .This implies a diffusive mean square displacement (MSD) ∼ Dt for t τ .However, the computation of MSD as shows a superdiffusive regime where ∆r 2 (t) ∼ t δ with power exponent 1.5 ≤ δ ≤ 2.0, before any transition to normal diffusion (see Fig. 2 (ab)).This is consistent with the behaviour observed in experiments of tracer diffusion in bacterial turbulence.
To further understand the droplet dynamics, we compute the normalized step autocorrelation function C ∆r (τ ) = ∆r(t) ∆r(t + τ ) / (∆r(t)) 2 [Fig.2(d of the background velocity field of the active bath, and associated integral length scale, v ≡ ∞ 0 C vv (r)dr (see Fig. 2 (e-f)).We find that as we decrease the radius of the droplet below v , the droplet executes relatively larger ∆r [Fig. 2 (c) and Fig. 1 (k)].Conversely droplets with size larger than v undergo negligible overall displacement in a given time.The variance of step size distribution drops as the size of the droplet is increased, and if R > v the variance drops to ≈ 1 [Fig. 2 (c)].We propose that droplets with R v experience relatively larger drifts and droplets with R v move relatively smaller distances because in the latter case the velocities at any two opposite ends of the droplet are expected to be decorrelated -resulting in negligible net advection of the interface.It is also to be noted that although the change in radius alters the integral time scale τ ∆r associated with the step autocorrelation function [Fig. 2 (d, inset)], it has negligible influence on v [Fig. 2 (e, inset)] implying that suspending a droplet has negligible effect on the velocity characteristics of the background active bath.
Our observations suggest that the random forces due to the active bath on the droplet are correlated and depends on the droplet size.Therefore, we assume the following form for the random force with a memory which decays exponentially in time [11,71], where τ F is a characteristic time scale, and Λ ≡ Dγ 2 /τ F , where γ is the coefficient of friction from the Langevin equation mr(t) = −γ ṙ(t) + F (t). Mean square displacement r(t) 2 can be obtained (see Supplementary for detailed derivations) and is depicted in Fig. 2 (c, inset).The crossover from a ballistic r(t) 2 ∼ t 2 to diffusive r(t) 2 ∼ t is observed.In the limit t τ F , the velocity autocorrelation is given as which gives the correlation time τ ∆r = [mk B T /τ F Λ] ≡ mk B T /Dγ 2 .As we increase R, the mass and inertial delay time increases, but is not the only factor behind the change in the dynamics as from simulations we see that τ ∆r decreases as we increase R (and m) of the droplet.This is because in the limit t τ F the velocity autocorrelation function have the form [11] v(t)v(0) = (2D/τ F ) exp [−t/τ F ].In this limit, the correlations in the random active force dominate over the inertial effect.Additionally, the value of R relative to v can have significant effect on the change in the crossover time.Smaller particles in active bath take too long to pass the crossover to diffusive MSD regime because τ ∆r increases -as shown through in Fig. 2 (c, inset).
Our results can be experimentally tested, for example using a setup of quasi two dimensional suspension of beads and bacteria on a soap film [11] or inside thin fluid films [72].Experiments on self-diffusion of tagged bacterium in a dense swarms grown on agar plates [73], or of amoeboid cells [74], already confirm existence of superdiffusive regime while it is hard to reach normal diffusive regime.For example, multi-potent progenitor cells exhibit superdiffusive regime which can be persistent even for hours [75].Here we have pin pointed how the size of the droplet and the properties of the active bath changes the spatial and temporal scales [11,[72][73][74][75].Our numerical study will also add to the better understanding of proposed analytical theories of diffusion of bodies in active systems [13,25,[38][39][40][41][42][43].
We acknowledge support from High Performance Computing facility at IISER Mohali.CS acknowledges support from the INSPIRE Faculty Fellowship of the Department of Science and Technology, India.

SUPPLEMENTAL MATERIAL: DYNAMICS OF A PASSIVE DROPLET IN ACTIVE TURBULENCE RESCALING OF GOVERNING EQUATIONS AND SIMULATION PARAMETERS
If we substitute the following reduced or rescaled variables/fields/properties/operators into the dimensional form of the governing equations and if we fix the following scales then the dimensional form of the governing equations can be rescaled into the following non-dimensional form Re ρ where denotes a non-dimensional variable/field/property/operator, and are non-dimensional groups signifying the ratio of inertial force to viscous force, active force to viscous force, friction force to viscous force, and interfacial tension force to viscous force, respectively.Implementation of the interfacial force f I is described in detail in [1].Numerical values of the non-dimensional parameters and non-dimensional properties in the active and passive phases are summarized in Table I.In addition to fluid properties in the active and passive phases, the table also summarizes other simulation parameters namely the system size, grid size, time step, and droplet radius.The grid and time steps are taken such that the advective as well as viscous time step conditions are well satisfied, namely [1] ∆t Numerical values of parameters and fluid properties adopted in the simulations.The above set of parameters are close to, for instance, the studies carried out by [2,3].
where ∆ is the grid step L x /M = L y /N , and subscripts max, min denote maximum and minimum values of the repective properties from the two phases.For simulation to remain stable, the factors C 1 , C 2 , C 3 have to be < 1, and to be conservative, we choose time and grid size combination such that these conditions are always satisfied.

MEAN SQUARE DISPLACEMENT AND VELOCITY AUTOCORRELATION FUNCTION
To incorporate the memory effect in the stochastic force in the Langevin equation, let us start from an exponentially correlated random force where τ F is the time scale which tells us the extent of the force memory, and Λ ≡ Dγ 2 /τ F with D and γ being diffusion constant and coefficient of friction respectively.For simplicity of the argument, only a component of the force is considered and the formulation can readily be extended to higher dimensions.The equation of motion for a passive particle in such an active bath can be written as Multiplying the equation of motion with x and v ≡ d dt x respectively, and taking the ensemble average, after few algebraic manipulations we obtain the equations for the mean square displacement (MSD) and the mean squared velocity, written as While the correlation function r(t)F (t) = 0 (force on the droplet by the active bath is independent of the position of the droplet), the correlation function v(t)F (t) is in general non-zero and carries information about the balance between the energy supplied by the active bath to the droplet and the energy dissipated by the droplet due to the drag offered by the bath -the fluctuation-dissipation mechanism.This correlation function can be obtained by multiplying F (t) on both sides of the trivial identity , and taking the ensemble average The second term on right hand side is zero (velocity at some earlier time has to be independent of the force at a later time), and we input d dt v(t ) from the equation of motion to obtain The first term on right hand side is again zero (velocity at some earlier time has to be independent of the force at a later time).Taking the limit s → ∞, and using exponentially correlated force from Eq. 16, we get Using Eq. 22 in Eq. 19, and assuming statistically stationary state for v(t) 2 ≡ k B T /m, we find that Here T is an effective temperature of the active bath.Using this in Eq. 18, we find the equation for the MSD for an exponentially correlated noise So we find that the time scale τ F of the exponentially correlated random force changes the crossover time of transition from ballistic to diffusive regime.It is important to note that there appears no such time scale in equation for MSD in case of delta correlated random force.Force time scale τ F can be related to velocity correlation time scale τ ∆r as follows.Let us multiply the Langevin equation with v(0) and then take the ensemble average, we get The last term is zero by the same argument that velocity at some earlier time has to be independent of the force at a later time.Thus we have This provides us an estimate of the correlation time τ ∆r associated with the velocity autocorrelation function, i.e.
as Λ ≡ Dγ 2 /τ F .The above result is valid in the limit t τ F , and it can be proved that in the limit t τ .

SPECTRAL ANALYSIS AND DISCUSSION
Figure 3 show the power spectrum of signals: (a) translational step size ∆r(t), (b) change in direction of motion ∆θ(t), and (c) droplet aspect ratio A(t), respectively.The indication from the analysis is that the translational step and the angular change in the direction of motion signals exhibit nearly 1/ω or pink noise.The spectrum of the aspect ratio or the deformation parameter A(t), however, exhibits nearly 1/ω 3 noise.Although its known that 1/ω noise is present in many natural, man made, and scientific systems [4][5][6][7], its underpinnings are debated over decades [8][9][10].Does above 1/ω noise has to do with self-organized states in our system?In addition, 1/ω 3 noise in the deformation or the droplet aspect ratio signal indicates to us that deformations somehow have distinct dynamics than the motion increments.This type of noise is usually observed in oscillators implying that the droplet deformation has some oscillatory feature.Figure 3  thus the case is different from standard Brownian increments as Brownian increments are independent.Similarly Fig. 3 (b) imlies that two consecutive ∆θ are not independent.

FIG. 1 .
FIG. 1. (a-b) Absolute velocity scaled by its root mean square value for system size 200 × 200 and 400 × 400 respectively.The droplet is shown with the solid yellow line.(c-d) Vorticity field for system size 200 × 200 and 400 × 400 respectively.(e)Average kinetic energy per unit mass with time.Gray symbols for individual realizations and the red line is ensemble average.Time from t = 0 to t ≈ 300, when the active turbulence is still developing, is neglected and a time window of 300 ≤ t ≤ 800 is considered for all the statistical results in this paper, including the above subfigures (a-d) which are taken at t = 500.(f) Time averaged energy spectrum with spectral exponent close to −4, a behavior also observed by[46,47].(g) Individual realizations of trajectories of the geometric center r(t) of passive nematic droplet of radius R = 15 and R = 45, and (h) droplet interface during typical realizations for R = 15 and R = 45.(i) Position of the droplet center relative to initial position, |r(t)| with time, exhibiting "run" and "stay" chracteristics along the droplet trajectories.(j) Ratio of the horizontal to vertical span of the droplet for R = 15, 30 and 45.Dynamics of translations as well as the aspect ratio slows down as the radius is increased.(k) Probability distribution function of the translational steps taken by the geometric center of the droplet.PDF(∆r) narrows down upon increasing the radius of the droplet.(l) PDF of discrete changes in the direction of motion of the droplet.PDF(∆θ) scales nearly exponentially for intermediate range of angles (∆θ ≈ 25 o to 65 o ) with heavy tailed rare events for > 90 o or close to complete reversal.
−1 [a • b/ab] where a = r(t + ∆t) − r(t) and b = r(t) − r(t − ∆t).Distribution of these two quantities are shown in Fig. 1 (k-l).The step length distribution (or equivalently the step velocity distribution) spreads upon decreasing the radius of the droplet [Fig. 1 (k)].Smaller droplets execute relatively larger translational steps.The change in direction of motion of the droplet is nearly exponentially distributed for intermediate range of angles (∆θ ≈ 25 o to 65 o ), and exhibits some rare events of > 90 o or close to complete reversal of the direction [Fig. 1 (h,l)].Assuming a functional form for the translational distribution f (∆r) = A e −α(∆r) 2 , the translational diffusion can be calculated as D = ∞ −∞ (∆r) 2 f (∆r)d(∆r)/2τ = √

FIG. 2 .00
FIG. 2. (a) Mean square displacement of droplets of radius R = 15, 30 and 45 with to = 0, and (b) zoom into the super diffusive regime with to = 300.(c) Variance of discrete step sizes taken by droplets of different radii.The variance decreases as the droplet radius increases; (c, inset) Mean square displacement shows crossover from ballistic to diffusive regime for a particle under exponentially correlated random force with time scale τF .The crossover depends on τF , or the velocity correlation time τ∆r ∼ m.(d) Step autocorrelation function C∆r(τ ) and (d, inset) associated integral time scale τ∆r ≡ ∞ 0 C∆r(τ )dτ .Smaller droplets exhibit relatively larger τ∆r, or increased memory of the history of steps taken.(e) Time averaged one-time two-point velocity correlation function of the velocity field in active nematic environment.The integral length scale v ≡ ∞ 0 Cvv(r)dr is marked with a vertical dashed line; (e, inset) v in the active bath changes negligibly upon changing the droplet radius.(f) Integral length scale with time during typical realizations.

aFIG. 3 .
Figure3show the power spectrum of signals: (a) translational step size ∆r(t), (b) change in direction of motion ∆θ(t), and (c) droplet aspect ratio A(t), respectively.The indication from the analysis is that the translational step and the angular change in the direction of motion signals exhibit nearly 1/ω or pink noise.The spectrum of the aspect ratio or the deformation parameter A(t), however, exhibits nearly 1/ω 3 noise.Although its known that 1/ω noise is present in many natural, man made, and scientific systems[4][5][6][7], its underpinnings are debated over decades[8][9][10].Does above 1/ω noise has to do with self-organized states in our system?In addition, 1/ω 3 noise in the deformation or the droplet aspect ratio signal indicates to us that deformations somehow have distinct dynamics than the motion increments.This type of noise is usually observed in oscillators implying that the droplet deformation has some oscillatory feature.Figure3(a) implies that two consecutive droplet position increments are not independent and