Elasticity-based polymer sorting in active fluids: A Brownian dynamics study

While the dynamics of polymer chains in equilibrium media is well understood by now, the polymer dynamics in active non-equilibrium environments can be very different. Here we study the dynamics of polymers in a viscous medium containing self-propelled particles in two dimensions by using Brownian dynamics simulations. We find that the polymer center of mass exhibits a superdiffusive motion at short to intermediate times and the motion turns normal at long times, but with a greatly enhanced diffusivity. Interestingly, the long time diffusivity shows a non-monotonic behavior as a function of the chain length and stiffness. We analyze how the polymer conformation and the accumulation of the self-propelled particles, and therefore the directed motion of the polymer, are correlated. At the point of maximal polymer diffusivity, the polymer has preferentially bent conformations maintained by the balance between the chain elasticity and the propelling force generated by the active particles. We also consider the barrier crossing dynamics of actively-driven polymers in a double-well potential. The barrier crossing times are demonstrated to have a peculiar non-monotonic dependence, related to that of the diffusivity. This effect can be potentially utilized for sorting of polymers from solutions in \textit{in vitro} experiments.

Here we study the dynamics of polymer chains in two dimensions (2D) in the presence of active Brownian particles (ABP) [51] by using Brownian dynamics simulations. The recent studies on swelling, collapse, and looping of actively-driven polymers [36][37][38][39][40] served as a starting point for the current investigation, with "active polymers" being a perspective research direction [6]. Such a 2D system is more relevant to the in vitro experimental setups, rather than to in vivo settings, as the former are frequently carried out in quasi-2D setups [25].
The polymer chain in the bath of ABPs is not in equilibrium and unusual behaviors can take place [36][37][38][39][40][41][42][43][44][45][46][47][48][49][50]. In this study we find that due to propelling forces of ABPs, the polymer dynamics is greatly facilitated. In particular, the polymer center of mass (COM) diffusivity shows a non-monotonic behavior as a function of both the chain length L and its bending stiffness κ. The polymer at maximal diffusivity has preferentially bent conformations, maintained by the balance of chain elasticity and propelling forces of ABPs. We also consider the barrier crossing dynamics of polymers in a double-well potential, finding that the crossing times are non-monotonous with L and κ, too. These results can potentially be utilized for separating polymers based on their length or stiffness [53] via using active fluids.
The paper is organized as follows. We introduce the model and simulation methods in Sec. II. The main results on the polymer diffusive behavior are presented in Sec. III and the barrier crossing dynamics of polymer chains is investigated in Sec. IV. Finally, we summarize and discuss our results in Sec. V. (Bottom) Configuration of a polymer (blue chain) in a bath of ABPs (red spheres). Here n=32, κ=360, and the packing fraction of ABPs is φ = 0.05. The figure was rendered using VMD [52]. The video files illustrating the chain dynamics are provided in the Supplementary Material.

