# Hydrodynamic interaction between particles near elastic interfaces

###### Abstract

We present an analytical calculation of the hydrodynamic interaction between two spherical particles near an elastic interface such as a cell membrane. The theory predicts the frequency dependent self- and pair-mobilities accounting for the finite particle size up to the 5th order in the ratio between particle diameter and wall distance as well as between diameter and interparticle distance. We find that particle motion towards a membrane with pure bending resistance always leads to mutual repulsion similar as in the well-known case of a hard-wall. In the vicinity of a membrane with shearing resistance, however, we observe an attractive interaction in a certain parameter range which is in contrast to the behavior near a hard wall. This attraction might facilitate surface chemical reactions. Furthermore, we show that there exists a frequency range in which the pair-mobility for perpendicular motion exceeds its bulk value, leading to short-lived superdiffusive behavior. Using the analytical particle mobilities we compute collective and relative diffusion coefficients. The appropriateness of the approximations in our analytical results is demonstrated by corresponding boundary integral simulations which are in excellent agreement with the theoretical predictions.

^{†}

^{†}preprint: AIP/123-QED

## I Introduction

The hydrodynamic interaction between particles moving through a liquid is essential to determine the behavior of colloidal suspensions Morrone, Li, and Berne (2012), polymer solutions Usta, Ladd, and Butler (2005); Wojciechowski, Szymczak, and Cieplak (2010), chemical reaction kinetics von Hansen, Netz, and Hinczewski (2010); Długosz et al. (2012), bilayer assembly Ando and Skolnick (2013) or cellular flows Popel and Johnson (2005); Misbah and Wagner (2013). As an example, hydrodynamic interactions result in a notable alteration of the collective motion behavior of catalytically powered self-propelled particles Molina, Nakayama, and Yamamoto (2013) or bacterial suspensions Wensink et al. (2012); Dunkel et al. (2013); Lopez and Lauga (2014); Zöttl and Stark (2014); Elgeti, Winkler, and Gompper (2015). Many of the occurring phenomena can be explained on the basis of two-particle interactions Guazzelli and Morris (2012) which in bulk are well understood. Some of the most intriguing observations, however, are made when particles interact hydrodynamically in the close vicinity of interfaces – a prominent example being the attraction of like-charged colloid particles during their motion away from a hard wall Larsen, Grier et al. (1997); Squires and Brenner (2000); Behrens and Grier (2001).

In the low Reynolds number regime hydrodynamic interactions between two particles are fully described by the mobility tensor which provides a linear relation between the force applied on one particle and the resulting velocity of either the same or the neighboring particle. In an unbounded flow, algebraic expressions for the hydrodynamic interactions between two Felderhof (1977); Kim and Mifflin (1985); Yoon and Kim (1987); Cichocki, Felderhof, and Schmitz (1988); Happel and Brenner (2012); Guazzelli and Morris (2012) and several Deutch and Oppenheim (1971); Batchelor (1976); Ermak and McCammon (1978); Ladd (1988); Ekiel-Jeżewska and Felderhof (2015); Zia, Swan, and Su (2015) spherical particles are well established. Experimentally, the predicted hydrodynamic coupling has been confirmed using optical tweezers Crocker (1997); Meiners and Quake (1999); Bartlett, Henderson, and Mitchell (2001); Henderson, Mitchell, and Bartlett (2002) and atomic force microscopy Radiom et al. (2015).

