Flow pumping by external periodic shear applied to a soft interface

Flow pumping in viscous fluids is of prime importance in micro-fluidic applications. Here we show that a single colloidal particle in front of a soft wall, manipulated by external means like an optical tweezer, can pump the ambient viscous fluid. The particle, moving back and forth parallel to the soft wall, can produce an averaged net flow in a direction perpendicular to the wall. Using a perturbative scheme, we present the results. Analysis show that this flow in terms of capillary number, scales as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {Ca}}^2$$\end{document}Ca2.

www.nature.com/scientificreports/ numerical study on the dynamics of a time-reversible squirmer illustrated that it is able to move near a weakly deformable interface 21 .
The rest of this article is organized in the following way. In "Model" section, we first introduce the model that we want to study. Then in "Governing equations" and "Perturbation analysis"sections, we present the governing equations and the details of perturbation scheme that we use to solve the equations. In "Steady state results"section, we present the flow pattern and hydrodynamic force distribution for an auxiliary steady state problem that will be used to obtain the solutions to our real problem. Numerical and scaling results for the pumping flow will be presented in "Average pumping flow " and "Dimensional analysis" sections, respectively.

Model
As shown in Fig. 1, consider a spherical micro-particle with radius a that is moving near a soft and initially flat wall. This flexible wall separates two immiscible fluid phases which are characterized by their densities and viscosities given by (η, ρ) and (η ′ , ρ ′ ) , respectively. The sphere is located in the phase that is denoted by unprimed variables. The soft wall can either be a simple interface or a physical membrane with respective elastic modulus. Here, for simplicity, we consider a soft interface with an isotropic surface tension denoted by γ . With the setup that is introduced here, we aim to consider the hydrodynamic flow pattern near an interface that is under a tangential periodic shear force. The micro-particle that is trapped by an optical tweezer, is the source of this shear force. As an approximation that can further simplify the problem, it is assumed that the rotational motion of the particle is balanced by a suitable external torque. Along this task, we let the trapped particle to move in a prescribed 1-dimensional periodic trajectory parallel to the interface.
To parametrize the geometry, as shown in Fig. 1, we choose a laboratory reference frame where the position of the interface, in its undeformed state, is given by z = 0 plane. As shown in the picture, both cylindrical variables (r, φ, z) and Cartesian variables (x, y, z) will be used to parametrize the problem. In this coordinate system and in terms of Cartesian variables, the position of the micro-particle is given by: where h denotes the distance between sphere and interface, x 0 stands for the amplitude of motion and ω shows the frequency of the motion. The instantaneous velocity of the sphere is given by U p (t) = U p cos(ωt)x with U p = ωx p . As a result of the shear forces exerted from moving sphere, the interface will be deformed. In the laboratory frame, time dependent deformation pattern of the interface can be denote by z = ζ(x, y, t).
In addition to the flow pattern that is generated in the fluid, the deformation profile z = ζ(x, y, t) is among the unknown variables that we aim to determine. In the following part we will present the necessary dynamical equations that will help us to calculate the unknown variables.

Governing equations
To investigate the dynamics of this system, we denote by u(x) and p(x) , velocity and pressure fields in the fluid medium. To distinguish the field variables at both sides of the interface, we will use the primed and unprimed symbols for top and bottom fluids, respectively.
To write the dynamical equations, we benefit the approximations that are essential for our micron-scale set up. Denoting the dimensionless Reynolds number by Re = ρU p a/η , we will assume that Re ≪ 1 . This will restrict the dynamics of the fluid to the viscosity dominated regime. In this limit, the Stokes and continuity equations will determine the flow filed. In terms of the dynamical variables, these equations read as: Similar equations will be held for primed variables. It should be noted that in general, the low Reynolds condition does not guarantee the steadiness of the flow. In our system for Re ≪ 1 , the flow pattern is assumed to reach an instantaneous steady state. As a result of this approximation, time does not appear in the Stokes equation explicitly. However as the sphere moves slowly in the medium, the flow will have a quasi-steady time evolution. (1) η∇ 2 u − ∇p = 0, ∇ · u = 0.