II. MODEL AND METHODS
We perform Brownian dynamics simulations of semiflexible polymers in the presence of active particles [2] on the plane in 2D. As a representation of active particles, we use a model of self-propelled ABPs, which can be man-made and used in experiments [6,51]. The position of the ith ABP at time t is described by the overdamped Langevin equation [6,51] where µ is the particle mobility, ∇ ≡x ∂ ∂x +ŷ ∂ ∂y , and ξ i (t) is the two-dimensional Gaussian white noise with unit variance in each dimension, ξ i (t) · ξ i (t ) = 2δ i,i δ(t − t ) and δ i,i is the Kronecker symbol. Potential U denotes the interaction potential between different ABPs, ABPs and polymer beads, and ABPs with external potential. We refer the reader to Sec. VII of Ref. [6] for the discussion of underdamped dynamics of active particles. Also, the recent examination of inertia effects in some anomalous diffusion processes is instructive [54].
The particle moves with a constant speed v a along the direction given by angle ψ This angle is subjected to rotational diffusion, as described by the rotational Langevin equation, where D r is the rotational diffusivity and ξ r is the Gaussian white noise with unit variance. The rotational diffusion leads to the decorrelation of particle velocity on the time scale of τ = 2/D r . For the case of spherical particles of diameter σ the value of D r is related to its translational diffusivity D t as [55] The strength of particle propulsion is measured in terms of the Péclet number The situation of v a = 0 corresponds to passive Brownian particle, studied previously in the context of macromolecular crowding in e.g. Refs. [56][57][58].
The polymer is modeled as the bead-spring chain of n monomers of diameter σ connected by harmonic springs with the corresponding potential where k is the spring constant and l 0 is the equilibrium bond length. Hereafter, the chain monomers are of the same size as the ABPs, see the Discussion section for the effects of ABP size. We choose the Hook's modulus as k = 10 3 k B T /σ 2 and l 0 = σ to prevent the crossing of ABPs by spring sections. The bending energy of the chain is given by where κ is the bending stiffness and θ i is the bending angle of the ith chain segment, see Fig. 1 (Top). For a given value of κ, the chain persistence length in two dimensions is then l p 2κl 0 /(k B T ). Note that the polymer behaves as a much softer chain in the presence of ABPs due to the enhanced fluctuations [39].
The effects of self-avoidance between different chain monomers and between ABPs and chain monomers are modeled by the Weeks-Chandler-Andersen (WCA) potential [59], for r i,j ≤ r cut and U WCA (r i,j ) = 0 for r i,j > r cut , where the cutoff distance is r cut = 2 1/6 σ. Here, r i,j = |r i − r j | is the inter-monomer distance and is the interaction strength. This potential corresponds to a polymer chain in a good solvent. To study the barrier crossing dynamics of the chain, we add an external double-well potential along the x-axis, Here the potential width, the distance from one potential minimum to the barrier, is ∆x = B/A and the barrier height is ∆U = B 2 /(4A). Then, the dynamics of the jth chain monomer is governed by the following overdamped Langevin equation An important parameter of the medium is the density of ABPs. To fix the density, we use the periodic boundary conditions with the square box of area L 2 . We choose L=60 that is larger than the typical size of the polymer, to prevent any artifacts of boundary conditions. The packing fraction of ABPs is defined as where N A is the number of particles and A A = π(σ/2) 2 is the surface area per ABP, see also Ref. [60]. We use hereafter φ = 0.05. At low packing fractions φ, the interaction between ABPs can be negligible. At high φ, in contrast, some clustering of ABPs and phase separation phenomena can take place [61][62][63][64][65], that is however beyond the scope of this paper.
To numerically integrate the equations of motion (2) and (10), we implement the stochastic Runge-Kutta algorithm [66]. We measure the length, the time, and the energy in units of σ, t 0 = σ 2 /D t , and thermal energy k B T , respectively. We set below the model parameters as σ = l 0 = 1, k = 10 3 , D r = 1, D t = 1/3, and = 1. The important length scales of the system are the chain length L, the persistence length l p , and the persistence length of the ABPs motion 2v a /D r . The main features of our results (shown below) will remain the same if we fix the ratio of those lengths. We use the integration time step ∆t = 2 × 10 −4 , so that in our plots, the simulation time of t = 1 corresponds to 5000 iteration steps of the evaluation scheme. Initially, the system is equilibrated for ∼ 10 6 steps and typically run up to ∼ 10 9 iteration steps.