The presence of an interface is known to drastically alter the hydrodynamic mobility. For a single particle, this wall-induced drag effect has been studied extensively over recent decades theoretically and numerically near a rigid Lauga and Squires (2005); Lorentz (1907); MacKay and Mason (1961); Gotoh and Kaneda (1982); Cichocki and Jones (1998); Franosch and Jeney (2009); Felderhof (2012); Padding and Briels (2010); De Corato et al. (2015); Huang and Szlufarska (2015), a fluid-fluid Lee, Chadwick, and Leal (1979); Berdan and Leal (1982); Bickel (2006, 2007); Bławzdziewicz, Ekiel-Jeżewska, and Wajnryb (2010a, b) or an elastic interface Felderhof (2006); Shlomovitz et al. (2013, 2014); Salez and Mahadevan (2015); Daddi-Moussa-Ider, Guckenberger, and Gekle (2016a, b); Saintyves et al. (2016). While rigid interfaces in general simply lead to a reduction of particle mobility, the memory effect caused by elastic interfaces leads to a frequency dependence of the particle mobility and can cause novel phenomena such as transient subdiffusion Daddi-Moussa-Ider, Guckenberger, and Gekle (2016a). On the experimental side, the single particle mobility has been investigated using optical tweezers Faucheux and Libchaber (1994); Dufresne, Altman, and Grier (2001); Schäffer, Nørrelykke, and Howard (2007), evanescent wave dynamic light scattering Holmqvist, Dhont, and Lang (2007); Michailidou et al. (2009); Wang, Prabhakar, and Sevick (2009); Kazoe and Yoda (2011); Lisicki et al. (2012); Rogers et al. (2012); Michailidou et al. (2013); Wang and Huang (2014); Watarai and Iwai (2014); Lisicki et al. (2014) or video microscopy Eral et al. (2010); Cervantes-Martínez et al. (2011); Dettmer et al. (2014). The influence of a nearby elastic cell membrane has recently been investigated using magnetic particle actuation Irmscher et al. (2012) and optical traps Shlomovitz et al. (2013); Boatwright et al. (2014); Jünger et al. (2015).

Hydrodynamic interactions between two particles near a planar rigid wall have been studied theoretically Squires and Brenner (2000); Swan and Brady (2007) and experimentally using optical tweezers Lele et al. (2011); Tränkle, Ruh, and Rohrbach (2016) and digital video microscopy Dufresne et al. (2000). Narrow channels Cui, Diamant, and Lin (2002); Misiunas et al. (2015), 2D confinement Bleibel et al. (2014) or liquid-liquid interfaces have also been investigated Zhang et al. (2013, 2014). Near elastic interfaces, however, no work regarding hydrodynamic interactions has so far been reported. Given the complex behavior of a single particle near an elastic interface (caused by the above-mentioned memory effect) such hydrodynamic interactions can be expected to present a very rich phenomenology.

In this paper, we calculate the motion of two spherical particles positioned above an elastic membrane both analytically and numerically. We find that the shearing and bending related parts in the pair-mobility can in some situations have opposite contributions to the total mobility. Most prominently, we find that two particles approaching an idealized membrane exhibiting only shear resistance will be attracted to each other which is just opposite to the well-known hydrodynamic repulsion for motion towards a hard wall Squires and Brenner (2000). Additionally, we show that the pair-mobility at intermediate frequencies may even exceed its bulk value, a feature which is not observed in bulk or near a rigid wall. This increase in pair-mobility results in a short-lived superdiffusion in the joint mean-square displacement.

The remainder of the paper is organized as follows. In Sec. II, we introduce the theoretical approach to computing the frequency-dependent self- and pair-mobilities from the multipole expansion and Faxén’s theorem, up to the 5th order in the ratio between particle radius and particle-wall or particle-particle distance. In Sec. III, we present the boundary integral method which we have used to numerically confirm our theoretical predictions. In Sec. IV, we provide analytical expressions of the particle self- and pair-mobilities in terms of power series, finding excellent agreement with our numerical simulations. Expressions of the self- and pair-diffusion coefficients are derived in Sec. V. Concluding remarks are made in Sec. VI.

## Ii Theory

We consider a pair of particles of radius suspended in a Newtonian fluid of viscosity above a planar elastic membrane extending in the plane. The two particles are placed at and , i.e. the line connecting the two particles is parallel to the undisplaced membrane. We denote by the center-to-center separation measured from the left () to the right () particle (see Fig. 1 for an illustration).

The particle mobility is a tensorial quantity that linearly couples the velocity of particle in direction to an external force in the direction applied on the same () or the other () particle. Transforming to the frequency domain we thus have (Kim and Karrila, 2005, ch. 7)

where Einstein’s convention for summation over repeated indices is assumed. The particle mobility tensor in the present geometry can be written as an algebraic sum of two distinct contributions

(1) |

where is the pair-mobility in an unbounded geometry (bulk flow), and is the frequency-dependent correction due to the presence of the elastic membrane. An analogous relation holds for .