Figure 1.
A soft interface separates two immiscible fluids and the shear exerted on this interface from a moving micro-particle, will deform it. The deformation pattern shown here, is pictorial and it does not necessarily shows a real pattern for a specific motion of a micro-particle. www.nature.com/scientificreports/ The quasi-steady approximation will allow us to consider a true steady state problem from which we can extract the results for our quasi-steady state problem. Let us define the steady state (instantaneous) problem with a micro-particle that is located at x s p = (0, 0, −h) with its surface velocity given by U p = U px . We denote the deformation pattern and velocity field of this problem by ζ s (x, y) and u s (x, y, z) , respectively. The solution to our real quasi-steady state problem at an arbitrary time t and in laboratory frame, can be written as: where, in the steady state solutions we have replaced the position and velocity of the particle by: U p = U p (t) and Boundary conditions. Boundary conditions for the above mentioned steady state problem, are given by: where the last equation represents the continuity of velocity field at the position of interface.
Regarding the energy stored in the interface, both the tension and gravitational energy, the components of the stress tensor, obey the following boundary conditions at the position of interface 16 : where n and t denote two unit vectors that are locally normal and tangent to the interface and (1/2)∇ ·n denotes the mean curvature of the interface. Density contrast is denoted by �ρ = ρ − ρ ′ . We note that in terms of the interface shape function H = z − ζ s , the normal vector can be written as: n = ∇H/|∇H| . Discontinuity in the normal stress that is reflected in the last part of the above equation, is the main mechanism that can initiate deformation in the membrane. Regarding the detail structure of this equation, three different forces namely, viscous, tension and gravity, compete in force balance. Relative strength of these forces can be given in terms of two dimensionless numbers: where, U p is a velocity scale that was defined before. Including the ratio between viscosities, = η ′ /η , these three dimensionless numbers can be used to determine the physical state of the system under investigation. For very small viscous tension, α −1 ≪ 1 and β −1 ≪ 1 , hence the membrane remains flat. We will use this observation in developing a perturbative procedure to solve the dynamical equations. It is worth mentioning that a new length scale, namely the capillary length can be defined as: ℓ c = γ /g�ρ . At scales much larger than this capillary length, the gravity dominates over the surface tension. Capillary length is a relevant length scale when the gravity and surface tension have the same order of magnitudes. In other regime where the effects of gravity are negligible, a dimensionless number so called the capillary number Ca = β −1 plays the role of a control parameter instead of capillary length. We present our main results in the latter regime.
To finish with the formulation, we need an equation that couples the dynamics of the interface and the fluid. We note that in the laboratory frame and as a result of no-slip boundary condition on the wall, the trajectory of a particular point of the interface should follow exactly the path of the fluid particles that are locally in contact with points on the wall. In terms of the fields in laboratory frame, we have n · u(x, y, z = ζ ) = ∂ t ζ / √ 1 + ∇ζ · ∇ζ , where ∂ t = ∂/∂t . In terms of the fields corresponding to the previously defined steady state problem, we will have: The set of boundary Eqs. (4), (5) and (7), will allow us to find complete information about the shape of interface and the velocity fields on both sides of the interface.