A. From superdiffusive to normal Brownian motion
We first consider the diffusive motion of a polymer chain in the presence of ABPs and no external potential, U ex = 0. From a long trajectory of the chain generated in simulations, we calculate the time-averaged MSD (tMSD) of the polymer COM [23], is the x-coordinate of the polymer's COM. Here, ∆ is the so-called lag time along the trajectory [23]. Moreover, the tMSD is averaged over an ensemble of N independent traces recorded in simulations, with N =40 for most of the findings presented below. The tMSD along the y-axis is naturally the same as along the x-axis in the absence of an external potential. Note that the ensemble of independent trajectories which is required for a satisfactorily smooth behavior of the tMSD is substantially smaller than that for the ensemble averaged MSD [23]. tMSD is therefore frequently used in single-particle tracking experiments, where often not so many but rather long traces are generated/available, see e.g. Refs. [23][24][25]. We find that for our system the tMSD is the same as the ensemble averaged MSD (not shown), which means the system is ergodic [23]. We thus use t for the lag time ∆ below, for simplicity of the notation. We recently demonstrated that the COM motion of a two-dimensional polymer chain is greatly enhanced as the activity of ABPs increases [39]. In this study, we focus on how the polymer COM motion depends on the chain length L = nσ and its bending stiffness κ. Active driving by ABPs renders the diffusion of the polymer COM superdiffusive/non-Brownian on intermediate time scales. In Figs. 2A and 2B, we show the tMSD of the polymer COM for different chain lengths and bending stiffnesses, respectively. We compute the time-local scaling exponent of tMSD as [23] We observe in Figs. 2C and 2D that the scaling exponent drops from the values α ≈1.6-1.8 at relatively short times, to less superdiffusive α values at intermediate times, to finally the normal diffusion behavior with α = 1 at very long times. Note that here, contrary to our recent study [39], we consider overdamped dynamics for both ABPs and polymer chains and thus one does not expect to see a ballistic regime of the tMSD even in the limit of short times. The noise in α(t) data in the long time limit is due to the worsening statistics as the lag time ∆ becomes comparable to the total time T in Eq. (11), see also Ref. [23].
To better characterize the diffusive behavior of the polymer COM, in particular the origin of the superdiffusive behavior at short times, we consider the velocity auto-correlation function (VACF) of the polymer COM, where V COM (t) is the COM velocity at time t. In Fig. 3, we show the VACF(t) for the case of n=32 and κ=360 at which the super-diffusive behavior is most pronounced. The decay of VACF(t) at short time is close to a powerlaw, ∼ t −β , with the exponent β 0.6. The powerlaw decay with the β < 1 indicates that the tMSD(t) of the polymer COM, which can be calculated by double integral of the VACF(t), increases super-linearly, in the time interval, consistent with our simulation results, see Fig. 2C and 2D. For the chains with a less pronounced super-diffusive behavior, the power-law decay regime of the VACF(t) becomes shorter or disappears  [3]. This question is, however, out of our scope of this study.
At long times, on the other hand, VACF(t) decays exponentially with the correlation time τ COM , which are shown in Fig. 4 (triangles, the right axis). Physically τ COM is the time at which the persistent chain motion start to decorrelate. The correlation time shows a nonmonotonic behavior as a function of the chain length n and κ. The value of τ COM is typically, except for the n = 8 case, much longer than the persistence time of the ABPs, τ = 2/D r = 2. In comparison, the VACF of the polymer COM in the absence of ABPs is delta-correlated as is the thermal noise. For the case of ring polymers in two dimensions filled by ABPs [41], the VACF is determined by that of the ABPs, independent of the chain elastic properties. The exponential decay of VACF at long times indicates that the tMSD(t) will increase linearly with time t in this domain, consistent with our simulation results, see Fig.2.
We extract the diffusivity of the polymer COM by fitting the long time limit, in the range t = [10 3 , 10 4 ], of the tMSD(t) with a linear function. We find that the diffusivity shown in Fig. 4 (circles, the left axis) also varies non-monotonically with the chain length n and stiffness κ. The optimal values of n* and κ*, which give rise to the maximum of the diffusivity, are coincident with those of the correlation times. This indicates that the non-monotonic behavior of the diffusivity is resulting from the non-monotonicity of τ COM . For the case of a fixed chain length (Fig. 4B), the diffusivity is proportional to the correlation time. This is related to the fact that the diffusivity of ABPs is proportional to the correlation time. However, for the case of varying chain length (Fig. 4A), the relation is more complicated. The reason is that in this case, not only the correlation time, but also the number of particles pushing the polymer chain, and hence the velocity of the polymer, is changing with the chain length. The non-monotonic dependencies we observe in Figs. 4 are fairly robust with respect to varying φ and v a , but the optimal chain length and bending stiffness depend on these model parameters.
As can be seen from visualizing the results of simulations at the optimal chain length n* or stiffness κ*, at these parameters the polymer captures the surrounding ABPs for longer time. Therefore, the superdiffusive interval of the polymer COM motion becomes more prolonged, see the α(t) dependencies in Figs. 2C and 2D. To quantify these observations, we consider various quantities such as end-to-end distance distribution, gyration radii of the polymer, and the Fourier spectrum of chain conformations. As we have shown in Ref. [39], the endto-end distance distribution and the radii of gyration of polymer significantly change in the presence of ABPs. However, we do not find any distinctive features that indicate the condition of maximum diffusivity of the polymer in these two quantities. Therefore, we present only the Fourier mode analysis of chain conformations in the next subsection.