For the determination of the particle mobility, we consider a force density acting on the surface of the particle , related to the total force by

which induces the disturbance flow velocity at point

(2) |

where denotes the velocity Green’s function (Stokeslet), i.e. the flow velocity field resulting from a point-force acting on . The disturbance velocity at any point can be split up into two parts,

(3) |

where is the flow field induced by the particle in an unbounded geometry, and is the flow satisfying the no-slip boundary condition at the membrane. In this way, the Green’s function can be written as

(4) |

where is the infinite-space Green’s function (Oseen’s tensor) given by

(5) |

with and . The term represents the frequency-dependent correction due to the presence of the membrane. Far away from the particle , the vector in Eq. (2) can be expanded around the particle center following a multipole expansion approach. Up to the second order, and assuming a constant force density, the disturbance velocity can be approximated by Kim and Netz (2006); Gauger, Downton, and Stark (2008); Swan and Brady (2007, 2010)

(6) |

where stands for the gradient operator taken with respect to the singularity position . Note that for a single sphere in bulk, the flow field given by Eq. (6) satisfies exactly the no-slip boundary conditions at the surface of the sphere, i.e. in the frame moving with the particle, both the normal and tangential velocities vanish. Using Faxén’s theorem, the velocity of the second particle in this flow reads Kim and Netz (2006); Gauger, Downton, and Stark (2008); Swan and Brady (2007, 2010)

(7) |

where denotes the usual bulk mobility, given by the Stokes’ law. The disturbance flow incorporates both the disturbance from the particle and the disturbance caused by the presence of the membrane. By plugging Eq. (6) into Faxén’s formula given by Eq. (7), the component of the frequency-dependent pair-mobilities can be obtained from

(8) |

For the self-mobilities, only the correction in the flow field due to the presence of the membrane in Eq. (3) is considered in Faxén’s formula (the influence of the second particle on the self-mobility is neglected here for simplicity Swan and Brady (2007); Gauger, Downton, and Stark (2008)). Therefore, the frequency-dependent self-mobilities read

(9) |

and analogously for .

In order to use the particle pair- and self-mobilities from Eqs. (8) and (9), the velocity Green’s functions in the presence of the membrane are required. These have been calculated in our earlier work Daddi-Moussa-Ider, Guckenberger, and Gekle (2016a) and their derivation is only briefly sketched here with more details in Appendix A.

We proceed by solving the steady Stokes equations with an arbitrary time-dependent point-force acting at ,

(10) | ||||

(11) |

where is the pressure field. The determination of the Green’s functions at is straightforward thanks to the system translational symmetry along the plane. After solving the above equations and appropriately applying the boundary conditions at the membrane, we find that the Green’s functions are conveniently expressed by

(12a) | ||||

(12b) | ||||

(12c) | ||||

(12d) |

where , with . Here denotes the Bessel function of the first kind of order . The functions , and are provided in Appendix A. It is worth to mention here that the unsteady term in the Stokes equations leads to negligible contribution in the correction to the Green’s functions Daddi-Moussa-Ider, Guckenberger, and Gekle (2016a), and it is therefore not considered in the present work.

The membrane elasticity is described by the well-established Skalak model Skalak et al. (1973), commonly used to describe deformation properties of red blood cell (RBC) membranes Eggleton and Popel (1998); Krüger, Varnik, and Raabe (2011); Krüger (2012). The elastic model has as parameters the shearing modulus and the area-expansion modulus . The two moduli are related via the dimensionless number . Moreover, the membrane resists towards bending according to Helfrich’s model Helfrich (1973), with the corresponding bending rigidity .

## Iii Boundary integral methods

In this section, we introduce the numerical method used to compute the particle self- and pair- mobilities. The numerical results will subsequently be compared with the analytical predictions presented in Sec. II.

For solving the fluid motion equations in the inertia-free Stokes regime, we use a boundary integral method (BIM). The method is well suited for problems with deforming boundaries such as RBC membranes Pozrikidis (2001); Zhao and Shaqfeh (2011). In order to solve for the particle velocity given an exerted force, a completed double layer boundary integral method (CDLBIM) Power and Miranda (1987); Kohr and Pop (2004) has been combined with the classical BIM Zhao, Shaqfeh, and Narsimhan (2012). The integral equations for the two-particle membrane systems read

