Active microrheology in corrugated channels

We analyze the dynamics of a tracer particle embedded in a bath of hard spheres confined in a channel of varying section. By means of Brownian dynamics simulations we apply a constant force on the tracer particle and discuss the dependence of its mobility on the relative magnitude of the external force with respect to the entropic force induced by the confinement. A simple theoretical one-dimensional model is also derived, where the contribution from particle-particle and particle-wall interactions is taken from simulations with no external force. Our results show that the mobility of the tracer is strongly affected by the confinement. The tracer velocity in the force direction has a maximum close to the neck of the channel, in agreement with the theory for small forces. Upon increasing the external force, the tracer is effectively confined to the central part of the channel and the velocity modulation decreases, what cannot be reproduced by the theory. This deviation marks the regime of validity of linear response. Surprisingly, when the channel section is not constant the effective friction coefficient is reduced as compared to the case of a plane channel. The transversal velocity, which cannot be studied with our model, follows the qualitatively the derivative of the channel section, in agreement previous theoretical calculations for the tracer diffusivity in equilibrium.


I. INTRODUCTION
Understanding transport of ions, molecules, cells and colloids in nano-and micro-fluidic devices is of primary relevance for its biological and technological applications. For example, the transport across synthetic [1][2][3] and biological [4,5] channels and pores is controlled by their shape, as well as by the effective interactions between channel walls and the transported objects. Similarly, in micro-and nano-fluidic circuitry the shape of the channel has been exploited to realize fluidic transistors [6] or diodes [7][8][9] and to control ionic [10,11] and electro-osmotic [12] fluxes.
At larger scales, the transport of colloids [13][14][15][16], polymers [17][18][19] and even active particles [20][21][22] has shown a sensitivity on the geometry of the confining channel. Interestingly, for hard sphere baths, polymer solutions and colloidal suspension the mutual interactions among peers can modulate the "bare" transport coefficient provided by the solvent. Therefore, in these scenarios the many body effect will play a relevant role in determining the effective transport performance.
Theoretical models have been proposed to describe the dynamics of a tracer particle, both in the active (forced) regime, or in the passive (unforced) case. For the active mode, the problem is typically reduced to one dimension (along the channel axis), the Fick-Jacobs equation, where the channel modulation enters as an effective potential (so-called entropic barrier) [23][24][25][26]. However, a channel with varying section induces anisotropic diffusion, which is only captured when the dynamics is studied in (at least) two dimensions, resulting in a diffusion matrix with nonzero out-of-diagonal terms. In equilibrium [14], it is found that the diffusivity in the longitudinal direction shows a maximum in the channel neck, whereas the transversal diffusion follows qualitatively the derivative of the channel section.
In this work, we study the dynamics of a forced tracer in a colloidal system confined in a corrugated channel with Langevin dynamics simulations and a theoretical model based on the Fick-Jacobs approach. In the simulations, the force has been varied covering the linear and the non-linear regimes while in the theoretical model only the small force regime can be studied. Measuring the steady tracer velocity, allows the determination of the longitudinal and transversal friction coefficients (one diagonal and one out-of-diagonal components of the friction tensor, respectively). Our results show that the effective friction experienced by the tracer particle is strongly affected by both, the geometry of the confining channel and the magnitude of the external force. Surprisingly our data show that effective friction can be reduced upon increasing the corrugation of the channel, i.e. a plane channel does not provide optimal transport. The linear regime at small forces is identified by the linear dependence of the tracer flux with the external force, and allows the application of results from equilibrium. For large forces, the tracer dynamics becomes increasingly dominated by the external force, with a small contribution from the channel corrugation.