B. Chain conformations via the Fourier modes
We analyze the bending modes of the chain in terms of the Fourier amplitudes of its bending harmonics. The polymer conformations are described in terms of the tangential angle θ(s), where s = [0, L] is the arc length. Following the method proposed in Ref. [67], the chain conformations are decomposed into the cosine modes, At equilibrium, each mode evolves independently and its variance is given by [67] In what follows we use the discrete approximation of this formula, see Ref. [67] for more details. We first show the variance of the Fourier modes a m in the absence of ABPs, see empty symbols in Fig. 5A. The results match well with the theoretical prediction of Eq. (15), shown as the solid lines, for not too stiff chains (l p ≤ L). As the chain becomes stiffer, the simulation results overestimate the theoretical values of a 2 m . The deviation is due to the additional "stretching fluctuations" of the chain in our model so that, as the spring constant k increases, the simulations results (not shown) become closer to the theoretical prediction. However, for the reasons of computational efficiency, we do not use larger spring constants here. In Fig. 5B we show variance of Fourier modes (filled symbols) in the presence of the ABPs. The variances are 1-2 orders of magnitude larger, depending on the mode m, compared to that of the chain at equilibrium in the absence of active particles. The enhancement, defined as the ratio of variance in the presence of ABPs to that in the absence of ABPs, is more significant for smaller mode number m as shown in the inset of Fig. 5B. We also find that the enhancement of the first two modes (m=1 and 2), which determine the large length scale of the polymer conformations, is the largest for the chain of optimal stiffness (κ=360). This finding indicates that chain conformational fluctuations are highly correlated with the enhanced diffusivity.
The enhancement of fluctuations in actively-driven systems was measured experimentally, among others, for microtubules in the presence of myosin motors [42]. This effect can be interpreted as an increase of the effective temperature in the system, known to facilitate the polymer looping kinetics [39]. In contrast to equilibrium systems, here the fluctuations are mainly due to collisions between ABPs and polymer chain, but the energy dissipation occurs via all length scales of the polymer. Thus, it is not surprising that the variance of a m does not follow the equilibrium scaling relation of Eq. (15). Note also that the distribution of a m amplitudes is always Gaussian in equilibrium systems [68]. Physically, this Gaussianity is due to the absence of correlations in the thermal noise.
In the presence of ABPs, however, the distributions of the Fourier amplitudes become strongly non-Gaussian, as we exemplify in Fig. 6 for the 1st and 2nd mode. We find that for small stiffness parameters κ the distribution p(a 1 ) is roughly uniform in a broad range and for very large κ the distributions p(a 1 ) and p(a 2 ) reveal a single peak. For intermediate chain stiffness values κ, the distribution p(a 1 ) becomes bimodal, which means the chain adopts preferentially bent conformations. The first Fourier mode corresponds to the half-period of the cosine function and each peak in p(a 1 ) correspond to the chains in the bent shapes of "⊂" and "⊃". Such polymer shapes are maintained by the balance between the elastic chain energy and the propelling force of ABPs (compare also to the conformations of actively driven fluid membranes [55]). This finding is consistent with the mechanism pro-posed in Ref. [40] and our analysis provides a quantitative evidence of the mechanism. For the higher modes (with m ≥ 2), we find that the distributions are unimodal, but also exhibit non-Gaussian features (the m = 2 case is shown in Fig.6B) To summarize, in this section we find that the polymer chains in the presence of ABPs reveal a non-monotonous diffusivity as functions of the chain length n and bending stiffness κ. At the optimal chain length or stiffness, the polymer has preferentially bent conformations maintained by its elasticity and propelling forces of ABPs. In the next section we show how this effect can be utilized for polymer sorting.