(13) |

where is the surface of the elastic membrane and is the surface of the two spheres. Here denotes the velocity on the membrane whereas represents the double layer density function on , related to the velocity of the particle via

(14) |

where are known functions Kohr and Pop (2004). The brackets stand for the inner product in the space of real functions whose domain is , and the function is defined by

with being the centroid of the sphere labeled upon which the force is applied. The single layer integral is defined as

and the double layer integral as

Here, is the traction jump, denotes the outer normal vector at the particle surfaces and is the force acting on the rigid particle. The infinite-space Green’s function is given by Eq. (5) and the corresponding Stresslet, defined as the symmetric part of the first moment of the force density, reads Kim and Karrila (2005)

with and . The traction jump across the membrane is an input for the equations, determined from the instantaneous deformation of the membrane. In order to solve Eqs. (13) numerically, the membrane and particles’ surfaces are discretized with flat triangles. The resulting linear system of equations for the velocity on the membrane and the density on the rigid particles is solved iteratively by GMRES Saad and Schultz (1986). The velocity of each particle is determined from (14). For further details concerning the algorithm and its implementation, we refer the reader to Ref. Daddi-Moussa-Ider, Guckenberger, and Gekle (2016a). Bending forces are computed using Method C from Guckenberger et al. (2016).

In order to compute the particle self- and pair-mobilities numerically, a harmonic oscillating force of amplitude and frequency is applied at the surface of the particle . After a brief transient time, both particles begin to oscillate at the same frequency as and as . The velocity amplitudes and phase shifts can accurately be obtained by a fitting procedure of the numerically recorded particle velocities. For that, we use a nonlinear least-squares algorithm based on the trust region method Conn, Gould, and Toint (2000). Afterward, the component of the frequency-dependent complex self- and pair-mobilities can be calculated as

## Iv Results

For a single membrane, the corrections to the particle mobility can conveniently be split up into a correction due to shearing and area expansion together with a correction due to bending Daddi-Moussa-Ider, Guckenberger, and Gekle (2016a). In the following, we denote by (“self”) the components of the self-mobility tensor, and by (“pair”) the components of the pair-mobility tensor. Note that for , and that .

### iv.1 Self-mobilities for finite-sized particles

Mathematical expressions for the translational particle self-mobility corrections will be derived in terms of . The point-particle approximation presented in earlier work Daddi-Moussa-Ider, Guckenberger, and Gekle (2016a) represents the first order in the perturbation series, valid when the particle is far away from the membrane.

#### iv.1.1 Perpendicular to membrane

The particle mobility perpendicular to the membrane is readily obtained after plugging the correction as defined by Eq. (4) to the normal-normal component of the Green’s function from Eq. (12a) into Eq. (9). After computation, we find that the contribution due to shearing and bending can be expressed as

(15a) | ||||

(15b) |

where the subscripts S and B stand for shearing and bending, respectively. The function is the generalized exponential integral defined as Abramowitz, Stegun et al. (1972). Furthermore, is a dimensionless frequency associated with the shearing resistance, whereas . Moreover, is a dimensionless number associated with bending. The functions , with are defined by

with

where and being the principal cubic-root of unity. The bar designates complex conjugate.

The total mobility correction is obtained by adding the individual contributions due to shearing and bending, as given by Eqs. (15a) and (15b). In the vanishing frequency limit, the known result for a hard-wall Swan and Brady (2010) is obtained:

(16) |

The particle mobility near an elastic membrane is determined by membrane shearing and bending properties.
We therefore consider a typical case for which both effect manifests themselves equally.
For that purpose, we define a characteristic time scale for shearing as together with a characteristic time scale for bending as Daddi-Moussa-Ider, Guckenberger, and Gekle (2016a). Then we take such that the two time scales are equal and can be denoted by .
In this case, the two dimensionless numbers and are related by .
The situation for a membrane with the typical parameters of a red blood cell is qualitatively similar as shown in the Supporting Information.
^{1}^{1}1See Supplemental Material at [URL will be inserted by publisher] for the frequency-dependent mobilities where typical values for the RBC parameters are used.