Perturbation analysis
Even though the governing equations at Re = 0 are linear in terms of velocity field, the boundary conditions at the position of interface are highly nonlinear with respect to the deformation profile, ζ s (x, y) . This nonlinearity in terms of ζ s , will result velocity profiles at both sides of the interface that are nonlinear in terms of the sphere velocity U p . Our strategy to overcome this nonlinearity, is to use some approximations. It is obvious that for either cases of a system with hard interface ( β −1 = Ca ≪ 1 ) or a system with large density contrast or large gravity ( α −1 ≪ 1 ), we expect a very small deformation pattern. Regarding this observation, we can define a small parameter denoted by ǫ = (α + β) −1 and set up a perturbative expansion of the required fields in terms of this small dimensionless quantity. Formal expansion of all fields read as: As it is expected, at zeroth order of ǫ , the interface remains flat and the equations are linear. To be mathematically precise, one should note that the above perturbative expansion should be written for dimensionless www.nature.com/scientificreports/ variables. Using characteristic length and time scales given by a and (a/U p ) and denoting the scales for pressure by ηU p /a and η ′ U p /a , we can transform the equations to dimensionless form. Following such transformation, small quantity denoted by ǫ , will arise naturally. However, to keep the notation as minimal as possible, we prefer to proceed with the dimensional form of the equations. The above perturbation analysis can work in all regimes but, as we stressed before, we will mainly concentrate on the regime where the capillary number Ca plays the main role. In this regime β ≫ α and ǫ ∼ Ca . In the following parts, we will present the equations and results up to the leading orders of ǫ.
Zero order equations. The equations for the zeroth order fields are given by: with similar equations for primed fields. The flow fields are subjected to the following boundary conditions: where l can be either x or y.
First order equations. Having in hand the zeroth order results for the flow fields, we can calculate the first order deformation by solving the following differential equation: After calculating the interface shape at the first order, we will be able to obtain the first order correction to the flow profile. At O (ǫ) , equations for the flow fields are given by: and corresponding boundary conditions are given by:

Steady state results
Having introduced the steady state equations for zeroth and first order fields, we can present their solutions. Our final goal in this part is to obtain the first order flow and the first order hydrodynamic force acting on the micro-particle.
Zeroth order results. So far, we have introduced the perturbative equations for the steady state problem.
As it is reflected in Eqs. (9) and (10), the zero order flow corresponds to the flow due to a moving sphere near a flat wall with vanishing normal velocity and continuous tangential velocity on the wall. This is a classical fluid dynamic problem with a general solution that can be obtained by using the Lorentz reciprocal theorem 22 . Using such method, one can obtain results that take into account the finite radius of the sphere. However, here and to avoid the complicated mathematical structure of the results, we restrict ourselves to a limiting case of a very small sphere and expand all the results in powers of small parameter (a/h). At the leading order of (a/h), the velocity field can be presented in terms of the flow corresponding to a point force f 0 = 6πηaU p . In terms of the Green's function, the velocity can be written as: 12 where r = (x, y, z + h) , R = (x, y, z − h) , and i, j, k can take values of 1, 2, 3 but α is restricted to 1 or 2. Using the Lorentz reciprocal theorem, the hydrodynamic force acting on the sphere, can be written as: 22 www.nature.com/scientificreports/ In terms of (h/a), this drag is plotted in Fig. 2a. As one can see, for = η ′ /η ≥ 2/3 , increasing the distance from wall will decrease the drag until it reaches to its asymptotic value, corresponding to the drag in a homogenous single-phase medium.
Note that the case with = 1 does not correspond to a single fluid system. To see the reason, we note that for the zeroth order problem, distribution of normal stress over the wall does not vanish in the limit of = 1 . This normal stress at the zeroth order, is the driving force that will deform the wall.
Normal stress distribution over the interface is another important quantity that we can calculate. The leading order normal stress difference at small (a/h), can be written as 22 : where R(r) = Ŵ 0 9h 2 ra 2 (h 2 +r 2 ) 5/2 . In the following part we will use this result to obtain the deformation pattern of the membrane.
First order deformation. After calculating the lowest order velocity profile and subsequently obtaining the lowest order stress tensor, we will be able to study the first order correction to the deformation pattern. According to Eq. (11), the discontinuity of the stress tensor (its normal component) can be used to calculate the deformation pattern. Using Eq. (11), the deformation pattern in Fourier space can be written as: Integrating over wave vector, we will arrive at the following result for the deformation field: For a case with very large capillary length ( ℓ c ≫ a or ǫ ∼ Ca ), the shape of the interface can be achieved as: Figure 3 shows an example of the interfacial deformation pattern for specific values of the parameters. As it is seen in this figure, the interface experiences a sink (rise) in front (back) of the moving sphere. The front sink observed in the interface, is due to the positive pressure produced by a moving sphere in its front.
First order flow. As we showed before, the zeroth order results do not take into account the deformation of the interface. Our main interest in this article is the first order correction to the velocity profile that takes into account the flexibility of the membrane. As we will see later, this contribution will give a net flow for a periodic and reciprocal motion of the sphere when we average it over time.
Eqs. (12) and (13) show that the first order problem is a stokes flow with vanishing body force that is subjected to an inlet normal velocity profile given at the position of undeformed interface ( z = 0 ). Moreover, for this flow, the tangential velocity and the components of the stress tensor are discontinuous at z = 0 and the velocity must vanish on the surface of the sphere. Having in hand the zeroth order flow and first order deformation pattern, the boundary normal velocity can be obtained. As a function of ρ and φ in the cylindrical coordinate system, this boundary velocity has the following form: (a) (b) Figure 2. (a) The zeroth order hydrodynamic force acting on the particle in terms of the distance to wall is plotted. This force is parallel to the wall. (b) The first order force is perpendicular to the wall and it is plotted in terms of the distance to wall. Parameters are U p = U px with U p > 0 and Ca = 0.1. www.nature.com/scientificreports/ where ρ 0 = ρ 2 + h 2 . Profile of this inlet normal velocity, is plotted in Fig. 4. As it is reflected from the this result, the inlet velocity points toward −ẑ direction, that tends to push and repel the particle from wall. Now, the Lorentz reciprocal relation can help us to find the first order velocity and stress tensor that are denoted by u s 1 and T s 1 , respectively. Let us consider a complementary problem with velocity and stress given by v and , respectively. The Lorentz relation reads as:

