Oscillatory active microrheology of active suspensions

Using the method of Brownian dynamics, we investigate the dynamic properties of a 2d suspension of active disks at high Péclet numbers using active microrheology. In our simulations the tracer particle is driven either by a constant or an oscillatory external force. In the first case, we find that the mobility of the tracer initially appreciably decreases with the external force and then becomes approximately constant for larger forces. For an oscillatory driving force we find that the dynamic mobility shows a quite complex behavior—it displays a highly nonlinear behavior on both the amplitude and frequency of the driving force. In the range of forces studied, we do not observe a linear regime. This result is important because it reveals that a phenomenological description of tracer motion in active media in terms of a simple linear stochastic equation even with a memory-mobility kernel is not appropriate, in the general case.

In the past three decades microrheology [1][2][3][4] has emerged as a useful technique for characterizing rheological properties of complex fluids, i.e. fluids that incorporate mesoscopic building units such as colloids, polymers or more complicated self-assembled structures. Unlike traditional rheology 5 which quantifies the viscoelastic properties of complex fluids by connecting stresses and (rates of) strain in the medium, in a typical microrheological [6][7][8][9][10][11][12] measurement one follows the trajectory of a probe particle (tracer) embedded in the medium under consideration. A microrheological study can be conducted in two ways: passive and active. Passive microrheology, in which one tracks the random motion of the tracer induced by fluctuations in the medium, has been employed to study a broad assortment of complex systems, spanning from colloidal suspensions [13][14][15] to biological matter 16,17 . In active microrheology, one either measures the mobility of the tracer particle driven by a static force [18][19][20] , or probes the frequency dependence of mobility by subjecting the tracer to an oscillatory external driving force 21 . Additional information about the material's response can be acquired by utilizing rotational microrheology of anisotropic probes 22 .
In this article we investigate the active microrheological response of a low Reynolds number suspension of active disks exhibiting self-propulsive motion [23][24][25][26] . They expend energy to propel themselves forward, and therefore constitute a nonequilibrium suspension, which displays a plethora of intriguing properties. For example, active particles can exert an active pressure on confining surfaces [27][28][29][30][31] while getting stranded at them 32,33 , which can be exploited for constructing rotational [34][35][36] and translational [37][38][39][40] ratchet motors powered by active particles as well as for cargo transport 41,42 and active self-organization [43][44][45] . Despite these dazzling features of active suspensions, their microrheological properties are still not sufficiently explored [46][47][48][49][50][51][52] . Previous theoretical studies have mainly focused on passive microrheology [53][54][55][56] and in particular probing the fluctuation-dissipation relations in active media [57][58][59][60][61] , or on a constant force active microrheology, either by relying on the low-density description based on the Smoluchowski equation 62 , or on numerical simulations for a wide range of suspension densities 63 . Here we use Brownian dynamics simulations to study the active microrheological response of active suspensions. Firstly, we extend the results of reference 63 by subjecting the tracer to a broad range of external forces, and report new results on the resulting nonlinear dependence of the tracer mobility on the driving force. Secondly, for the first time in literature, we report an oscillatory active microrheological study of active suspensions. We reveal that the frequency-dependence of tracer mobility can be described by a Lorentzian form in the low frequency domain.
Our setup is depicted in Fig. 1: We study a tracer particle of radius R, immersed in a two dimensional suspension of active disks, and driven by an external force F e (t) . The active disks have a fixed self-propulsion speed v A , diameter σ , mobility µ A and they perform persistent motion within a characteristic time τ R = D −1 R , where D R is their rotational diffusion constant. They interact among themselves and with the tracer via purely repulsive steric forces. The "bare" mobility of the tracer in a fluid free from active disks is then µ T = µ A σ/2R . Neglecting hydrodynamic interactions, the motion of the tracer and active disks is described with a set of overdamped stochastic equations, which we solve numerically. The equations of motion and details of the numerical integration www.nature.com/scientificreports/ scheme are presented in the Methods section. The suspension of active disks is described by two dimensionless parameters: the Péclet number Pe = σ v A /D and the total area fraction occupied by disks φ = Nσ 2 π/(4A) ; here N is the number of disks, D is their translational diffusion constant, and A is the total area of the simulation box. The Péclet number compares the time t D = σ 2 /D it takes an active disk to diffuse its own length with the time t S = σ/v A it takes the disk to swim the same distance. For Pe ≫ 1 the active motion dominates over the diffusive transport, while for Pe ≪ 1 the disk behaves as a passive particle exhibiting diffusive motion. Throughout this study we set φ = 0.12 , so that motiltiy-induced phase separation [64][65][66] does not occur, and select several values of Pe ranging from the limit of suspensions of passive disks, Pe = 0 , to highly active disks, Pe = 240. The tracer is driven through the active medium with the external force F e (t) acting along the x-direction with the aim to probe the microrheological properties of the suspension. The collisions with active disks influence the mobility of the tracer in the direction of applied force: instead of having the bare mobility µ T the tracer now acquires an effective mobility µ.
Firstly, we examine the case in which a constant external force, F e (t) = F = const is applied on the tracer. We define the static mobility where v is the average speed of the tracer in the x-direction in the steady state. We introduce a dimensionless parameter f = σ F/(2Rf A ) , where f A = v A /µ A can be interpreted as the force to stop an active disk from moving. It is easy to see that the tracer driven by a force f = 1 , in a fluid free from disks, moves with a speed v = v A . It is expected therefore that the speed of the tracer moving in an active suspension should be much smaller than v A for f ≪ 1 , and conversely much larger than v A for f ≫ 1 . According to our results these two regions are separated by a relatively narrow crossover region located in the vicinity of f ≈ 1 . Indeed, we find that the mobility µ of the tracer in the active suspension displays a nonlinear dependence on f below this crossover value ( f 1 ), and it is amplified with respect to the mobility of a tracer driven through the bath of passive disks. On the other hand, for larger forces ( f 1 ) the tracer mobilities in passive and active baths match each other, and we find that they are independent of f in this region.
Secondly, we study a tracer driven by a harmonic external force F e (t) = F sin(ωt) , where ω is the frequency of oscillations, and F now denotes the force amplitude. We find that the velocity v =ṽ sin(ωt + φ) of the tracer, averaged over many independent realizations, oscillates with the same frequency ω as the driving force, with the phase angle φ ≈ 0 ; on the other hand, the velocity amplitude depends on both the amplitude and the frequency of the driving force, ṽ =ṽ(F, ω) . As a consequence of this the dynamic mobility, depends on ω and F as well. For small values of driving force amplitude ( f 1 ) we find that the dynamic mobility decreases rather quickly with ω up to some limiting value of ω and that our results fit well to a Lorentzian. The Lorentzian spreads out for larger values of f. It turns out, however, that for sufficiently large values of f (depending on setup conditions) the mobility µ changes its behavior-it is increasing rather than decreasing with ω . Demonstrating that the dynamic mobility significantly depends on frequency for small driving forces but also on the driving forces itself is the main result of our study.
(1) µ = �v�/F, www.nature.com/scientificreports/ Recently, it has become common in literature to describe the motion of a tracer in an active bath via a set of effective coarse-grained stochastic equations with a time-independent 37,40,51 effective friction coefficient. As our results suggest, it is hardly possible to get a satisfactory description in terms of linear stochastic equations in general. They also suggest that in the limit where a linear response is applicable one has to include a memorymobility kernel 51,56,67 instead of simply a constant mobility.
The rest of the article is organized as follows. The "Results" section is partitioned into two segments presenting constant force and oscillatory force microrheology of suspensions of active disks, respectively. We elaborate the implications of our oscillatory force results and offer our conclusions in the Discussion section. Finally, the Methods section presents equations of motion of the tracer and active disks and provides details of their numerical integration.