In Fig. 2 , we show the particle scaled self-mobility corrections versus the scaled frequency , as stated by Eqs. (15a) and (15b). The particle is set a distance above the membrane. We observe that the real part is a monotonically increasing function with respect to frequency while the imaginary part exhibits a bell-shaped dependence on frequency centered around . In the limit of infinite frequencies, both the real and imaginary parts of the self-mobility corrections vanish, and thus one recovers the bulk behavior. For the perpendicular motion we observe that the particle mobility correction is primarily determined by the bending part.

A very good agreement is obtained between the analytical predictions and the numerical simulations over the whole range of frequencies. Additionally, we assess the accuracy of the point-particle approximation employed in earlier work Daddi-Moussa-Ider, Guckenberger, and Gekle (2016a), in which only the first order correction term in the perturbation parameter was considered. While this approximation slightly underestimates particle mobilities, it nevertheless leads to a surprisingly good prediction, even though the particle is set only one diameter above the membrane.

#### iv.1.2 Parallel to membrane

We proceed in a similar way for the motion parallel to the membrane. By plugging the correction from the Green’s function in Eq. (12b) into Eq. (9) we find

(17a) | ||||

(17b) |

where we defined

The well-known hard-wall limit, as first calculated by FaxénFaxén (1922); Swan and Brady (2010), is recovered by considering the vanishing frequency limit:

(18) |

The mobility corrections in the parallel direction are shown in Fig. 2 . We observe that the total correction is mainly determined by the shearing part in contrast to the perpendicular case where bending dominates.

### iv.2 Pair-mobilities for finite-sized particles

In the following, expressions for the pair-mobility corrections in terms of a power series in will be provided. To start, let us first recall the particle pair-mobilities in an unbounded geometry. By applying Eq. (8) to the infinite space Green’s function Eq. (5), the bulk pair-mobilities for the motion perpendicular to and along the line of centers read (Kim and Karrila, 2005, p. 190)

(19) |

and are commonly denominated the Rotne-Prager tensor Rotne and Prager (1969); Ermak and McCammon (1978). Note that the terms with vanish for the bulk mobilities when considering only the first reflection as is done here. The axial symmetry along the line connecting the two spheres in bulk requires that and that the off-diagonal components of the mobility tensor are zero. Physically, the parameter only takes values between 0 and as overlap between the two particles is not allowed. In this interval, the pair-mobility perpendicular to the line of centers is always lower than the pair-mobility , since it is easier to move the fluid aside than to push it into or to squeeze it out of the gap between the two particles.

Consider next the pair-mobilities near an elastic membrane. By applying Eq. (8) to Eqs. (12a) through (12d), we find that the corrections to the pair-mobilities can conveniently be expressed in terms of the following convergent integrals,

(20a) | ||||

(20b) | ||||

(20c) | ||||

(20d) |

where and

The terms involving and in Eqs. (20a) through (20d) are the contributions coming from shearing and bending, respectively. Due to symmetry, for .

For future reference, we note that each component of the frequency-dependent particle self- and pair-mobility tensor can conveniently be cast in the form

(21) |

where indices and superscripts have been omitted. Here denotes the scaled bulk mobility (cf. Eq. (1)), and the integral term represents either shearing or bending related parts in the mobility correction. Note that and are real functions which do not depend on frequency. Moreover, or for the shearing related parts and for bending such that , .

In the vanishing frequency limit, i.e. for both taken to zero we recover the pair-mobilities near a hard-wall with stick boundary conditions, namely

(22a) | ||||

(22b) | ||||

(22c) | ||||

(22d) |

in agreement with the results by Swan and Brady Swan and Brady (2007).

In Fig. 3 we plot the particle pair-mobilities as given by Eqs. (20a) through (20d) as functions of the dimensionless frequency for . We observe that the real and imaginary parts have basically the same evolution as the self-mobilities. Nevertheless, two qualitatively different effects are apparent from Fig. 3: First, the amplitude of the normal-normal pair-mobility in a small frequency range even exceeds its bulk value. This enhanced mobility results in a short-lasting superdiffusive behavior as will be described in Sec. V.