Scientific Reports
where, short-hand notation for matrix product is used as: a · Ŵ · b = ij a i Ŵ ij b j . Note that, as before, primed variables are reserved for all fields at the upper fluid (on top of the membrane).
We aim to calculate the velocity field at an arbitrary observation point x 0 , located in the lower fluid. We consider a complementary flow that corresponds to a point force located at the observation point x 0 . This complementary problem is subjected to free-slip velocity at z = 0 and it vanishes at infinity. For this complementary problem, we have ∇ · � = −f 0 δ(x − x 0 ) , where f 0 is a given force strength exerted on the fluid. One should note that our main problem (first order fields) and the complementary problem are subjected to different boundaries as shown in Fig. 5.
To be more precise, the velocity fields of our complementary problem in upper and lower fluid, v and v ′ , satisfy the following conditions: where l can be either x or y. Exact solution to this complementary problem, using the well known image method, is presented in 12 . Note that this velocity profile corresponding to a point force was presented before in Eq. (14).  www.nature.com/scientificreports/ We integrate the Lorentz relation, Eq. (21), over the domain of the lower fluid as shown in Fig. 5a and use the divergence theorem to reach the following relation: where we have used the relevant boundary conditions and also benefited the fact that for our main problem ∇ · T s 1 = 0 (no body force exerted). Following a similar procedure and integrating the Lorentz relation over the domain of upper fluid, we will arrive at: Combining Eqs. (24) and (23), we will have: where u s 1 and T s 1 have been defined in Eq. (13) and = ′ − . It should be noted that t · ·ẑ = 0 then, in the second integral of the above equation, the z component of the velocity, u s 1z will appear, and it is given in Eq. (13). Furthermore, as the normal component of v in the third integral vanishes, then T s 1zt will appear in the equation and it is given in Eq. (13).
Having in hand the exact solution to the complementary flow, partially presented in Eq. (14), we can use the above equation and numerically calculate the velocity profile. Choosing point force in directions parallel and perpendicular to the wall, we can obtain the parallel and perpendicular components of the velocity.
Dealing with the first, second and the third integrals that appeared in the above equation, is straightforward. The main difficulty comes with the fourth integral that needs the distribution of first order tension over the sphere as an input function. In other words, the above equation can be considered as an integral equation for the first order flow pattern. Far from the sphere and as an approximation, we can neglect the effects of this integral and proceed with the first terms.