Results
In the first part of this section we present our results obtained for the tracer driven through a suspension of active disks by a constant external force, while the second part shows our findings for the tracer guided through the suspension by a harmonic force.

Constant force microrheology.
Taking the average of Eq. (8) from the Methods section, which describes the tracer dynamics, and using Eq. (1), we find that the effective mobility of the tracer driven by a constant force is where the sum goes over all active disks in the bath and F T ix is the x-component of the force exerted on the tracer by the i-th disk. We first perform an average over 100 independent simulation runs and then a time average over the steady-state motion of the tracer (see Methods section for details). This double averaging is denoted by . . . . In Fig. 2a (solid lines) we show the temporal fluctuations of the tracer speed v, obtained by averaging over independent simulation runs only, for several external force amplitudes f. One observes a clear increase in tracer speed v with increasing driving force f. For the same set of external forces f we also show (black dashed lines) the tracer speed v in a bath of passive disks ( Pe = 0 ). For values of driving force f 1 , including the narrow crossover region located around f ≈ 1 , we find that the tracer moves with a larger speed in an active than in a passive bath, and thus has a higher effective mobility, as clearly illustrated in Fig. 2b. On the other hand, for larger external forces, f 1 , the tracer moves with approximately equal speed in the passive and active baths. Thus, in this regime the tracer speed increases linearly with f and the mobility µ does not depend on f. In addition, we find that the mobility µ versus f data fall roughly on top of each other for baths with different but large activities Pe = 80 and Pe = 240 . This is due to using the rescaled external the characteristic force scale of the active bath. For the range of external forces considered, the tracer mobility in the passive bath does not change with f, as shown in Fig. 2b. This seems to be in agreement with the findings presented in Ref. 15 in the limit of low bath densities.
To develop some understanding for the observed behavior, in Fig. 3 we plot the average density of disks in the vicinity of the tracer driven through passive and active baths. In passive baths the tracer leaves behind a channel free from disks, which is qualitatively similar for all external forces f that we have considered. For f ≥ 0.2 the tracer moves so fast that the diffusing passive disks cannot immediately fill the space behind the tracer (see Supplementary Movie 1, displaying a tracer driven through a passive bath with f = 0.2 ). Thus, the tracer pushes the www.nature.com/scientificreports/ passive disks in front of it forward, which in turn reduces its mobility with respect to µ T , while leaving behind it a trace without disks (cf. Fig. 3, top row). We find the tracer speed to increase linearly with external driving f, leading to a constant mobility µ for the range of forces considered, Fig. 2b. In contrast to the passive bath, the tracer moving through the active bath with Pe = 80 displays a much richer behavior with changing f in the same range, as the bottom panels in Fig. 3 demonstrate. For large force f = 5 , the tracer moves faster than the active disks, v > v A , and consequently leaves a wake behind it; this is illustrated in Movie 2. The shape of the wake is somewhat different compared to the passive bath since the active disks move ballistically into the wake, which hence assumes the shape of a cone (compare top and bottom right panels of Fig. 3). Nevertheless, the motion of the tracer is qualitatively similar in these two cases. Indeed, the tracer moving with a large speed does not distinguish whether the disks it encounters on its front are passive or active, meaning it has the same effective mobility µ in both types of baths as demonstrated in Fig. 2b. As the driving force is gradually reduced towards f = 1 , by entering the crossover region the speed of the tracer becomes roughly equal to the active disk speed, v ≈ v A , and the active disks start to catch up with the tracer, allowing them to push it from behind. In the density profile of Fig. 3, middle bottom panel, this manifests itself by the disappearance of the void at the rear side of the tracer. For even lower external force, f = 0.2 , the active disks move significantly faster than the tracer, and accumulate more at the rear than at the front side of the tracer, left bottom panel in Fig. 3 and Movie 3. Thus, they also push the tracer forward. This explains the increase of the mobility with decreasing f and is investigated in more detail in the following.
The force i F T ix = F F + F R , with which active disks push against the tracer, can be split into two components F F and F R . They are exerted by the active disks along the x-direction on the front half and rear half of the tracer, respectively. We can rewrite Eq. (3) to get . Note that f R has the sign of the x-component of tracer's velocity, while f F has the opposite sign ( f R > 0 , f F < 0 ). In Fig. 4a we plot the absolute values |�f F �| and |�f R �| as functions of the external force f for an active bath with Pe = 80 . Obviously, the front force has a higher magnitude than the rear force, |�f F �| > |�f R �| . As can be inferred from Fig. 4a, for f 1 the force exerted on the tracer front scales linearly with f, while the force acting on the rear side of the tracer is vanishingly small, |�f R �| ≈ 0 . This behavior agrees well with our previous observation of constant tracer mobility µ for large external force (cf. Fig. 2b). One can also notice that |�f F �| further decreases with f but in a nonlinear fashion in the region f 1 . On the other hand, |�f R �| gradually increases with decreasing f in the same region, which corresponds to the onset of active disk accumulation at the rear-side of the tracer as observed in Fig. 3. Taken  Fig. 5. The time evolution of the velocity v(t) has been calculated by averaging over 100 independent simulation runs. The velocity profiles were then fitted to a simple form v(t) =ṽ sin(ωt + φ) , with the velocity amplitude ṽ and phase shift φ being the fit parameters. We have been able to get very good fits with φ ≈ 0 for all amplitudes f and frequencies ω examined in our study, which clearly has to be assigned to the active or nonequilibrium nature of the bath particles performing an overdamped motion. On the other hand, we  The variation of the dynamic mobility µ(ω) with frequency ω in a suspension having Pe = 80 is calculated using Eq. (2) and presented in Fig. 6. Depending on the magnitude of the amplitude f, we distinguish two qualitatively different regimes. In the case of small and moderate amplitudes ( f 7 for the case Pe = 80 ) we find that µ(ω) decreases with ω up to approximately ωτ R = 3π . In particular, in this region µ(ω) changes quite quickly for the case of small amplitudes ( f = 0.2, 0.5 ), but the change slows down as f increases gradually ( f = 1.5, 2.5, 5, 7 ). Nevertheless, in this region all of our mobility data sets fit quite well to the shifted Lorentz form: µ ∞ + µ 0 /[1 + (ωτ ) 2 ] , where the fit parameters µ ∞ , µ 0 and τ depend on the amplitude f. We find that τ −1 , which provides a measure of the characteristic width of the Lorentz form, increases with f (see Fig. 6). For larger frequencies ( ωτ R > 3π ) it seems that µ slowly increases with ω , somewhat faster for larger values of f.
As one can also notice from Fig. 6, the overall shape of the curve µ(ω) changes considerably by increasing f (compare e.g. the cases f = 0.2 and f = 7 ). It would be interesting therefore to explore the regime in which the magnitude of the driving force significantly exceeds the self-propulsion force of active disks. It is expected, namely, that in the limit of large force ( f ≫ 7 for Pe = 80 ), the dynamic mobility µ(ω) should be alike to that one finds in the passive bath (note that in the passive bath µ(ω) is an increasing rather than decreasing function of ω , as illustrated in Fig. 7). It is however rather difficult to perform numerical simulations in this limit because in that case one has to deal with a very small simulation time step t . Instead, one can consider a tracer driven by a force of smaller amplitude f, but immersed in a bath having a lower Pe number. For Pe = 30 , two examples are given in Fig. 7 for amplitudes f = 0.5 and f = 2.5 . One can see that in the case f = 2.5 , the mobility µ(ω) is a nondecreasing function of ω . For f = 0.5 the mobility decreases in the low frequency domain, and after that, in the region of larger frequencies µ(ω) increases. However, interestingly the variation in the probed frequency range is rather small. For the sake of comparison, in Fig. 7 we also present the mobility for the bath with Pe = 80 , as well as for the case of a passive bath ( Pe = 0 ). In contrast to the case of active bath,the mobility of the tracer in a passive bath jumps up quickly to a plateau value µ = µ T in the low frequency region. This can be explained by the fact that a tracer oscillating with a higher frequency moves practically in a region without obstacles, since the time t R = R 2 /4D = 12τ R it takes a passive disk to diffuse the tracer radius R is larger than the period T = 2π/ω of tracer oscillations. Conversely, when t R < T the disks have enough time to diffuse into the wake created by the moving tracer, and thus they act as obstacles, which leads to a decrease of the mobility. Then the crossover region between these two regimes can be roughly estimated by using the condition t R ≈ T , which gives ωτ R ≈ π/6 . As one can verify, this estimate indeed falls into the crossover region, see Fig. 7. www.nature.com/scientificreports/