Secondly, for the components and in Fig. 3 we find that, unlike the self-mobilities, shearing and bending may have opposite contributions to the total pair-mobilities. For the component this implies the interesting behavior that hydrodynamic interactions can be either attractive or repulsive depending on the membrane properties. This will be investigated in more detail in the next subsection.

### iv.3 Perpendicular steady motion

A situation in which hydrodynamic interactions are particularly relevant is the steady approach of two particles towards an interface, such as e.g. drug molecules approaching a cell membrane, reactant species approaching a catalyst interface, charged colloids being attracted by an oppositely charged membrane, etc. For hard walls, it is known that hydrodynamic interactions in this case are repulsive Dufresne et al. (2000); Squires and Brenner (2000); Swan and Brady (2007) leading to the dispersion of particles on the surface. Near elastic membranes, the different signs of the bending and shear contributions to the pair-mobility in Fig. 3 point to a much more complex scenario including the possibility of particle attraction.

The physical situation of two particles being initially located at and suddenly set into motion towards the interface is described by a Heaviside step function force . Its Fourier transform to the frequency domain reads Bracewell (1999)

Using the general form of Eq. (21), the scaled particle velocity in the temporal domain is then given by

(23) |

where is a dimensionless time. At larger times, the exponential in Eq. (23) can be neglected compared to one. In this way, we recover the steady velocity near a hard-wall.

In corresponding BIM simulations, a constant force of small amplitude towards the wall is applied on both particles in order to retain the system symmetry. At the end of the simulations, the vertical position of the particles changes by about 8 % compared to their initial positions .

In Fig. 4 we show the time dependence of the vertical velocity which at first increases and then approaches its steady-state value. Figure 4 shows the relative velocity between the two particles: clearly, the motion is attractive for a membrane with negligible bending resistance (such as a typical artificial capsule) which is the opposite of the behavior near a membrane with only bending resistance (such as a vesicle) or a hard wall.

In order to illustrate more clearly for which wall and particle distances a repulsion/attraction is expected we show in Fig. 5 the pair-mobility correction for the shear and bending contributions in the plane. To reduce the parameter space and to bring out the considered effects most clearly, we consider the idealized limit . In this limit, the contributions become independent of the elastic moduli since directly implies that meaning that even infinitesimally small shearing and bending resistances would make the membrane behave identical to the hard wall. This unphysical behavior is remedied in a realistic situation where a small bending resistance will lead to a correspondingly large time scale and thus to a long-lived transient regime as given by Eq. (23) and shown in the Supporting Information. Therefore, the contours shown in Fig. 5 faithfully represent the behavior of membranes with small bending (Fig. 5 ) or small shear (Fig. 5 ) resistance. The corresponding equations can be found in Appendix B.

By equating Eq. (49d) to zero and solving the resulting equation perturbatively, the threshold lines where the shearing contribution changes sign are given up to fifth order in by

(24) |

Eq. (24) is shown as circles in Fig. 5. The bending contribution in Fig. 5 always has a positive sign corresponding to a repulsive interaction similar as the hard wall.

Similar changes in sign are observed for for shear and for bending. The corresponding contours are given in the Supporting Information. Their physical relevance, however, is less important than for shown in Fig. 5 as the effects may be overshadowed by the bulk values of the mobilities (which is zero only for ).

## V Diffusion

The diffusive dynamics of a pair of Brownian particles is governed by the generalized Langevin equation written for each velocity component of particle as Kubo (1966)

(25) |

A similar equation can be written for the velocity components of the other particle . Here, denotes the particles’ mass, stands for the time-dependent two-particle friction retardation tensor (expressed in kg/s) and is a random force which is zero on average. By evaluating the Fourier transform of both members in Eq. (25) and using the change of variables together with the shift property in the time domain of Fourier transforms we get

(26) |

where and are the Fourier-Laplace transforms of the retardation function defined as

(27) |

and analogously for .

In the following, we shall consider the overdamped regime for which the particles are massless . Solving Eq. (26) for the particle velocities and equating with the definition of the mobilities,

(28a) | ||||

(28b) |

leads to expressions of the mobilities in terms of the friction coefficients:

where the brackets [ ] are dropped out for the sake of clarity. Similar as for the mobilities, the self- and pair components of the retardation function are denoted by and , respectively. Note that so that as required by symmetry.