First order force.
To calculate the first order hydrodynamic force exerted on the sphere, F s 1 , the Lorentz relation can be used. Here we choose a complementary flow that corresponds to the flue due to a spherical particle translating with velocity V at point x p . This complementary problem satisfies the stokes equation ( ∇ · = 0 ) with following boundary conditions: where l can be either x or y. Note that a perturbative solution to this complementary problem (up to leading orders of a/h), is presented in Eq. (14). The left-hand side of the Lorentz relation (Eq. 21) vanishes then, we will have: Integrating this equation over the domain of the lower fluid, will give us: Repeating a similar procedure over the domain of upper fluid and combining the result with the above relation, we will have: www.nature.com/scientificreports/ As n points into the fluid, − A s T.n is the force exerted from the particle to the fluid. As one can see, the first order hydrodynamic force acting on the particle has a component perpendicular to the wall.
This perpendicular force is quadratic in terms of U p . In Fig. 2b, the first order contribution to the force is plotted in terms of distance from wall and for various values of . This repulsion (from wall) force will increase for larger values of . Expanding this force for small a/h, the following result can be presented: In terms of U p , the above quadratic force can be expressed as: −(9/2γ )πaη 2 U 2 p (a/h)Ŵ 0 . In Fig. 6, we summarize our results and compare it with the results of a mobility problem, presented elsewhere 16 . Part (a) shows the results of a mobility problem where, a particle under the action of external force f moves near a soft wall. As a result of interface flexibility, the particle will achieve a parallel velocity v � ∝ f � and a perpendicular velocity v ⊥ ∝ f 2 � . The results of our resistivity steady state problem in presented is Fig. 6b, where for a parallel velocity v , parallel force f � ∝ v � and perpendicular force f ⊥ ∝ v 2 � will act on the particle. The direction of first order flow that is quadratic in terms of U p is shown in this figure. This pumping flow points toward −ẑ direction.

Average pumping flow
After finding the steady state solutions to the problem of a moving sphere near a wall, we can use Eq. (3), and find the time dependent flow pattern as: Due to the quasi-steady approximation that we have used, the time dependence in the results is only due to the time evolution of x p (t) or subsequently U p (t) . For our sinusoidal time evolution of the back and forth motion, the zeroth order result for the velocity profile and the hydrodynamics forces acting on the sphere, those are proportional linearly to the velocity, will vanish on average. In zeroth order approximation the membrane remains flat. The first order corrections take into account the deformation of the membrane and these corrections are quadratic in particle's velocity. Such terms will have finite contribution when we average over time. Denoting by T = 2π/ω , the period of back and forth motion, the average velocity profile in a given point can be calculated as: as we expected, only the first order correction contributes to the average velocity. Figure 7 shows the averaged streamlines far below the position of the moving particle. As one can distinguish two important features of the streamlines include its right-left symmetry and the direction of averaged flow. Right-left symmetry is a fact that we expected from the beginning. Interface is isotropic and back and forth u(x, y, z, t) = u s (x − x p (t), y, z). Our resistivity problem where, a particle is constraint to move parallel to a soft wall with a given velocity v . The hydrodynamic force acting on this particle will have both parallel and perpendicular components proportional to v and v 2 , respectively. In this resistivity problem, to keep a fixed distance from the wall, we should apply an external force that is perpendicular to the wall with f e = −f ⊥ to the particle. This hydrodynamic pushing force is a result of a repelling (pumping) flow that is back scattered from the wall. Pumping flow that is the nonlinear part of the velocity profile (first order flow) is shown in this figure. www.nature.com/scientificreports/ motion of the particle is symmetric as well. There is no reason to break the left-right symmetry. But, the up-down symmetry is broken, evidently from the geometrical constraints. Reflected in Fig. 7, the direction of the averaged flow for = 0 is perpendicular to the interface in a way that it tends to repel the particle from interface. To have an intuition about the mechanism behind this repulsing flow, one can see from Fig. 4, the boundary condition for the first order flow. Actually this is the back flow corresponding to the zeroth order problem that needs to satisfy the boundary conditions correctly. As one can see, the inlet flow is concentrated on top of the particle and it points toward the particle. Propagation of this boundary flow, will repel the particle from boundary. Figure 6b also helps to obtain feeling on the pumping flow. As one can see from this figure, independent of the direction of parallel motion (up or down in the figure), a nonlinear contribution to the perpendicular force will arise that is always repelling. The pumping flow is solely due to this repelling force.
Variation of the averaged pumping velocity is studied in Fig. 8. The pumping velocity is plotted in Fig. 8a in terms of the separation from wall. Increasing , will increase the strength of pumping flow. Fig. 8b shows the role of micro-particle position in the pumping flow. As one can see, micro-particles that are near to the wall can produce more stronger perpendicular flow. The parallel component of the flow just below the particle is plotted in part (c). As expected from the beginning the left-right symmetry is reflected in this figure.