Discussion
In this article we studied the microrheological properties of 2d suspensions of active disks. For this purpose we considered the motion of a tracer particle immersed in the active suspension. The tracer was driven either by a constant or an oscillatory force. In the case of a constant driving force, we found that the mobility of the tracer significantly decreases with growing magnitude of the external force f for weaker strengths ( f 1 ), while it approaches an approximately constant value for stronger forces ( f 1 ) [see Fig. 2b], which agrees with the tracer mobility in a suspension of passive particles. We stress that we did not observe a plateau in the mobility in the region of small f values, meaning that a linear response approach is not applicable in our system. Thus, the active suspension exhibits a highly nonlinear behavior for weak forces, f 1 . As a consequence of this, a simple coarse-grained description 37,40,46,51 of tracer dynamics in active suspensions in terms of an effective linear stochastic equation is not applicable in our case. It is interesting to note that in the region of very small forces ( f 0.01 ) a mobility plateau has been observed in active suspensions in the low density limit 62 . This is in contrast with our observations for active suspensions, at least for Pe ≥ 80 and the range of f values we explored so far. We note that much more numerical efforts are needed to get a reliable estimate of the mobility for very small values of f, due to very large fluctuations of the tracer velocity. For this reason the existence of a region in which the system displays a linear response remains open for now.
Importantly, we find that the dynamic mobility of a tracer driven by an oscillatory force is a quite complex quantity due to the ability of the active disks to push against the tracer. The mobility shows a strong frequency dependence up to moderate force amplitudes ( f 2.5 in the case of Pe = 80 ). In the range of frequencies ωτ R ≤ 3π the mobility decreases with ω and its frequency dependence can be well described by a Lorentzian that spreads out with f, as presented in Fig. 6. Neglecting the weak increase of the mobility observed at high frequencies, this means that the tracer motion can be described by an effective Langevin equation with an exponential memory-mobility kernel but only provided that the system's response is linear for very small f. Our aim in the future is to provide a more detailed numerical study of the system in this region. Let us note that in our study we fixed the ratio 2R/σ = 8 . It is possible, however, that a linear regime might appear more clearly by changing the ratio 2R/σ or some other bath parameter.
In Ref. 47 the active microrheology experiments performed in a bacterial bath do not show any frequency dependence in the dynamic mobility. We attribute this to a bacterial density of φ = 3 × 10 −3 , which is much smaller than φ = 0.12 used in our work. Thus, in this article we demonstrate that at moderate densities the dynamic mobility depends on frequency and shows a highly nonlinear behavior in the driving force. With our simulation work we provide a clear orientation where to search for this behavior in the available parameter space when performing active microrheology experiments of active suspensions of artifical or biological microswimmers.