II. SIMULATION DETAILS
FIG. 1: Snapshot of the system with Lz = Lx = 6 a and Ly = 40 a. The red particle is the tracer. The external force is parallel to the channel mid-plane (y-axis), as shown.
In the simulations, a system of quasi-hard particles is considered. All particles undergo microscopic Langevin dynamics, which for particle j reads [27]: where m j is the particle mass, and the terms in the right hand side correspond to i) the interaction between particles i and j, and with the confining wall, ii) the friction force with the solvent, proportional to the particle velocity, iii) Brownian force, and iv) the external force, which acts only on the tracer (labeled as j = 1). The Brownian force, f (t), is random in time, but its intensity is linked to the friction coefficient, γ 0 , via the Fluctuation Dissipation theorem: where k B T is the thermal energy [27]. The external force F ext = F ext e y is constant and parallel to the y-axis. The use of Langevin microscopic dynamics to approximate colloidal dynamics is convenient because it gives direct access to the instantaneous velocity v j = dr j /dt, which is also used in the equation of motion, and the algorithm is more stable than pure Brownian dynamics due to the inertial term (if time step is smaller than the momentum relaxation time, m j /γ 0 , collisions are treated as in Newtonian dynamics, while single particle diffusion is obtained at longer times).
The system is confined in the z-axis between two walls, as shown in the snapshot in Fig. 1. The shape of the walls is defined by where L z is the mean separation between the walls, A is the amplitude of the corrugation and λ is its wavelength. The system has periodic boundary conditions in the XY plane (with dimensions L x × L y ). The particle-particle and particle-wall interactions are quasi-hard: respectively, where a is the particle radius, r is the center-to-center distance and d is the minimal distance from the particle center to the wall. This is calculated expanding the cosine in Eq. 2 in power series of 2πA/λ, up to second order; our simulations are thus valid for small amplitudes and large wavelength. The approximation was validated comparing the density profile of a single particle with the theoretical prediction.
In our simulations, all particles (including the tracer) have the same mass, m, and radius, a, and the volume fraction of the system is fixed in all cases, φ = 0.20. The corrugation amplitude is A = 1 a and its wavelength is λ = 20 a, and the simulation box has dimensions L x = L z = 6 a and L y = 2λ = 40 a (see the snapshot in Fig. 1). With these parameters, the system contains N = 69 particles. Despite the low number of particles, there are no significant finite size effects as shown by the results of simulations with L x = 12 a. The length, mass, and energy units are a = 1, m = 1, and k B T = 1. The solvent friction coefficient is γ 0 = 10 √ mk B T /a. The equations of motion of the tracer and bath particles were integrated with the Heun algorithm [28], with a time step of δt = 0.0005 a m/k B T . In this algorithm, the friction force is integrated analytically in the time interval δt.
The tracer is pulled with the constant force F ext and dragged through the channel. Due to the periodic boundary conditions, it travels through the simulation box several times. The tracer position and velocity distributions in the Y Z-plane are recorded in the stationary regime for different forces and analyzed below. Note that the instantaneous velocity, used in the velocity distribution, is well defined since the inertial term is kept in the Langevin microscopic dynamics. With this setup, the non-diagonal component of the diffusion tensor D Y Z can be determined from the transversal component of the tracer velocity: For comparison purposes, a planar channel has also been simulated with the same particle density and other characteristics; the width of this channel has been set equal to the mean width of the corrugated channel, i.e. L z = 6 a. Some results of this planar channel, as well as for the bulk system with the same density, are given below.