IV. BARRIER CROSSING OF POLYMERS IN ACTIVE FLUIDS
Here, we examine the barrier crossing problem for the actively-driven polymers in a double-well potential. The barrier crossing dynamics of polymers in equilibrium media were considered in a number of recent studies [69][70][71][72][73]. The crossing times were shown to be strongly dependent on the properties of conformational rearrangement of the polymer in an external potential. In the absence of ABPs, the crossing times of the polymer chains can be rather long, even if the potential barriers are rather low, ∆U ≤ 1k B T , because all chain monomers need to cross the barrier at the same time. In the presence of ABPs with large v a values, however, the chain can cross the barrier much shorter time due to the enhanced fluctuations, see Fig. 5.
Here we consider both the polymers and the ABPs to be subjected to an external potential U ex (x) acting along the x-axis, as expressed in Eq. (9). Depending on the ratio of the potential width and polymer gyration radius, different barrier crossing scenarios can realize [69][70][71][72][73]. We have chosen the potential width ∆X = 30σ is larger than the typical polymer size, so that the polymer COM is placed in either left or right well of the potential. The barrier height also needs to be carefully chosen; if it is too high, the crossing time can be enormously long, but if it is too low, the polymer will move freely in the potential. Here we use ∆U = 6k B T to make the barrier substantial even in the presence of ABPs, but not too high so that a sufficient number of crossing events occurs during our simulation time. In comparison, with the same barrier height but with less active ABPs (smaller v a ) or in the absence of ABPs, the barrier crossing events happened extremely rarely (not shown).
As mentioned before, the distribution of ABPs in the external potential U ex (x) can deviate strongly from the Boltzmann distribution, and the general form of this distribution is not known for a given system. Here, the external force acting on the ABPs is much weaker than the active force, namely so that the distribution of ABPs is barely affected by the potential. On the other hand, the polymer chains are confined in one of the potential minima. We track the COM coordinate X(t) of the polymer which shows a hopping dynamics between the two minima of the potential. We define the crossing time, T cr , as the mean first passage time of the polymer COM from one potential minimum to the other. Figure 7 shows the dependence of T cr on the chain length nσ and bending stiffness κ. The crossing time shows a minimum both at a certain polymerization degree and bending stiffness of the polymer.
Previously, it was shown that the polymer barrier crossing time in equilibrium can be a non-monotonic function of the polymer length or bending stiffness [69,[71][72][73]. In those studies, the polymer barrier crossing dynamics was mapped onto a one-dimensional barrier crossing processes by considering only the COM coordinate and the remaining degrees of the freedom were taken into account by the effective free energy (also known as a potential of the mean force) of the polymer COM. The barrier height of the free energy can be a non-monotonic function of the chain length or stiffness [69,[71][72][73], due to the chain conformational changes, which is the reason of non-monotonic behaviors of the crossing time. Following the approach of Ref. [73], we obtain numerically the effective free energy by using the Boltzmann inversion of the distribution function P COM (X) along the x-coordinate, namely The effective free energy also exhibits a double-well potential (not shown) and the barrier height of the potential dU is shown in the insets of Fig.7. The effective barrier height is much smaller than the real barrier height ∆U and shows a non-monotonic behavior as a function of n or κ. Following the Kramers' barrier crossing theory [75], we assume that the crossing time of the polymer COM scales as with the effective diffusivity D of the polymer COM and the effective barrier height dU . The non-monotonous behavior of both dU and D gives rise in a dramatic nonmonotonous behavior of the crossing time. In comparison, in equilibrium only the effective barrier height dU can show a nonmonotonous behavior, see Refs. [69,[71][72][73].

V. DISCUSSION AND SUMMARY
Active fluids are inherently out of equilibrium those are of relevance for a number of living systems. The physical understanding of those systems is, however, far from being complete [4,6]. Even some basic properties, for example, the distribution of the active particles in external potentials U ex is only known in some simple setups [27,28].
Here we numerically studied the dynamics of the actively-driven semiflexible polymers in two dimensions. We found that the ABPs are accumulated in the concave regions of the chain, which results in superdiffusive motion of its COM at intermediate times. At long times, the diffusive motion of the polymer COM becomes normal, but with the diffusivity which was much higher than for the motion without active driving. The diffusivity revealed a maximum versus the polymer length and bending stiffness. The chain at the optimal length or optimal stiffness had preferentially bent conformations, as we have demonstrated examining the chain conformations in the Fourier modes. This occurs when the polymer elastic force and the propelling forces of the ABPs are balanced.
As an application of the nontrivial behavior of the polymer diffusion in active fluids, we also considered the polymer dynamics with the ABPs in the presence of a double-well potential. We found that, as the activity of ABPs increases, the crossing time is shown to be greatly decreased. The crossing time also showed a non-monotonous dependence with the chain length n and bending stiffness κ. This is because of a non-monotonous behavior of the diffusivity and the effective barrier height versus n and κ. This suggests that the polymer chains can be separated from mixtures based on their length or bending stiffness [40], those are important for a number of practical applications [53]. This scenario might be possible experimentally, for example for sorting of stiff biopolymers such as microtubules or actin filaments immersed in a fluid of active colloidal particles or bacteria [6] in combination with microfabrication channel [11]. In experiments, it will be important to choose proper system parameters; the activity of the fluids and the potential barrier should be comparable. The former can be controlled, for example by varying the energy source in the fluids [3,25], and the latter, by designing the geometry of the channel [11].
Note that in our study we considered a single chain in a simulation box, hence the polymer-polymer interactions were absent. For a mixture of many chains, the polymerpolymer interactions could change the crossing dynamics. However, if the polymer density is not very high, the effects of interactions should be minor, not affecting our main findings and trends. Our simple setup with the close-contact potentials neglects also the long-ranged hydrodynamics interactions [1,6,74]. The latter can govern, among others, some energy transfer reactions and tune collective effects in actively-driven systems, such as those in a 2D diffusion of micron-sized spheres driven by swimming bacteria [25]. Nevertheless, we expect the main features of our findings to stay valid in real systems, particularly when the hydrodynamic effects can be accounted for via a renormalized friction coefficient.
Finally, we have considered that the ABPs are of the same size as the chain monomers. For the case of bigger ABPs, we found that the "capturing" of ABPs by the polymer is not possible and the non-monotonous behavior of the diffusivity of the polymer COM disappears (results not shown). For smaller ABPs, since the rotational diffusivity scales with the diameter of the particle as D r ∼ σ −3 , the persistence length of the ABPs' motion decreases very rapidly as ∼ σ 3 . In this case, the distribution of ABPs can be mapped onto the Boltzmann-like distribution, but with a higher "effective" temperature [27], and the non-monotonous diffusive behavior of the polymer chain will disappear. The typical size of selfpropelling colloidal particles is in the range of 0.1-10 µm and the length of the biopolymers such as microtubules can be up to 10 µm long. Therefore, it would be possible to choose the proper experimental parameters that could demonstrate the validity of our main findings experimentally.