Dimensional analysis
Here we use simple dimensional arguments to obtain the order of magnitude for the fluid velocity created by a moving micro-particle behind the interface. Let us assume that the particle velocity is given by U(t) = U p cos(ωt)x and it is located in a fixed position in a distance h, below the interface. This small sphere can be replaced by a time dependent point force given by f (t) = 6πηaU p cos(ωt)x . We want to estimate the fluid velocity in a position just below the point force, with a distance H ≫ h from the wall. As a result of right-left www.nature.com/scientificreports/ symmetry, the flow averaged over time, will not have any component parallel to the membrane. Perpendicular component of the flow needs more inspection.
To obtain the flow velocity, we follow the same steps as we have defined in the perturbation analysis. In order to estimate the quantities, we can distinguish two different limiting regimes. First regime corresponds to the case where, the gravitational and capillary effects have the same order of magnitudes. In this case α ∼ β and capillary length ℓ c , is a relevant length scale in the system. Second regime corresponds to the case where, gravitational effects are negligible ( α ∼ 0 ). In this regime, the capillary number Ca = β −1 , is relevant dimensionless number that takes into account the deformability of the interface.
We begin by the first regime and estimate the instantaneous deformation of the membrane. We denote by and δ , the lateral and normal deformations of the interface, respectively. Looking to the left hand side of Eq. (11), we see that the capillary length is the natural length scale for the lateral deformation, � ∼ ℓ c . Using the same equation and Eq. (16), and replacing by , the variable r that appears at the right hand side, we see that: Our previous perturbative formalism, showed that the first order velocity field is subjected to inlet velocity field given on the undeformed membrane. The inlet profile on the membrane can be estimated from Eq. (13), with a typical term given by u 1z (z = 0) ∼ u 0 ∂ x ξ 1 : To calculate the first order velocity field, we use a typical term from Eq. (23) (Lorentz integral), as: where f 0 is a point force in the observation point and the stress tensor corresponding to this point force scales as � ∼ (f 0 /H 2 ) . Estimating the domain of integration by , this will give us the following average velocity: where ψ is a dimensionless function with ψ(0, 0, 0) ∼ 1 . As it is reflected from above relation, the velocity is nonlinear in terms of f and it averages to a finite value with �f (t) 2 � � = 0.
For the second regime with α = 0 , we can see that h is a natural length scale that limits the lateral deformation so that ∼ h . Using Eq. (11), we see that δ ∼ Ca × a . Now following the same procedure as we used for the first regime, we will arrive at the following result: where φ is a dimensionless function with φ(0, 0) ∼ 1 . Inset in Fig. 8a www.nature.com/scientificreports/ expected to see in the parallel direction. But we expect to obtain a nonlinear contribution to the perpendicular force that can mediate an average flow in that direction.