III. MODEL
In order to capture the dynamics of the driven tracer under confinement we extend the Fick-Jacobs approximation to the generic case of interacting systems. For simplicity's sake we restrict to the case of channels that are translational invariant along the x-direction.
The overdamped dynamics of the density of a non interacting system, ρ(x, y, z, t), is described by the Smoluchowski equation:ρ (x, y, z, t) = D∇ · [∇ρ(x, y, z, t) + βρ(x, y, z, t)∇W (x, y, z)] where β −1 = k B T and W (x, y, z) accounts for the geometrical confinement and for all conservative forces acting on the particles, as will be described below in detail. When the channel section is varying smoothly, ∂ y h(y) ≪ 1, then the probability distribution can be approximated by [29,30]: with where 2L x is length of the channel in the x direction (perpendicular to the force) and h 0 its mean width. Integrating Eq.(4) in dx and dz leads to with D(y) the local tracer diffusion coefficient [14,30]. a. Non-interacting systems In the absence of interactions among the particles, W contains contributions only due to the geometric confinement, interactions with the walls (ψ) and the applied external forces (F ext ) and can be written as: b. Interacting systems In the mean field approach of the Fick-Jacobs framework, particle interactions can be accounted for by terms that are quadratic in the density field. In this case, can be understood as the sum of the potential of mean force experienced by the tracer due to the interaction with its peers (via the interaction kernel W, first term in Eq.(9), the interactions with the walls (second term in Eq.(9) and the external force (last term in Eq.(9). In the steady state , ∂ t P = 0, Eq. 7 becomes a linear first order differential equation, with solution: where D 0 is the tracer bulk diffusion coefficient. Imposing the normalization of the particle probability density and periodic boundary conditions, we can determine J and Π: In the following we are interested in the case in which just one particle (the tracer) experiences the action of the external force. Accordingly, the density distribution of the "passive" particles is barely affected by the "active" motion of the tracer. Hence, for mild values of the external force, the effective potential can be expressed conveniently as where W 0 can be approximated using the equilibrium distribution, leading to with carrying all the information about the interactions among colloids and with the wall at equilibrium.
Since F 0 depends on the mutual interactions, it is hard to provide an analytical prediction. In our analysis, F 0 (y) will be obtained from the simulations, fitting the tracer position distribution from the simulations when no external force is applied. This result for F 0 (y) is then used for finite forces to calculate both the tracer density and velocity. Using the proposed splitting of F , Π can be rewritten as: which is more useful in the following derivations.
A. Small forces In the limit of weak forces, F ext ≪ k B T /a, the flux reads where we have used the definition of Π (Eq. (13)) and the normalization condition, Eq.(11). Substituting Eq.(16) into Eq. (18) and introducing the local equilibrium density profile leads to To proceed, we expand ρ eq z (y) and D(y) around the solution for a plane channel: where ǫ is a small parameter encoding the deviation of the channel section from the plane case. In particular we assume that the modulation of the channel section does not affect the volume of the channel and so the total mass inside the channel is not affected by the modulation. This implies that Similarly we assume that the density, ̺ of the bath of hard spheres can be expanded with L/2 −L/2 dy̺ 1 (y) = 0. This implies that the diffusion coefficient of the tracer (that is determined by ̺ [14]) can be also expanded Accordingly, to leading order in ǫ we get The second factor on the rhs of Eq. (24) is the modulation in the friction coefficient induced by the geometry. Interestingly, when the diffusion coefficient is independent on the geometry (D 1 = 0), as it is for an ideal gas, Eq. (24) predicts that for small modulations the effective friction coefficient is independent on the shape of the channel. However, for larger modulations, for which higher order in ǫ should be considered, Eq. (20) predicts that the flux across a corrugated channel is always smaller than in a flat channel and hence the effective friction is increased. In contrast, for interacting systems (for which D 1 = 0), when D0 dy > 0 Eq. (24) predicts that the flux is enhanced as compared to the flat channel and hence the effective friction coefficient is decreased by the confinement.

IV. RESULTS
We start showing the tracer position distribution, and then we move to the analysis of its dynamical properties. In both cases, the simulation results are compared with the theoretical model. Given that the model is expected to fail for large forces, the regimes of small and large forces are presented separately.

A. Tracer position distribution
The tracer position distributions in the channel, for different values of the external force, are presented in Fig.  2. For small and intermediate forces, the tracer distribution mimics the bath equilibrium density. For small forces the maximum probability to find to the tracer is close to the wall (see top right panel), but there is an increasing probability in the bottleneck, and, less noticeable, in the mid-plane of the channel. Eventually, the tracer density strongly deviates from the equilibrium profile for large forces, as shown in the bottom panels of Fig.2.
In order to analyze quantitatively these distributions, the tracer density is studied in different planes, see Fig. 3. Since the cross section of the channel varies, the tracer density is studied in slabs of width ∆z = 1a. The upper panel of Fig. 3 shows the differences between the tracer density profiles when the channel narrows or widens, for a mild external force F ext a = k B T . Similar to the density of bath particles (data not shown) the tracer density is larger close to the walls, with the effect being more pronounced in the neck.
The lower panel of Fig. 3 shows the tracer position distribution in the widest section of the channel, for different applied forces. As mentioned above, for weak forces, F ext a ≪ k B T , the tracer accumulates at the walls, similarly to bath particles. However, upon increasing the force, the tracer density decreases close to the wall and increases in the channel mid-plane (around z = 0). As the velocity of the tracer increases, its translocation time across subsequent bottlenecks becomes smaller than the diffusion time along the transverse direction and the tracer cannot explore the wider parts of the channel, getting effectively confined about the mid-plane of the channel. It is worth studying also an equivalent, uniform and planar channel (Fig. 4). In this case, the tracer position distribution grows close to the wall and in the mid-plane for increasing forces (and decreases for intermediate values of z). Namely, in the planar channel, the tracer moves preferentially close to the walls for large forces, contrary to the case of the corrugated channel.
In order to compare with our model, we focus on the density profiles integrated in the XZ plane. In the simulations, a slab of width ∆y = 1 a, parallel to the XZ plane is used to calculate the average: where V is the integration volume, and P z is a constant introduced to normalize the tracer density, λ/2 −λ/2 ρ z (y)dy = 1. This is presented in Fig. 5 for different values of the external force in the corrugated channel, including the bath density distribution (F ext = 0). In the latter case, the integrated density is modulated by the channel, with the minimum in the channel neck. For increasing forces, the modulation decreases and displaces to the right, indicating that the channel is explored less efficiently in the transversal direction when the force increases. Interestingly, even the weakest applied force, F ext = 0.1 k B T /a, provokes noticeable deviations of the tracer distribution with respect to the equilibrium bath density.
From the equilibrium density profile, F 0 can be estimated and, once plugged into Eq.(15), used to predict the tracer density profile under the action of an external force. The theoretical results, obtained by assuming D 1 = 0, show a good agreement with the simulations for small forces, but deviations are observable for large y above F ext = 0.5k B T /a, due to the accumulation of errors in the numerical integrations implied in the theoretical calculation (see Fig.5). For forces above F ext = 1k B T /a, the theoretical calculation does no longer predict the tracer behavior. This underlines the fact that for weak forces the transverse probability distribution is weakly affected by external force and retains its equilibrium form. In contrast for F ext > k B T /a the transverse distribution is strongly affected by the external force and the velocity obtained by using the equilibrium transverse distribution does not match the one obtained from simulations.

B. Effective friction with the bath
We focus now on the tracer dynamics, analyzing its longitudinal and transversal motion (parallel and perpendicular to the external force, respectively). The distribution of the longitudinal component of the instantaneous tracer velocity in planes perpendicular to the external force is shown in Fig. 6. Interestingly, the velocity is constant (within the statistical noise) both in the neck and the widest section, but it varies close to the wall when the channel cross-section is changing. When the channel narrows, the longitudinal velocity is smaller close to the wall, whereas it increases when the channel widens. Additionally, it can be seen in the figure that close to the channel mid-plane, the velocity does not show any dependence on its location within the channel. In the equivalent planar channel, the velocity is almost constant (increasing slightly in the mid-plane and close to the wall).
In order to analyze the impact of the channel constriction on the tracer dynamics, we average its longitudinal velocity in slabs perpendicular to the external force, v y (y) = 1 V ρ(x, y, z)dx dz V v y (x, y, z)ρ(x, y, z)dx dz (26) Fig. 7, presents v y from simulations and theoretical predictions as the magnitude of the applied force varies. In order to compare the impact of the channel on the motion of the tracer, we normalize its velocity with the one of an isolated tracer in the bulk, v 0 = F ext /γ 0 . The ratio gives us direct information on the tracer longitudinal diffusivity since v y /v 0 = D Y /D 0 , where D Y stands for the average of the local longitudinal tracer diffusivity, D Y , over the transverse channel section. The average velocity along the channel for very small forces has a maximum in the narrowest point, but changes notably for increasing force. For large forces, the profile is almost flat, shifting to larger velocities. The comparison with the theory (thin lines) is possible only for small forces, and the agreement is semi-quantitative. It is worth noting that the contribution from interactions between particles and with the wall are encoded in F 0 , which is obtained from the tracer density profile in equilibrium. The theory captures nicely both the displacement of the maximum to larger positions, and the decrease of the maximum. For large forces, the contribution of the corrugation can be disregarded and the model predicts a constant velocity profile growing linearly with the force, which is similar to the simulation results. It must be also mentioned that the simulation results are compatible with the theoretical expectation from [14], where the full diffusion tensor of particles (without external force) in a corrugated channel are calculated. In that case, the diagonal term in the direction of the channel has a maximum in the neck of the channel.
Despite the strong variations of the velocity with the position in the channel, the flux must be homogeneous, as expected in the stationary regime. This is confirmed in Fig. 8, where the flux has been calculated as J(y) = ρ(y) v y (y) (27)   and normalized with the stationary velocity of the single particle, F ext /γ 0 , in bulk. Within the statistical noise, for small forces the flux is proportional to the force, corresponding to the linear regime shown in Eq. 18, and increases for large external forces.
The averaged flux provides a robust means to extract the effective friction coefficient experienced by the tracer particle: Fig. 9 compares the dependence of γ eff on the external force for a corrugated channel, an equivalent planar channel and in bulk, the latter for a system with N = 1000 particles. Overall, the three sets of data follow the same generic trend as that reported for bulk systems [31], with a low-force plateau, force thinning, and (apparently) a high-force plateau. In the limit of small forces, F ext < 1 − 2k B T /a, the friction is almost constant for both channels and for the bulk case, identifying the linear response regime. In this regime the presence of the confining walls induces a larger friction in both channels as compared to the bulk. For large forces, on the other hand, the plateau at high forces is above 1 in the three sets of simulation data, while the theory sets it at γ/γ 0 = 1. The origin of this discrepancy was first shown by Squires and Brady in the bulk, within the theoretical framework of the two-particle Smoluchowski equation [32]. Comparing the three cases, the planar channel always shows a higher friction coefficient than the corrugated one [33]. Interestingly, a crossover is observed for larger forces, when the friction of the corrugated channel is smaller than that in the bulk system. As shown above, this stems from the effective confinement of the tracer in the central region of the channel, decreasing the friction experienced by the tracer in the corrugated channel.
Comparing the reduction in the effective friction coefficient shown in Fig. 9 to Eq. (24) (and assuming linear response) pinpoints that the contribution of D 1 is crucial in determining the flow, even in the regime of very small forces. Indeed, Eq. (24) shows that in order to reduce the friction coefficient by means of corrugating the channel, D 1 = 0 is a necessary condition. However, in the same regime, Fig. 7 shows that the theoretical model provides quantitatively reliable predictions even assuming D 1 = 0. Hence our data show that the different observables can display quite a different sensitivity to D 1 .
Finally, we analyze the transversal component of the velocity, giving the non-diagonal component of the diffusion tensor, D Y Z . In the planar channel this component vanishes identically (not shown), while this is not the case for the corrugated channel, as shown in Fig. 10.
The depicted transversal velocity, normalized by the corresponding tracer velocity in the bulk, deviates from zero more significantly close to the walls when the channel widens or narrows (following the wall direction). Close to the mid-plane, and both to the neck and the maximal aperture positions, this velocity component vanishes. This result cannot be discussed within the simple theory model used above, and we must turn to the full model presented previously [14]. The theoretical results (for the unforced tracer) indeed predict this behavior close to the wall. Fig. 11 presents the transversal velocity averaged in slabs perpendicular to the external force, as studied above for the density and longitudinal component of the velocity. (To avoid a vanishing average due to the channel symmetry, v z has been defined as v z = v ·n, wheren is the unit vector in the vertical direction pointing from the mid-plane to every point -upwards in the upper channel half, and downwards in the lower one-). Again, to compare different forces, the velocity is normalized with the longitudinal velocity of an isolated tracer (i.e. the ratio of the transversal diffusivity to the diffusion constant of the single particle). For small forces (upper panel), in the linear regime, the results collapse onto a master curve, which follows, qualitatively, the derivative of the wall. For large forces, on the other hand, the behaviour changes; the effect disappears with increasing forces, indicating that the tracer motion is increasingly confined in the transversal direction and confirming that the tracer does not explore the full section of the channel.

V. CONCLUSIONS
The dynamics of a tracer pulled through a colloidal system confined in a corrugated channel has been analyzed. The tracer is pulled with a constant force, and the whole range of forces, has been studied. The results are compared with a simple one dimensional model based on the Fick-Jacobs approximation, but the results from the full model, studied previously, have also been considered; the latter predicts that for the unforced tracer the diffusion tensor has non-zero out-of-diagonal terms.
Our simulations confirm these predictions in the limit of small forces, and show that the linear response regime extends up to F ext ∼ 1 − 2 k B T /a. The tracer longitudinal velocity has a maximum in the neck of the channel, whereas the transversal component is non-zero and has a maximum where the channel cross-section varies more strongly. Likewise, in this region both the longitudinal and transversal velocities (or local diffusion constants) vary close to the walls, while they remain essentially constant in the rest of the channel.
The theoretical model describes the tracer dynamics effectively, fitting the contribution from particle-particle and particle-wall collisions in the equilibrium case (F ext = 0), and using this result for finite forces. The results for the tracer density and longitudinal diffusivity agree almost quantitatively with the simulation results within the linear regime. Outside this, the calculations break down and cannot provide reliable results.
For larger forces, the tracer is confined to a narrow region parallel to the channel axis, set by the minimum crosssection of the channel. The longitudinal component of the velocity in this region is almost constant, as the channel section is not explored, and the transversal component becomes increasingly small. As a result, the effective friction experienced by a tracer pulled with a large force in the corrugated channel is smaller than in the bulk with the same density. This region falls out of the theoretical model developed here.