Methods
We study a system of N interacting active disks in two dimensions which move with a constant speed v A and have mobility µ A . Their translational and rotational dynamics are described by a set of coupled overdamped Langevin equations Here r i denotes the position vector and n i = (cos θ i , sin θ i ) the unit orientation vector of the i-th disk. The disks are subjected to Gaussian noises of zero mean, �ξ i (t)� = 0 and �η i (t)� = 0 , and unit variance, www.nature.com/scientificreports/ �ξ α i (t)ξ β j (t)� = δ αβ δ ij δ(t − t ′ ) and �η i (t)η j (t ′ )� = δ ij δ(t − t ′ ) , and have translational diffusivity D and rotational diffusivity D R . They interact among themselves via purely repulsive pairwise forces F A ij = −∇ r i V σ (r i − r j ) which stem from the Weeks-Chandler-Andersen (WCA) potential The potential introduces a characteristic distance σ , where it assumes the strength ε . One can roughly interpret this distance as the diameter of the active disks. The force −F T i is connected to the passive tracer, which we introduce now.
A passive tracer disk of radius R and mobility µ T = µ A σ/2R is immersed in the active bath and driven by an external force F e (t) = F e (t)e x . We investigate two different settings: a time-independent driving force F e (t) = F and a harmonic force F e (t) = F sin(ωt) , where ω is the oscillation frequency. The position vector r T of the tracer changes according to The i-th active disk pushes the tracer with a force F T i = −∇ r T V d (r T − r i ) , with d = R + σ/2 the characteristic interaction distance between them, and V d is also a WCA potential as introduced in Eq. (7). Therefore, a force of equal magnitude and opposite sign is included into Eq. (5). We are interested in measuring the effective mobility µ of the tracer along the x-direction, which is primarily influenced by its interactions with the surrounding active disks. This allowed us to neglect thermal motion of the tracer in Eq. (8).
We use σ as the unit of length, the characteristic reorientation time of an active disk τ R = D −1 R as the unit of time, and k B T as the unit of energy. Assuming that D = D R σ 2 /3 , the mobility of an active disk takes the value µ A = 1/3 σ 2 /(τ R k B T) . Disks are placed in a box of area A and are subjected to periodic boundary conditions. The active bath is described by two dimensionless parameters: the Péclet number Pe = σ v A /D and the area packing fraction of disks φ = Nσ 2 π/(4A).
We set N = 10, 000 , φ = 0.12 , ε/k B T = 100 , 2R/σ = 8 and consider baths characterized by various Pe numbers, including the limiting case of passive disks, Pe = 0 . The amplitude F and frequency ω of the external force are varied according to the protocol explained in the "Results" section. Equations (5)-(8) are numerically integrated using an Euler scheme with a time step of �t = 10 −5 τ R . We first ensure that the tracer in the active bath has reached a steady state in the absence of a driving force. Then the external force is introduced, and data are collected in simulations running over a time period of 40τ R and averaged over 100 independent simulation runs, unless stated otherwise. The size of symbols is chosen such that it matches the size of error bars where applicable (Figs. 2b, 4, 6 and 7).